© ESO, 2012

1. Introduction

One of the most complicated phases of the evolution of low- and intermediate- mass stars is the thermally pulsing-asymptotic giant branch (TP-AGB) stage. The physical conditions that are sampled by models of this phase are quite extreme and the recurring thermal instability of the helium-burning shell requires codes to be robust on short timescales during the thermal pulse cycle. It is common for convergence problems to arise during the study of such stars, although rarely are these problems discussed in the recent literature. The failure to converge usually occurs during the later thermal pulses, when the envelope mass is relatively low due to mass loss through the AGB phase. We are currently using the Monash stellar evolution code monstar (Campbell & Lattanzio 2008) to study the evolution of super asymptotic giant branch (super-AGB) stars and have encountered such a problem.

When the envelope mass is below ≈ 2 M_⊙, the code’s solution typically tends to a configuration with a negative gas pressure at the bottom of the convective envelope. We find that this convergence problem can be delayed by slowly increasing α, the ratio of the convective mixing length to the pressure scale height. Our models are originally evolved with α = 1.75. Typically, with a slightly higher α, a few more thermal pulses could be modelled until the same convergence problem occurs again. This works until the envelope mass is in the region of 0.5−1 M_⊙, when increasing α will no longer avoid the problem. A typical value of α ~ 4 is needed to evolve the model up to this stage. Previously, Herwig (2001) has also suggested that a higher α value (of 3) helps to ease post-AGB numerical simulations. Recipes of adopting higher values of α during the later stages of the AGB evolution are also mentioned in Miller Bertolami & Althaus (2006), Kitsikis (2008) and Weiss & Ferguson (2009).

In this paper we seek a physical explanation for such behaviour and find that the radiation pressure is so large that it supplies all the pressure required for hydrostatic equilibrium. In other words, the local luminosity exceeds the Eddington limit. This type of instability has previously been discussed by Wood & Faulkner (1986) who suggested that the envelope could be ejected and the star would then evolve to the planetary nebula. Eddington limit luminosities as the cause of the numerical problems during low-mass envelope thermal pulses has already been identified by Lawlor & MacDonald (2003)¹ and Miller Bertolami et al. (2006).

The stability of stellar models during the late thermally pulsing AGB phase has also been considered by Han et al. (1994), Wagenhuber & Weiss (1994), Wood & Vassiliadis (1993), Renzini et al. (1992), Tuchman (1984), Barkat & Tuchman (1980), Tuchman et al. (1978), Wood (1973), and Paczyński & Ziółkowski (1968). The transition from AGB to post-AGB evolution has been considered, from the observational point of view, by Engels et al. (2009).

In this paper, we discuss in detail the cause of the instability in the models and whether real stars do encounter such an instability in nature. We also discuss possible scenarios for when this instability may occur and try to determine whether the envelope would be ejected based on the results of our models. The next section of the present work gives a physical description of the instability, considering the cases when energy transport is either dominated by radiation or cases when convection is efficient. In the third section a hypothesis for the cause of the instability based on the existence of an iron element opacity peak is proposed and tested. In the fourth section we consider the possible fate of the star after the instability. Finally we present the main conclusions and discuss our work.

2. Physical description of the instability

The computation of the late stages of massive TP-(super)AGBs is hampered by the occurrence of a type of instability that causes evolutionary models to fail to converge and can affect our understanding of the most advanced phases of the evolution and eventually the resulting white dwarf remnants. This instability is characterised by the development of a zone at the base of the convective envelope in which the density and the gas pressure, P_gas tend to zero and, therefore, only radiation contributes to the total pressure in the equation of state. Eventually the models converge to solutions with local P_gas below zero near the base of the convective envelope (BCE) and the codes finally crash. As we will show in the next subsection, this is equivalent to having local luminosity values higher than the Eddington luminosity. The temperature in the stellar zone in which the instability develops is typically 1 − 3 × 10⁵ K and thus much higher than the values found in the region where hydrogen and helium are partially ionized. Hence the instability we are concerned with is not directly related to the instability described by Wagenhuber & Weiss (1994). These authors identified the local decrease in efficiency of the energy transport due to the recombination of hydrogen as the main cause of the instability that quenched the computation of AGB evolution.

We have noticed the presence of the instability described above when using the Mount Stromlo and Monash Stellar Evolution code (MONSTAR), but from Siess (priv. comm.) we know that the stellar evolution code STAREVOL (Siess 2010), presents the same convergence problems during late massive AGB and super-AGB evolution. Both codes and their respective super-AGB models have been described by Doherty et al. (2010). It is worth noting that such convergence problems are common. Many evolution codes experience numerical problems that could be related to this instability. The MONSTAR code (Campbell & Lattanzio 2008) and the STAREVOL code (Siess 2007, 2010) and the STARS code (Stancliffe & Jeffery 2007; Lau et al. 2009) have convergence problems during the late massive AGB or super-AGB evolution.

We have also checked that whether this instability occurs is insensitive to most input physics of the code. For example, we have used different implementations of the OPAL opacities, of the treatment of mixing (with or without time independent mixing), or of the mass-loss prescription: Reimers (1978), Vassiliadis & Wood (1993), or Bloecker (1995b). The models presented here are modelled with the Vassiliadis & Wood (1993) mass-loss prescription. These make only minor differences to the occurrence of the instability. Unless one employs a very high mass-loss rate, such as the one used by Kovetz et al. (2009), the instability cannot be avoided². In fact different input physics does change the envelope mass when the instability begins, but they can not prevent it from occuring. The opacity calculation could on the other hand, have a huge effect on whether this instability occurs, as shown in a later section. By removing the Fe opacity peak, it is possible to avoid the instability. Computations with old Los Alamos(LAOL) opacities may not suffer this instability e.g. Bloecker (1995a). In fact, the main improvement of the new opacities (OPAL/OP) tables is the treatment of Fe and iron-group opacities in the range of temperatures 10⁶ − 10⁵ K (Iglesias & Rogers 1996a). This could explain why instabilities are more common with the adoption of new opacities. This is consistent with our identification of the iron opacity peak as the cause of the instability in the models discussed in this paper.

Fig. 1 Density, temperature, radius, luminosity, κ and β, the ratio of gas pressure to total pressure in the instability region. The model profiles represented in solid lines correspond to the times labelled in Fig. 3 and in particular thick lines represent the last converged models. Dashed lines represent the final progress of the instability.

In terms of numerical details, we have tried varying the obvious things, such as increasing the temporal or spatial resolution, but these make little difference. We may delay the instability for a few pulses but not eliminate it. However, we discovered that increasing the mixing length parameter α enables the code to continue to converge for a little longer. This was also discovered by Herwig (2001), Miller Bertolami & Althaus (2006), Kitsikis (2008) and Weiss & Ferguson (2009). In this way we can push the evolution for a few more thermal pulses, although inevitably the convergence problem re-appears. A further increase in α again delays the problem. This is a hint to the underlying physics but we are not free to increase α arbitrarily and without limit!

The instability is manifested in the mesh points at the base of the convective envelope. Physical conditions at this time are shown in Figs. 1 to 5 for a 8.5 M_⊙ star of solar metallicity. The MONSTAR code uses temperature and pressure as dependent variables and the iterations soon produce a T and P corresponding to such a high radiation pressure, P_rad that the implied gas pressure P_gas = P − P_r is negative, as in indicated in Fig. 1. From Fig. 4 we can see that the instability starts at the boundary between radiative and convective zones. As the instability develops, the density at the bottom of the convective envelope decreases (Fig. 1) because the model star is expanding, as evidenced by the rapid increase in radius around that region (Fig. 1).

Fig. 2 Upper panel: evolution of luminosity associated to H-burning (dotted line) and He-burning (solid line); middle panel: evolution of convection; lower panel: evolution of the surface radius, all during the second last thermal pulse. Relevant structure profiles of models labelled a to f appear in Fig. 8.

Fig. 3 Upper panel: evolution of luminosity associated to H-burning (dotted line) and He-burning (solid line); middle panel: evolution of convection; lower panel: evolution of the surface radius, all during the last thermal pulse. Relevant structure profiles of models labelled a to f appear in Fig. 9.

Wood & Faulkner (1986), section IIIb, while developing an analysis of hydrostatic sequences of planetary nebulae, justified the existence of the TP-AGB instability. They computed models that would produce the nuclei of the PNe between 0.60 and 0.89 M_⊙, and realised that, at some point in the late evolution of TP-AGB stars, the radiation pressure, P_rad, becomes so large that its contribution accounts practically for the total pressure at the base of the convective envelope. So P_gas is forced to be zero. These authors suggest that, when this phenomenon takes place in a real star, a hydrodynamic expansion of the stellar envelope occurs and that convergence problems develop in the models as a consequence of the lack of existence of a hydrostatic solution.

This is consistent with the instability we are studying in this work. In particular, β, the ratio of gas pressure to total pressure at the base of envelope, drops as the stars evolve along the last thermal pulses, with our last successfully converged models yielding values of β below 0.001. The behaviour of P_gas tending to (or falling below) zero is accompanied and, in fact, is equivalent to local luminosity values reaching (or exceeding) the Eddington luminosity. We will develop this further in the next subsections.

2.1. The importance of the convective efficiency

The case where the energy is carried by radiation alone has already been presented in Wood & Faulkner (1986). The situation is more complicated when convection is present because it adds a second energy transport mechanism. The case of vanishingly low convective efficiency corresponds to the radiative case. We are interested in how convective efficiency modifies the behaviour of the star.

In a region where convection is efficient, the local luminosity can exceed the Eddington limit without encountering the above instability. This occurs frequently during the lives of stars. The core helium flash is one example, where the luminosity can be 10⁸ L_⊙ or more, greatly exceeding the local Eddington value. Instead convection sets in when the local luminosity exceeds the Eddington limit and the highly efficient convection carries most of the energy and the instability does not develop. Similarly, during AGB thermal pulses themselves, in the helium burning region we have similarly high values and again the efficient convection helps to keep β positive even though the local luminosity exceeds the Eddington value.

In a region where convection is efficient, most of the energy is transported by convection instead of radiation. We derive below a corresponding Eddington Limit for a convective region and show that the value depends on the efficiency of convection. In regions where convection is inefficient the local Eddington limit still applies because most energy is transported by radiation.

Fig. 4 Relevant structure and composition parameters for the model labeled c in Fig. 3. Upper panel: actual gradient (thin solid), radiative gradient (dotted line) and β (thick solid). Middle panel: local luminosity values (dotted line), Eddington luminosity (thin solid) and Eddington luminosity corrected with the effect of convection (thick solid line). Lower panel: hydrogen (solid line), helium (dotted line) and carbon (dashed line) composition profiles.

Fig. 5 Relevant structure and composition parameters for the model labeled f in Fig. 3. Upper panel: actual gradient (thin solid), radiative gradient (dotted line) and β (thick solid). Middle panel: local luminosity values (dotted line), Eddington luminosity (thin solid) and Eddington luminosity corrected with the effect of convection (thick solid line). Lower panel: hydrogen (solid line), helium (dotted line) and carbon (dashed line) composition profiles.

We begin with the equation of pressure equilibrium. (1)where, P is the pressure, r is the radius, G is the gravitational constant and M_r is the mass within the radius.

The energy transfer equation is given by (2)where T is the temperature, κ is the opacity, ρ is the density and L_r is the local luminosity. In the radiative case, ∇ = ∇_r and ∇_r is defined as .

The above equation can be rewritten in the form (3)where (4)is the gradient if the star is radiative and ∇ is the actual gradient of the star. By definition the radiative pressure , so we can obtain (5)Combining with we obtain the differential equation (6)As we can see from the above equation, a jump in local luminosity or opacity will cause the local radiation pressure to increase much more rapidly than the local pressure. It is therefore possible that the radiation pressure completely dominates the total pressure.

In particular, for the base of the convective envelope (BCE), the stellar region in which we are interested, L_r and M_r can be well approximated by the total luminosity of the star, L, and its core mass M_c, respectively. If we assume the envelope mass is small and take M_r as the core mass M, then the luminosity L is a constant and equals to the value at the base of the convective envelope. By integrating over a narrow region around the BCE and using a local average value for the opacity ⟨ κ ⟩ , we can express the ratio of radiation pressure to total pressure at the bottom of the envelope as follows: (7)where , which is the average of the opacity over the region.

From the above equation, we note that β < 0 is equivalent to , where . In convective region where convection is efficient ∇ is very nearly equal to the adiabatic gradient (∇_a) and is much smaller than ∇_r, and hence the equivalent Eddington limit is much larger and the radiation pressure does not dominate the total pressure. However, in regions where convection is inefficient, the value of ∇ approaches the value of ∇_r, and the above equation approaches the equation , that is, the Eddington limit converges back to the classical value in a radiative region.

The actual value ∇ depends on the efficiency of the convection and must satisfy the relation ∇_r > ∇ > ∇_a, where ∇_a is the adiabatic gradient. But it is important to keep in mind that radiation and convection can coexist in the same region of a star and both can be relatively efficient. It is only in some particular regions, such as the stellar centre, in which convection fully dominates and radiation can be ignored.

Kippenhahn & Weigert (1994, p. 51) derive the value of ∇ based on mixing length theory. The most important quantity is (8)where l_m is the mixing length and H_P is the pressure scale-height. The quantity U can be taken as the ratio of the radiative “conductivity” to the convective “conductivity”. If U is large, convection is ineffective and cannot transport a substantial fraction of the luminosity. From the equation, also in Kippenhahn & Weigert (1994)(9)where ∇_e is related to the temperature changes of the moving elements. We see that, in order for the equation to have solution when U is larger, ∇ ≈ ∇_r and hence most of the energy is transported by radiation.

Because U is inversely proportional to ρ², convection is very efficient in the dense central part of the star but not necessarily in the envelope where the density is low. In our model, U is of the order of 100 at the bottom of the convective envelope, and according to Kippenhahn & Weigert (1994) this value of U implies radiation dominates over convection for the energy transport. Therefore, ∇ ~ ∇_r and the actual Eddington limit is close to the radiative value. When the local luminosity exceeds this the instability occurs.

Notice that this naturally gives us the explanation for why increasing the mixing-length can delay the onset of the convergence problems. We see that U is inversely proportional to . So by increasing α we increase the value of l_m, decreasing the value of U, and forcing the convection to be more efficient. A smaller U implies more energy is transported by convection. Hence ∇ decreases and the effective Eddington luminosity increases beyond the radiative value. Thus by increasing α we enable more efficient convection to carry the energy away and the instability is avoided.

It is interesting to see how the above argument is reflected on an actual computation of a full evolutionary sequence. The next few paragraphs compare our previous analysis with the results of an actual computation of the final stages of a TP-(super)AGB star. We have considered a 8.5 M_⊙ star of solar metallicity, from its main sequence, until the instability is reached after 152 thermal pulses. At this time our model star has 2.71 M_⊙, of which 1.17 M_⊙ correspond to the core and the remaining 1.54 M_⊙ constitute the envelope.

Figure 2 shows the temporal evolution of luminosity, convection and radius during the second last thermal pulse. The upper panel shows the temporal evolution of the luminosity associated to H-burning, L_H, and He-burning, L_He, the middle panel shows the evolution of the BCE and the convective shell associated with the He-flash and the lower panel shows the variation of the surface radius with time. The labelled vertical lines correspond to selected models whose structure profiles are shown in Figs. 8 and 9.

Figure 3 shows the same as the former figure, for the case of the last thermal pulse. As we can see, when L_He decays after the flash L_H increases and tends to recover its interpulse value. But, opposite to what happens in former thermal-pulses, L_H reaches a maximum of about 10^-3 L_⊙ and then decreases again. This is because the stars tries to expand around the region where the gas pressure is approaching zero. As a result, the region around the burning shell cools, and the hydrogen burning shell cannot reignite.

Figure 4 shows a section of the mass profile of the star near the BCE during the last converging thermal pulse (at time labelled c in Fig. 3). Figure 5 shows the same as Fig. 4 for the time labelled f, that is, for our last converged model. The upper panel of both figures represent the actual gradient, ∇ (thin solid line), and the radiative gradient, ∇_r (dotted line). Convection dominates the energy transport in the regions where ∇_r ≫ ∇. The thick solid line represents the value of . As we can see β approaches zero at the bottom of the convective envelope. The middle panels represent the local luminosity, L_r (dashed line), the Eddington luminosity computed for radiative regions (thin solid line) and the Eddington luminosity computed for convective regions (thick solid line), that is, multiplied by the factor , L_Edd,conv – in red. In the convective region, the local luminosity could be higher than the Eddington luminosity, but lower than L_Edd,conv, as we have justified above. However, in the last converged models (see Fig. 5), L_r ~ L_Edd,conv and β ~ 0 in the convective zones and we encounter the convergence problem. Note that the instability is actually occuring in a very narrow region of about 1.5 × 10^-3 M_⊙ near the base of the convective envelope.

Note that both the hydrogen and the helium burning shells are active at the time of the evolution labelled c. But in the very last converged model, labelled f, the hydrogen burning shell is completely switched off. Finally, for a reference, the lower panels represent the abundance profile of H, ⁴He and ¹²C.

In a summary, our analysis of the convergence limits in our 1D-stellar evolutionary code allows us to explain a number of consequences that we actually encounter when we perform full evolutionary computations of super-AGB stars. The code fails to converge when the envelope masses are still relatively large (even above 2 M_⊙). The divergence, associated with the increase of luminosity in the intershell region above L_Edd, occurs in narrow region near the base of the convective envelope. Failure of convergence is temporarily solved by increasing the mixing length parameter through an increase in α. This increase ultimately allows an increase in the actual Eddington luminosity – conveniently modified by a factor and, therefore, convergence for higher values of luminosity in the intershell region.

Fig. 6 Left upper panel: mass coordinate versus density near the κFe peak. Left lower panel: mass coordinate versus temperature near the κFe peak. Right panel: opacity profile versus mass. The insert corresponds to the κFe peak.

Fig. 7 First panel: Log Pgas; second panel: Log T; third panel: κ; fourth panel: ϵbind versus the logarithm of radius, Log R near the zone of the instability for some selected models prior to the divergence. The arrows indicate the direction in which time increases.

3. Test of the Fe-peak opacity hypothesis

One possible cause for the instability that stops the evolution of our model massive AGB and super-AGB stars is the effect of the iron opacity peak (κ_Fe peak). The κ_Fe peak is not a unique feature in OPAL tables (Iglesias & Rogers 1996b). This feature can also be found in opacities from the Opacity Project (OP). Seaton & Badnell (2004) has a detailed comparison between OPAL and OP opacities. Typically, the Fe-bump occurs at a temperature around 2.0 × 10⁵ K, this feature being more pronounced as the metal content increases. The Fe-bump occurs at a higher temperature in OP than in OPAL (Jeffery & Saio 2006).

The effect of this κ_Fe peak in the WR stars, studied by Petrovic et al. (2006), was to decrease the efficiency of energy transport near the local opacity maximum. For stars more massive than about 15 M_⊙ and close to the Eddington limit such a decrease in the efficiency of energy transport would lead to local values of the luminosity above L_Edd in the zones just below the opacity peak. These authors consider that such situation would lead to a significant outwards acceleration and the inflation of the stellar envelope. Even though the stars we are considering here are in the intermediate-mass range and are in a different evolutionary stage, we also find the κ_Fe peak (see Figs. 6 and 7) at a temperature of t_peak = 1.6 × 10⁵ K, so also close to the value found by Petrovic et al. (2006). Furthermore our last converged models also show a fast increase in the surface radius from 1000 R_⊙ to about 1700 R_⊙ in a time interval of 10 years approximately. Therefore, our model star also experiences inflation. Gräfener et al. (2012) have recently studied inflation in Wolf-Rayet stars and luminous blue variables and established a relation with the topology of the κ_Fe peak.

In our case the κ_Fe peak develops after the helium flash and advances inwards, both in mass and radius coordinates, as we can see in panel 3 of Fig. 7. In this figure we can see the profiles of temperature, gas pressure, opacity and binding energy versus the logarithm of the radius, for some selected models between the last helium flash and the time of code divergence. As we have explained, the opacity peak is associated to a certain temperature T_peak. This means that an advance inwards of the κ_Fe peak is directly connected to an advance inwards of T_peak and, therefore, to the cooling and expansion of the affected zone of the star. The existence of κ_Fe peak reduces drastically the efficiency of energy transport and thus favours the trapping of the energy generated at the helium-burning shell. This implies a local increase of internal energy u_int and, as can be seen in panel 4 of Fig. 7 an increase in binding energy ϵ_bind, where ϵ_bind = ϵ_grav + u_int, and ϵ_grav is the gravitational energy. u_int increases at a zone that is expanding while cooling is caused by recombination, and this process is fed-back, that is, more cooling favours more recombination that implies more energy for expansion and cooling and so on.

Eventually the densities between the helium-burning shell and the BCE decrease so much and so fast that P_gas shows an inversion, that is, an increase with increasing radius that we can see in the first panel of Fig. 7. This can be identified as the start of the instability that finally leads to the crash of the code. The fact that the instability is associated to a decrease in the efficiency in energy transport is supported by the fact that the first model presenting P_gas inversion corresponds to the time at which the inner convective shell disappears.

This explanation is consistent with Sect. 2. On one hand, the zone of instability identified as the zone of reduced density and P_gas – and even P_gas inversion, coincides in mass with the zone of instability where β tends to zero. Furthermore, the peak in binding energy also coincides with the zone where L_r > L_Edd,conv. On the other hand, the instability described in this section can be delayed, not only by artificially removing the opacity peak κ_Fe, but also by increasing the α parameter. Both increasing α or decreasing κ would allow increasing the efficiency of energy transport and the effect of blocking of energy would be milder.

The local maximum in the opacity κ_Fe might be responsible for a behaviour similar to the one caused by the κ-mechanism in Cepheids and RR-Lyrae. For these pulsating stars, the hydrogen or helium ionisation regions are responsible for a local opacity increase and the loss of efficiency in energy transport. This causes an energy excess below the opacity maximum that, eventually, can be transformed into work of expansion in the envelope. Once the envelope expands, density and temperature decrease and the degree of ionisation also decreases, the star contracts again and the process repeats.

Fig. 8 First panel: logarithm of temperature, second panel: specific binding energy ϵbind, third panel: specific internal energy uint and fourth panel: logarithm of opacity profiles for a few selected models along the second last thermal pulse, represented in Fig. 2.

As our code is not hydrodynamic, we cannot follow the evolution further. Nevertheless it seems reasonable to assume that, as in the massive stars studied by Petrovic et al. (2006), the energy accumulated below the κ_Fe peak will eventually be released and transformed into work of expansion, so the (super)-AGB envelopes will also inflate. During this phase of inflation the mass-loss rates will be correspondingly high and therefore an enhanced superwind is likely to develop.

We have tested the hypothesis of the κ_Fe peak being responsible for the divergence of our code by artificially removing this peak. We have used Figs. 6 and 7 and inspected the data files of the structure profiles of the models prior to divergence. Using this information we have imposed that the opacity for temperatures between 1.5 × 10⁵ K and 3 × 10⁵ K and densities between 10^-11 and 5 × 10^-7 g/cm³ must have a constant value equal to the opacity at 3 × 10⁵ K and density 5 × 10^-7 g/cm³. By imposing this constant value, we avoid the instability. This suggests the Fe peak is a major reason for this instability, but we must take into account that the κ-peak is real and, unless the energy transport is much more efficient than what we expect from the mixing length theory, the peak does exist and is likely to lead, as we have seen, to the inflation and, perhaps, even to the disruption of the star.

4. What will happen to the star?

Wood & Faulkner (1986) emphasise that the instability that affects TP-AGB stars with degenerate cores more massive than 0.89 M_⊙ results from the disappearance of a hydrostatic solution to the stellar structure equations. Hence they postulate that the resulting hydrodynamic event will remove the stellar envelope.

We have estimated the outward velocity of the envelope by comparing the radii at the same mass mesh points between the last few converging models. We found that almost throughout the whole envelope, the shell velocities exceed the local escape velocities. But it is important to note that, because the code crashes during the ejection (as seen from the expansion of the radius of the star), the actual velocity during ejection could be even be higher than the one we are estimating. This problem, in fact, should be tackled by using hydrodynamical calculations.

Fig. 9 First panel: logarithm of temperature, second panel: specific binding energy ϵbind, third panel: specific internal energy uint and fourth panel: logarithm of opacity profiles for a few selected models along the last thermal pulse, represented in Fig. 3.

We have also approached the problem of the outcome of the instability by considering the binding energy of the envelope. Figure 8 shows the profiles of Log T, specific binding energy ϵ_bind = u_int + ϵ_grav, specific internal energy u_int and opacity Log κ for the models of the second last thermal pulse labelled in Fig. 2, and using the same color code. Figure 9 shows the same for some selected models of the last thermal pulse, as labelled in Fig. 3. The specific internal energy u_int has been computed taking into account the terms associated to the ideal gas, to radiation, to ionization and dissociation and to electron degeneracy, as in Han et al. (1994). Wherever ϵ_bind has a high negative value the stellar shells are strongly bound and wherever the binding energy approaches zero, the shells are loosely bound. Positive values of the binding energy correspond to the case in which the shells would have enough energy to escape, which can be seen in Figs. 7 and 9. The positive values in ϵ_bind are located between the top of the helium-burning shell (1.1700 M_⊙ approximately), and with the base of the convective envelope (1.1706 M_⊙ approximately) in Fig. 9.

The most representative features of the internal energy profiles (third panels of Figs. 8 and 9) are the the local maxima that correspond to the zones where energy is accumulated due to the κ_Fe peak, as we have explained in the former section. As we can see, these peaks practically coincide in mass point with the drop in temperature at the base of the convective envelope and the jump in ϵ_bind.

The second last thermal pulse is characterised by the fact that the mass point of the base of the convective envelope and, therefore, also the u_int peaks and ϵ_bind and κ peaks are located practically at the same mass point before the helium flash and until the inner convective shell disappears. Along the last thermal pulse, though, these masses reach values closer to the centre. Finally, when the base of convection reaches the mass point M = 1.17056 M_⊙ and temperatures between 0.1 MK and 1.5 MK, the opacity peak narrows but reaches a diverging value. Due to this u_int diverges, and so does ϵ_bind. This is the point at which the it is not possible to proceed further with our calculations.

Another important feature that distinguishes the second last thermal pulse from the last one is the variation of the surface radius (see Figs. 2 and 3). The second last thermal pulse behaves in a standard way. When the helium flash develops the HeBS expand and pushes the HBS that also expands and cools down. This leads to the switching off of the HBS and, without its energy supply gravity dominates and the stellar envelope contracts. The stellar radius only recovers and increases its value when the helium flash is over and the luminosity provided by hydrogen burning is again higher than the luminosity due to helium burning. On the other hand, during the last thermal pulse the stellar radius remains practically constant. Therefore the energy provided by helium burning is not absorbed in the upper HeBS and intershell region. Instead, it is sufficient to sustain the star against gravity. This situation is stable as long as the inner convective shell is present. But once this disappears while helium burning is still strongly active, a runaway process occurs that leads to the increase of the stellar radius on dynamical time scales and, eventually, to the instability described in this work and to the divergence of the calculations.

From this simple approach, we can derive the same conclusion as Wood & Faulkner (1986), that is, the envelope will be ejected shortly after the instability is reached. In any case, we would like to point out the limitations of our hydrostatic approach to a problem that would be better treated with hydrodynamical codes. Besides, it is also important to realise that the envelope masses we are considering in this work −≲2 M_⊙, are much higher than the ones considered in Wood & Faulkner (1986). Therefore, the amount of work necessary to eject completely the stellar envelopes of massive AGB stars is much higher and, even though our panels showing the binding energy profiles seem to point to the disruption of the star, we cannot discard an eventual fallback of the envelopes. We have calculated the envelope binding energy for the last computed models as Meng et al. (2008), and obtained a result of −3 × 10⁴⁵ erg. This value is negative and, therefore technically the envelope still remains bound. But this value should be considered in perspective. First, the total binding energy of the star is of the order of 10⁴⁹. Therefore the envelope mass (about half the mass of the star at this point) has only between 10^-4 and 10^-3 of the total binding energy, so in practice the envelope is very loosely bound to the core. Second, we can see that the binding energy becomes less negative during the last computed models. Thus the envelope becomes less bound with time, until the evolution stops and we cannot proceed further, but this trend seems to be a good hint of the future behaviour of the star. A hydrodynamic code would be more suitable to accurately model this ejection process.

If our estimates showing that most of the envelope (about ≲2 M_⊙) will be ejected are correct, the outcomes of our model stars might be relatively massive planetary nebulae, with some hydrogen-rich matter surrounding a central star of around 1 M_⊙. However, it is not clear for how long this post-ejection system could be observed as planetary nebula because it could move very quickly across the HR-diagram to the white dwarf cooling track (Vassiliadis & Wood 1994). The evolution of a 1 M_⊙ post-AGB remnant may be too fast to light a massive PN, and even if all massive AGB stars undergo the ejection mechanism, they are unlikely to be observed as a PN with a massive nebula. Yet, it is interesting to note that the observational counterpart of such a system could indeed exist. Planetary nebula N66 contains about 0.6 M_⊙ of hydrogen-rich matter in the nebula, and the luminosity of the post-AGB star corresponds to a core mass of about 1.2 M_⊙ (Hamann et al. 2003). Indeed, Vassiliadis (1996) has already suggested that this planetary nebula could be caused by the radiation ejection mechanism. The very low frequency of N66-like PN with a large nebula could be consistent with the fact that not all AGB stars going through such an ejection mechanism would be observed as a PN.

5. Discussion and conclusions

In this work we consider the problem of the instability of massive AGB and super-AGB stars and identify the reasons why stellar evolution codes break when envelope masses are still relatively massive – as massive as 2 M_⊙ for the case of super-AGB stars. If this instability is actually encountered by real stars and if it immediately causes their total disruption then it would affect the expected yields of intermediate-mass stars. Furthermore, the results of our analysis may have an effect on the determination of the mass limit for the formation of electron-capture supernovae, as the presence of the instability we have described might quench the evolution of the star toward core masses above Chandrasekhar mass. Such an analysis might be particularly important in the case of extremely metal-poor super-AGB stars, in which computational time is still a problem and whose final fate determination has been often affected by very simple extrapolations.

The problem of the instability of stars at the end of the AGB phase has been considered by Wood & Faulkner (1986). These authors realised that, near the maximum of a late helium flash during AGB evolution, local values of the radiation pressure became so high that they were enough to balance gravity. Under these circumstances the gas pressure tended to zero and the local luminosity (L_r) became higher than the Eddington luminosity (L_Edd). Therefore the model stars departed from hydrostatic equilibrium and their evolution was quenched. We have computed the evolution of massive AGB stars and realised that the zones at which L_r > L_Edd occurred were at the base of the convective envelope. Therefore we have generalised the results by Wood & Faulkner (1986) for the case in which the departure of hydrostatic equilibrium developed in the presence of convection and have generalised an expression of L_Edd for convection-dominated environments.

Opacities have frequently been suggested as responsible for the instability. For instance, in the classic works by Baker & Kippenhahn (1965) and Baker (1966), it is justified that the conditions for the stability of static, radiative layers in gas spheres only depend on the local thermodynamical state of these layers – Baker’s one-zone model. Thus, if the constitutive relations – equations of state and Rosseland mean opacities – are specified, the stability conditions can be evaluated without specifying further properties of the layer. In the case of solar-composition gas spheres the κ-mechanism can even work in regions where H₂ dissociation and H-recombination are occuring, and such regions are much more extended than those of the second ionization zone that drives Cepheïd pulsations. Wagenhuber & Weiss (1994) also pointed to H-recombination as responsible for the termination of the AGB phase.

We show that, in the case of AGB and super-AGB models at the late stages of their thermally pulsing phase a similar instability occurs, but it is caused by the Fe-peak in the opacity at the base of the convective envelope. This phenomenon was already described by Petrovic et al. (2006) for the case of massive stars. As in Wood & Faulkner (1986), the radiation pressure completely dominates the gas pressure and density at the convective envelope converges to zero. When this instability occurs depends on the treatment of convection. A more efficient scheme of convection or overshooting could delay or even avoid its occurrence. In particular, the value of the mixing length parameter α determines when the instability occurs. This point is very relevant, as it has recently been suggested that α should vary during the evolution of the star – Meakin & Arnett (2007). Moreover, it is possible that mixing length theory cannot adequately describe the convection of the actual star – Stancliffe et al. (2011). Alternatively, if we could find some observational constraints on these models, we might be able to use these to discriminate between different convection theories or different values of α for this phase of the evolution.

We have seen that for massive AGB stars or super-AGB stars with zero-age main-sequence masses 7−10 M_⊙ the instabilities occur when the envelope mass is 1 − 2 M_⊙. The instability could occur for initial masses as low as 2.5 M_⊙ in the solar metallicity case. The stars encountering the instability in this study are all hot bottom burning AGB stars. Therefore, the envelope ejected should be nitrogen rich. The envelope mass when instability occurs decreases with the masses of the AGB. For low mass AGB stars, whether the instabilities occur depend on the choice of the mixing length parameter α and treatment of convection. It is fair to conclude the instability does occur for massive AGB stars and super-AGB stars unless mixing length theory is fundamentally flawed.

Let us speculate briefly on the behaviour of the stellar material in the region of the instability. Our estimates of the velocity of the convective envelope are greater than the escape velocity. Furthermore, when we consider the binding energy profiles for our last converged models we can see that binding energy reaches positive values. This suggests that almost all of the envelope could be ejected. In fact the radius of the whole star increases very rapidly. It is thus possible that the AGB stars will enter post-AGB phase right after the instability. As to the effects of our analysis on the computed yields of AGB and super-AGB stars, we can expect that the ejection phase, if it occurs, will be very short compared to the whole TP-(super)AGB phase and that, at this point, the nuclear reactions do not play an important role. The nucleosynthesis of the AGB stars will be truncated due to the ejection and the yield of the s-process isotopes will be lower compared to the alternative assumptions that thermal pulses will keep occuring.

If ejection does occur after the instability, the TP-(super)AGB evolution would be truncated much earlier than in a model star that underwent “normal” evolution. This would prevent the occurrence of electron-captured supernovae. Compared to the result of Gil-Pons et al. (2007), Siess (2007) and Poelarends et al. (2008), which deduce the possibility of electron-captured supernovae based on the assumption of regular thermal pulses and extrapolation on core growth rate and mass-loss rate, the actual ranges for the occurrence of electron-captured supernovae could be much narrower in terms of both mass and metallicity. For example, our models show that electron-capture supernovae is very unlikely in solar metallicity. On the other hand, if the instability merely leads to a higher mass-loss rate instead of the ejection of the whole envelope, then electron-capture supernovae may still be possible even though the instability is recurrent. The effect of the instability is certainly a factor that cannot be discarded when considering the possibility of electron-captured supernovae. On the other hand, at lower metallicities, the opacity jump due to the peak would be reduced because of the lower abundance of iron, thus the instability might be reduced and the ejection occurs much later or not occur at all for lower metallicity models. In such a case, if the instability described by Wagenhuber & Weiss (1994) cannot operate either, it might possible to find a critical maximum metallicity for electron-capture supernovae to occur. Evolution of metal-poor super-AGB stars will be modeled in future work.

As a final note, our hydrostatic code could not model through the whole episode of the instability. Although our estimate shows the envelope is likely to be ejected, there are a few processes that cannot be modelled properly, for instance, the situation when the envelope expands and the density converges to zero at the bottom of the convective envelope. As a result, the density of the envelope is inversed from the bottom and induce the Rayleigh-Taylor instability, that could facilitate mixing of material and energy transport. Thus the instability may be overestimated by our code. However, modelling such a process is beyond the capability of the Monash evolutionary code. Modelling on a 3D hydrodynamic code such as Djeuhty

(Dearborn et al. 2006) may be needed in order to study this instability properly.

Acknowledgments

This work was supported by the Australian Research Council through grants DP0877317 and DP1095368. H.B.L. would like to thank Norbert Langer for interesting discussion on instability of massive stars. H.B.L. would like to thank Christopher Tout for comments on energy transport. P.G.P. would like to thank the kindness and hospitality of the SINS group at Monash and Jordi Ortiz Domenech for his constant support. We would like to thank the anonymous referee for suggestions improving our manuscript.

References

All Figures

	Fig. 1 Density, temperature, radius, luminosity, κ and β, the ratio of gas pressure to total pressure in the instability region. The model profiles represented in solid lines correspond to the times labelled in Fig. 3 and in particular thick lines represent the last converged models. Dashed lines represent the final progress of the instability.
In the text

	Fig. 2 Upper panel: evolution of luminosity associated to H-burning (dotted line) and He-burning (solid line); middle panel: evolution of convection; lower panel: evolution of the surface radius, all during the second last thermal pulse. Relevant structure profiles of models labelled a to f appear in Fig. 8.
In the text

	Fig. 3 Upper panel: evolution of luminosity associated to H-burning (dotted line) and He-burning (solid line); middle panel: evolution of convection; lower panel: evolution of the surface radius, all during the last thermal pulse. Relevant structure profiles of models labelled a to f appear in Fig. 9.
In the text

Fig. 4 Relevant structure and composition parameters for the model labeled c in Fig. 3. Upper panel: actual gradient (thin solid), radiative gradient (dotted line) and β (thick solid). Middle panel: local luminosity values (dotted line), Eddington luminosity (thin solid) and Eddington luminosity corrected with the effect of convection (thick solid line). Lower panel: hydrogen (solid line), helium (dotted line) and carbon (dashed line) composition profiles.

In the text

Fig. 5 Relevant structure and composition parameters for the model labeled f in Fig. 3. Upper panel: actual gradient (thin solid), radiative gradient (dotted line) and β (thick solid). Middle panel: local luminosity values (dotted line), Eddington luminosity (thin solid) and Eddington luminosity corrected with the effect of convection (thick solid line). Lower panel: hydrogen (solid line), helium (dotted line) and carbon (dashed line) composition profiles.

In the text

	Fig. 6 Left upper panel: mass coordinate versus density near the κFe peak. Left lower panel: mass coordinate versus temperature near the κFe peak. Right panel: opacity profile versus mass. The insert corresponds to the κFe peak.
In the text

	Fig. 7 First panel: Log Pgas; second panel: Log T; third panel: κ; fourth panel: ϵbind versus the logarithm of radius, Log R near the zone of the instability for some selected models prior to the divergence. The arrows indicate the direction in which time increases.
In the text

	Fig. 8 First panel: logarithm of temperature, second panel: specific binding energy ϵbind, third panel: specific internal energy uint and fourth panel: logarithm of opacity profiles for a few selected models along the second last thermal pulse, represented in Fig. 2.
In the text

	Fig. 9 First panel: logarithm of temperature, second panel: specific binding energy ϵbind, third panel: specific internal energy uint and fourth panel: logarithm of opacity profiles for a few selected models along the last thermal pulse, represented in Fig. 3.
In the text