© The Authors 2026

1. Introduction

The turbulent transport of heat, angular momentum, and chemical species stands among the most important unsolved problems in the modeling of stellar interiors (e.g., Arnett et al. 2015; Salaris & Cassisi 2017; Xiong 2021; Anders et al. 2022; Sibony et al. 2023). Convective energy transport, as well as mixing processes in stellar astrophysics, are usually described with the help of the so-called mixing length theory (MLT; Prandtl 1925; Biermann 1932; Vitense 1953; Böhm-Vitense 1958). This theory, together with the so-called Schwarzschild instability criterion (Schwarzschild 1906), provides a simple algorithm to treat convection locally at any given point in the star. Despite its simplicity, it shows a good agreement with observations. In its standard and most adopted form, MLT does not take into account the impact of chemical gradients on the development of convective motions (Garaud 2020), nor does it describe the more subtle double-diffusive mixing processes, such as thermohaline mixing (Kippenhahn et al. 1980) and semiconvection (Kato 1966; Langer et al. 1983).

Grossman et al. (1993) and Grossman & Taam (1996), henceforth referred to as GNA, developed a theory that considers the Rayleigh-Taylor (RT) instabilities as well as the double-diffusive mixing processes that can develop in the presence of chemical gradients. This is a nonlocal description of convection equipped with its own stability conditions and successfully reproduces most known results in stellar convection (Wachlin et al. 2011). However, GNA does not consider the possible electronic degeneration and nonideal effects in the stellar medium and in evolved stars, where the deep interior becomes degenerate, this description can be insufficient as well. Moreover, the mathematics involved in GNA are far more complicated and numerically unstable than standard MLT; thus, it becomes a problematic issue when introducing this theory into stellar evolution codes.

Several works have extended MLT to deal with convection in the presence of background chemical gradients in the past. For example, Umezu & Nakakita (1988) extended the MLT to study the impact of positive (stabilizing) background chemical gradients in the development of convective cores in massive stars. More recently, Castro-Tapia et al. (2024) adopted a similar model to explore convective and nonconvective motions above the crystallization front in white dwarfs. In particular, their work focuses on the final steady state that develops for given fluxes of both heat and chemical elements and they derived chemical gradients for their solutions for given chemical fluxes. In this work, we adopt a similar treatment to that of both works, essentially using the same mathematical prescription, but with a different scope and setting of the relevant physical conditions. As in Umezu & Nakakita (1988), we are interested in the development of convection in the presence of background chemical gradients, yet we are mostly interested in the destabilizing effects of negative chemical gradients that may develop in the late stages of low-mass stars and their consequences for the final chemical stratification and pulsations of pre-white dwarfs and white dwarfs. Consequently, we are mostly interested in the derivation of well-behaved equations that can be included, at minimal numerical cost, into stellar evolution computations.

Asymptotic giant branch (AGB) stars are evolved objects with initial masses under 8–10 M_⊙ (Kippenhahn et al. 2013; Herwig 2005) that have already ended their helium core burning phase. They are evolved red giant stars, with an inert carbon-oxygen core surrounded by an inner helium burning shell and a more external hydrogen burning shell. The later phase of this two-shell burning configuration is not quiescent, but rather dominated by the development of events commonly termed thermal pulses. Thermal pulses are quasiperiodical thermal runaway events that develop in the helium burning shell and lead to the formation of an inner convective zone in the region located in between the two burning shells (usually named the intershell). This short burst expands the outer layers, cooling the hydrogen burning shell, resulting in a diminution of the surface luminosity. This is followed by a longer period of the system slowly contracting again and going back to the initial configuration. Each pulse leaves behind a strong stratification in the carbon-oxygen profile, which could produce inversions of the chemical gradient leading to the development of thermohaline mixing or, alternatively, Rayleigh-Taylor instabilities and the development of chemically driven convection. This situation makes AGB stars an ideal astrophysical location for testing extensions of the MLT where chemical gradients play an important role.

Mixing processes on the AGB shape the internal structure of their descendants: white dwarfs (WDs) and pre-WDs (Althaus et al. 2010; Córsico et al. 2019). Among them, one particularly interesting population for the present study is the group of GW Vir variable stars, named after the prototype of the class PG 1159-035 (McGraw et al. 1979; Oliveira da Rosa et al. 2022). These stars are hot, luminous stars, characterized by He-, C-, and O-rich atmospheres and by the occurrence of multiperiodic luminosity variations attributed to nonradial gravity (g) pulsation modes. GW Vir stars are of great interest for this work for two main reasons. First, GW Vir stars are known to display large numbers of periods corresponding to large numbers of normal oscillation modes. The presence of a large number of normal modes allows for the use of asteroseismological techniques to probe the internal chemical structure of these stars. Secondly, these stars departed the AGB only tens of thousands or even only thousands of years ago (Miller Bertolami 2024). Consequently, many of the physical processes that alter the chemical stratification in older WDs, such as chemical diffusion and gravitational settling (Althaus et al. 2010), did not have enough time to act in these stars. This combination of a recent departure from the AGB and the large number of pulsation modes observed in their lightcurves, makes GW Vir stars excellent candidates to test the impact of convective and nonconvective mixing processes on the AGB. Asteroseismology provides an ideal tool to unveil the internal chemical structure of stars, by comparing the observed periods with the periods predicted by theoretical models (Córsico et al. 2019). Moreover, the total number of observations, quantity of detected modes, and quality of observations in GW Vir pulsators has strongly increased with the advent of space-based observatories (Córsico et al. 2021; Oliveira da Rosa et al. 2022; Calcaferro et al. 2024a,b)

This work is organized as follows. In Section 2, we describe the MLT extension we developed to include chemically driven convection and thermohaline mixing. In Section 3, we describe how we included this new set of equations into a stellar evolution code to compute the evolution of AGB and post-AGB stars. We apply this prescription, compare it with the standard MLT approach and with GNA theory, and analyze the differences that arise in the chemical structure of the same star under the different prescriptions for convective (and nonconvective) motions. In Section 4, we explore the asteroseismological and observable consequences of these different chemical evolutions. Finally, in Section 5, we summarize our conclusions and discuss the future outook.

2. A simple extension to mixing length theory

Currently, MLT is the most common and used model of convection in stellar astrophysics given to its simplicity and its good agreement with respect to observational data. Conceptually, it proposes the existence of macroscopic mass elements (or “blobs” of matter) that move and transfer heat to their surroundings. Mass elements are thought to travel a given distance (the so-called mixing length) before dissolving to the surroundings.

Convection theories are commonly used in tandem with instability criteria, which determine in the regions where turbulence develops and transports heat and/or chemical species. Radial displacements are stable against the temperature stratification (Kippenhahn et al. 2013) if¹

$$\begin{array}{c}\hfill \mathrm{\nabla }<{\mathrm{\nabla }}_{\text{e}}+\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu },\end{array}$$(1)

where

$$\begin{array}{c}\hfill \mathrm{\nabla }\equiv (\frac{dlnT}{dlnP}{)}_{\text{s}},{\mathrm{\nabla }}_{\text{e}}\equiv (\frac{dlnT}{dlnP}{)}_{\text{e}},{\mathrm{\nabla }}_{\mu }\equiv (\frac{dln\mu }{dlnP}{)}_{\text{s}},\end{array}$$(2)

with the subindex “e” referring to the displaced element of stellar matter and “s” referring to the surroundings, respectively, and

$$\begin{array}{c}\hfill \phi \equiv (\frac{\partial ln\rho }{\partial ln\mu }){|}_{T,P},\delta \equiv -(\frac{\partial ln\rho }{\partial lnT}){|}_{P,\mu },\end{array}$$(3)

where T is the temperature, P the pressure, and μ the mean molecular weight. It is noteworthy that even when Eq. (1) is fulfilled, turbulence might arise. A clear example of this was demonstrated by Kato (1966) who shown that, under the right conditions, heat diffusion can lead to the development of an oscillatory behavior with an increasing amplitude (i.e., semiconvection).

Usually, the values of the actual temperature gradient of the stellar material (∇) and convective elements (∇_e) are not known beforehand and Eq. (1) cannot be used as a stability criterion². Consequently, it is a standard practice to determine convective regions by checking which regions of the star cannot be stable and transport heat solely by radiation and conduction. Moreover, since ∇_e is also not known beforehand, the most common practice is to check in which regions mass elements displaced adiabatically become dynamically unstable (∇_e = ∇_ad). Under these assumptions, we can obtain the so-called Ledoux stability criterion (Ledoux 1947; Sakashita & Hayashi 1959) via

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\text{rad}}<{\mathrm{\nabla }}_{\text{ad}}+\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }.\end{array}$$(4)

The radiative gradient ∇_rad is defined as the temperature gradient that would be required if the heat were to be transported only by radiation (or conduction). Moreover, if we neglect the impact of the chemical gradients (as the standard application of MLT does after), we then have the Schwarzschild stability criterion (Schwarzschild 1906), namely,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\text{rad}}<{\mathrm{\nabla }}_{\text{ad}}.\end{array}$$(5)

In MLT, the velocity of the convective blob (Kippenhahn et al. 2013, see also Appendix A) is calculated as

$$\begin{array}{c}\hfill v^2=g\delta (\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}})\frac{l_m^2}{8H_P},\end{array}$$(6)

with g being the gravitational acceleration, l_m represents the ML of MLT and H_P is the scale height of pressure, defined as

$$\begin{array}{c}\hfill H_P\equiv -\frac{\mathit{dr}}{dlnP}·\end{array}$$(7)

In MLT, we aim to obtain the values for ∇ and ∇_e and, with that objective, the following two equations can be derived:

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_{\text{e}}-{\mathrm{\nabla }}_{\text{ad}}& =2U\sqrt{\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill {(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}})}^{3/2}& =\frac{8U}{9}({\mathrm{\nabla }}_{\text{rad}}-\mathrm{\nabla }),\hfill \end{array}$$(9)

with

$$\begin{array}{c}\hfill U\equiv \frac{3ac{T}^3}{c_P{\rho }^2\kappa l_m^2}\sqrt{\frac{8H_P}{g\delta }}·\end{array}$$(10)

In this expression, a is the radiation constant, c the light velocity, c_P the specific heat at constant pressure, ρ the density, and κ the opacity. In MLT equations, when the Schwarzschild criterion is considered and the usual quantities as ∇_rad and ∇_ad are given, we can obtain unique real solutions for ∇ and ∇_e and the convective velocity, v, can be obtained on their basis. However, if we consider the Ledoux criterion instead of a Schwarzschild, a region that is Ledoux unstable but Scharzschild stable³ appears. The new unstable region corresponds to a Rayleigh-Taylor (RT) instability.

We can obtain expressions similar to Eqs. (8) and (9) without neglecting the impact of the chemical gradients to describe the convective motions. The modified equations for ∇, ∇_e, and v are (see Appendix A) as follows,

$$\begin{array}{cc}& ({\mathrm{\nabla }}_{\text{e}}-{\mathrm{\nabla }}_{\text{ad}})(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }{)}^{1/2}=2U(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}),\hfill \end{array}$$(11)

$$\begin{array}{cc}& (\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }{)}^{1/2}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}})=\frac{8}{9}U({\mathrm{\nabla }}_{\text{rad}}-\mathrm{\nabla }),\hfill \end{array}$$(12)

$$\begin{array}{cc}& v^2=g\delta (\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu })\frac{l_m^2}{8H_P}·\hfill \end{array}$$(13)

This set of equations is equivalent to Eq. (25) of Umezu & Nakakita (1988) and Eq. (20) of Castro-Tapia et al. (2024). Neither of these works used the derived mixing velocities and temperature gradients as part of stellar evolution computations but to study specific processes. On the one hand, Umezu & Nakakita (1988) studied to which extent convection was possible in the presence of positive (stabilizing) chemical gradients. They concluded that in the presence of finite thermal diffusion (regardless how trivial), the instability criterion derived under the assumption of adiabatic perturbations overestimates the instability region. On the other hand, Castro-Tapia et al. (2024) studied the expected steady state that develops in the interior of a white dwarf, under the assumption of given values of the heat and chemical fluxes. Our objective, however, is to use Eqs. (11) and (12) as part of a stellar structure and evolution code to compute the mixing and heat transport in the presence of chemical gradients caused by nuclear burning. We refer to this simple extension as MLT♯ and its derivation is described in Appendix A.

Having obtained equations for convective motions considering Ledoux criterion and the effect of chemically driven convection (RT instabilities), we still need to consider one more scenario. First, we know that a negative ∇_μ will favor mixing due to the tendency of the heavier layer to sink under the lighter one. However, regions with ∇_μ < 0 and Ledoux stable are possible when the heavier element is sufficiently hotter than the light layer below. While it is dynamically stable, over time, thermal diffusion will act, cooling the blob, and the element will then sink and mix with the surrounding environment. This is the regime of thermohaline mixing and this process is significantly slower than dynamic convection. Nonetheless, the timescales are still much shorter than typical stellar evolution timescales and, thus, thermohaline mixing becomes relevant, for example, in stars on the red giant branch (RGB; Wachlin et al. 2011) and during the thermal pulses on the AGB.

It is noteworthy that in this model, Eqs. (11) and (12) give solutions beyond the limit established by the Ledoux criterion, entering the regime of thermohaline mixing, as long as $\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }>0$ and ∇_rad < ∇_ad. However, Castro-Tapia et al. (2024) showed that when assuming inefficient convection, the mixing velocity of the extended MLT model in the sub-Ledoux regime takes the same form as in the parametrization given by Kippenhahn et al. (1980) for thermohaline mixing. Below the critical line defined by Eq. (4), we find that Eqs. (11) and (12) might have multiple solutions, ranging from an almost fast-adiabatic convection to a slow thermohaline behavior with ∇ ≃ ∇_rad. In these circumstances we decided for the most conservative approach and decided to adopt the thermohaline solution for all cases below the critical line given by Eq. (4) to deviate the least from previous works. We will employ then, for this region of the parameter space, the following expression given by Kippenhahn et al. (1980) for the mixing velocity,

$$\begin{array}{c}\hfill v_T=\frac{6ac{T}^3}{c_P{\rho }^2\kappa }\frac{\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }}{(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{ad}})}·\end{array}$$(14)

For the temperature gradient, in this case, we consider ∇ = ∇_rad.

As stated before, both RT and thermohaline instabilities are considered stables under Schwarzschild criterion (SC) and, then, standard MLT will not mix regions that are neither RT, nor thermohaline.

The extension of the standard MLT that we present in this paper, which we refer to as MLT♯ is defined by the following prescriptions: (i) standard MLT when the regions are Schwarzschild unstable; (ii) Eqs. (11), (12), and (13) when the region is Schwarzschild-stable, but Ledoux-unstable (RT instabilities); and (iii) thermohaline mixing given by Eq. (14) when the layer is Ledoux-stable, but ∇_μ < 0. In the next section, we describe how we applied the MLT♯ model in AGB stars and we compare the results with using standard MLT and GNA prescriptions.

3. Mixing processes in the asymptotic giant branch

After the helium core-burning phase, stars with initial masses of M_in < 8–10 M_⊙ move to the AGB in the Hertzprung-Russell (HR) diagram (Kippenhahn et al. 2013). These stars posses an inert carbon(C)-oxygen(O) core surrounded by two concentric shells where helium and hydrogen burn. The advancing He-burning shell becomes progressively thinner, leading to an unstable condition. This unstable setup will lead to the development of thermal instabilities in the He-burning shell. These thermal instabilities, known as thermal pulses (TP), lead to the formation of inner convective regions where the material is periodically mixed. Each thermal pulse will leave behind a chemical stratification with a sharp O abundance peak (Fig. 1, blue lines) at the outer regions of the CO core (Herwig 2005; Kippenhahn et al. 2013).

Fig. 1. Chemical abundances of He (dotted lines), C (continuous lines) and O (dashed lines) of stellar models of initial mass Min = 1 M⊙ (left panel), 1.5 M⊙ (middle), and 3 M⊙ (right) at the end of the AGB phase. The blue lines correspond to the case where MLT was used, the orange ones are for MLT♯, and the black ones are for GNA.

These peaks in the O abundance could lead to inversions in the chemical gradient and then those layers in the star will be either RT unstable or unstable against thermohaline processes. As stated before, these regions are not mixed in standard stellar evolution computations (MLT + SC). This scenario is then ideal for testing the development of RT or thermohaline instabilities and to study the impact of chemically driven convection (CDC) and thermohaline mixing in evolved stars. Moreover, the existence of variable stars immediately after this stage (the GW Vir stars mentioned in the introduction) opens up the possibility to test these predictions and, in the future, calibrate the intensity of those mixing processes.

To explore the impact of mixing processes induced by negative chemical gradients, we carried out simulations for three different initial masses: M_in = 1 M_⊙,  1.5 M_⊙,  and 3 M_⊙ using LPCODE (Miller Bertolami 2016; Althaus et al. 2020). Since we are interested only in the physical processes that take place during the AGB phase⁴, all the prior evolution was carried using standard MLT and Schwarzschild criterion to avoid numerical discrepancies. When the stars were located into the AGB, we split the sequences into three different paths: one with standard MLT, one with MLT♯, and a third with GNA. The GNA test was done with the purpose of comparing the outcome of the simple MLT♯ with that of a significantly more complex (and, thus, more complicated and numerically problematic) theory. Moreover, in this paper, we limit our analysis to the outer regions of the core and the intershell. Possible mixing processes at the innermost regions of the CO core immediately after the end of core He burning will be addressed in future works. We note here that both GNA and the current version of MLT♯ (which considers $\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }$ and not the complete Ledoux term B) are inadequate when strong degeneracy of the electron gas and nonideal interactions become important.

In all the sequences, we considered an initial metallicity of Z = 0.015 and mixing length of l_m = 1.822. The rest of our choices for the input physics are the same as in Miller Bertolami (2016).

As stated before, during each thermal pulse, a discontinuity in the carbon-oxygen chemical profile is created, which is peak-shaped as we can see in Fig. 1 (see the blue lines). However, this does not directly imply an inversion of the chemical gradient in the region, as long as a remnant abundance of helium capable of stabilizing those layers is present. Only after the helium in this region is burned ∇_μ becomes negative and MLT♯ and GNA mix the elements and smooth the profiles. MLT never mix those regions, as explained before. In Fig. 2 we show a comparison of the chemical evolution using MLT and MLT♯ during the same TP for the case of M_in = 1.5 M_⊙. We find that the timescale for the mixing of an O-peak in the simulations is roughly similar to the duration of a single TP cycle. From Fig. 2, we can see that the mixing in the ninth TP for the sequence with M_in = 1.5 M_⊙ occurs in less than 50 000 years.

Fig. 2. Top two panels: Surface luminosity (top) and helium luminosity (bottom) during the ninth thermal pulse of the sequence with Min = 1.5 M⊙. Four stages, labeled as a, b, c, and d, are highlighted, which correspond to the chemical stratifications shown in the four bottom panels. Bottom panels: Formation of the O peak during the thermal pulses. Panels b and c show that the O peak is formed in the same way using both prescriptions, while panel d shows how the additional instabilities considered in MLT♯ lead to the mixing of the O peak formed in the previous TP. We also highlight the position of the eighth and ninth O peaks.

Regarding the effects of the MLT♯ and GNA prescriptions in the cases studied here (M_in = 1 M_⊙,  1.5 M_⊙,  and 3 M_⊙), we can see some differences between them. In the M_in = 1 M_⊙ sequence, the star was characterized by four TP before leaving the AGB. This value did not allow for a complete burning of He in the region where the O-peaks were located and, consequently, these regions remained stable. Had the star prolonged its AGB phase (e.g., with a slower rate of mass loss), the remaining He would have been burned and the mixing induced by the negative chemical gradients would have started. The sequences with M_in = 1.5 M_⊙ and M_in = 3 M_⊙ were characterized by 11 and 16 TP, respectively, more than enough to burn the He tail and mix the carbon-oxygen profiles. Furthermore, as the temperature of the helium remnant region rises with the mass of the star, they need fewer TP to start the process. In the present computations, the M_in = 1.5 M_⊙ sequence started the mixing in the fourth TP and the M_in = 3 M_⊙ sequence burned the remaining helium tail as early as in the second TP. As the burning shell advances, the remaining He tail is burned and mixing erases the O-rich peaks. At any given time, this means that only the O-rich peaks created during the last TP (or last two TPs) will still be present in the interior of the model (see middle panel in Fig. 1). Consequently, after the star departs from the TP-AGB only the later TP-induced O peaks will still be present in the model.

We can make one last observation regarding the M_in = 3 M_⊙ sequence. While in the other two sequences the three chemical profiles have the same overall shape (except for the peaks in MLT) and they have, for example, the same place and height for the maximum in the carbon abundance; however, this is not the case for the more massive star. The GNA profile appears to be shifted towards the exterior. This shift is a direct consequence of the stabilizing effect of the positive chemical gradient at the bottom of the H-rich envelope that weakens the intensity of the third dredge up episodes. Such a stabilizing effect is only taken into account in GNA simulations, where the instability criterion is self-consistently based on Eq. (1). We note that our simulations do not consider any kind of convective boundary mixing at the bottom of the convective envelope. Such nonlocal mixing process is known to favor third dredge up and would erode the chemical transition and remove any stabilizing effect of ∇_μ at the bottom of the H-rich envelope (e.g., Herwig 2000; Straniero et al. 2006; Miller Bertolami 2015, 2016; Wagstaff et al. 2020). In our setup with M_in = 3 M_⊙, this resulted in the third dredge-up starting one TP earlier and also in stronger dredge-up episodes in MLT and MLT♯ sequences with respect to GNA ones.

4. Pulsational properties of PG1159 stellar models

In the previous section, we showed that the choice of different prescriptions for the treatment of stellar convection can lead to significant differences in the final chemical profile of a stellar model. As the intensity of those mixing processes are badly understood, we might wonder whether the predictions from GNA and MLT♯ prescriptions can be tested in real stars. One way to test these chemical structures is by studying GW Vir stars with asteroseismological techniques (see Córsico et al. 2019, and references therein).

As mentioned in the introduction, GW Vir are pulsating pre-WDs and WDs with PG1159 spectral types (Werner & Herwig 2006). The vast majority of PG1159 stars are thought to have formed after a very late thermal pulse (VLTP)⁵. This process takes place once the star has departed from the AGB and is entering the WD cooling stage (see Miller Bertolami 2024, for a review on different evolutionary channels). To produce appropriate stellar models of PG1159 stars, compute their pulsational properties and test the impact of different chemical profiles, we need to continue the evolution of the sequences computed in the previous section through their departure from the AGB and through a VLTP. In this first exploratory study, we have focused on the three M_in = 1.5 M_⊙ sequences computed with the three different convection theories. We have tuned the mass loss at the end of the AGB in a way so that the model departs at the right time for a VLTP to occur. We then followed the evolution through the He-flash and the consequent proton ingestion event. Figure 3 shows the resulting evolution on the HR diagram. Following the VLTP, the star evolves back to the giant region of the HR diagram and it later contracts again to the WD cooling track, now as a H-deficient PG1159 star. At this stage, we took three different snapshots of the stellar structure and computed their adiabatic g-mode periods (crosses in Fig. 3). Those snapshots correspond to log T_eff = 5.1, before the knee in the HR diagram, the point of maximum effective temperature at log T_eff = 5.25 and, finally, log T_eff = 5 after the knee, at the entrance of the WD cooling phase.

Fig. 3. Evolutionary track for the Min = 1.5 M⊙ star. The VLTP is highlighted as well as the three points (marked with crosses), where the pulsational calculations were made.

Adiabatic pulsations of our PG1159 models were computed using an updated version of LP-PUL (Córsico & Althaus 2006). Among other quantities, LP-PUL computes the pulsation periods Π_k, ℓ = 2π/σ_k, ℓ, with k being the radial order and ℓ the harmonic degree of a given normal mode. In this work, we limited our analysis to dipolar g-modes (ℓ = 1). To get the values of periods as precise as possible, we have employed about 5000 mesh-points to describe background stellar models. In pulsating stars, in general, the structure of the g-mode period spectrum is extremely sensitive to the radial profile of the squared value of the Brunt-Väisälä (or buoyancy) frequency (N). In order to accurately compute the value of N in the core of PG1159 stars, we used the treatment described by Tassoul et al. (1990), Brassard et al. (1991), where the Brunt-Väisälä frequency is given by

$$\begin{array}{c}\hfill N^2=\frac{g^2\rho }{P}\frac{{\chi }_T}{{\chi }_{\rho }}({\mathrm{\nabla }}_{\text{ad}}-\mathrm{\nabla }-B).\end{array}$$(15)

In this expression, B is the so-called Ledoux term (Tassoul et al. 1990), expressed as

$$\begin{array}{c}\hfill B=-\frac{1}{{\chi }_T}{\displaystyle \sum _1^{M-1}}{\chi }_{X_i}\frac{dln{X}_i}{dlnP},\end{array}$$(16)

while the terms χ_ρ, χ_T, and χ_Xi are

$$\begin{array}{c}\hfill {\chi }_T=(\frac{dlnP}{dlnT}{)}_{\rho ,\{X_i\}},{\chi }_{X_i}=(\frac{dlnP}{dln{X}_i}{)}_{\rho ,T,\{X_{j\ne i}\}},{\chi }_{\rho }=(\frac{dlnP}{dlnT\rho }{)}_{T,\{X_i\}}.\end{array}$$(17)

Figure 4 shows the internal chemical stratification of our computed PG1159 models for the main nuclear species (upper panels), together with the values of the Ledoux term (middle panels), and the squared Brunt-Väisälä frequency (lower panels). The figure emphasizes the role of the chemical interfaces on the shape of the Brunt-Väisälä frequency profile. In fact, each chemical transition leads to distinctive features in N², which are eventually responsible for the mode-trapping properties of the model (Córsico et al. 2019). In the inner part of the core (q = −log1 − m/M_★ < 0.5), there are several peaks in B (and N²) that are related to changes in the mass fractions of C and O created immediately after the end of the He-core burning stage. In the outermost part of the CO core (0.7 < q < 1.5), we can see the abrupt features shaped during the TP-AGB phase. They are specially pronounced in the case of the model computed under the standard assumptions (MLT + SC, blue lines in Fig. 4). In this case, we can see how the O-peaks created during the thermal pulses lead to the formation of clear spikes in N². These spikes are not present in those models computed with either the GNA or MLT♯ approaches, where these regions are mixed by either RT instabilities or thermohaline processes. We expect these spikes to lead to clear mode-trapping features of g-modes.

Fig. 4. Chemical profiles (top), Ledoux term B (middle), and the logarithm of the squared Bränt-Väisälä frequency (bottom) for a PG1159 model with a mass of M = 0.58 M⊙ (Min = 1.5 M⊙) at three different stages after the VLTP. As in Fig. 1, blue lines correspond to MLT, orange lines to MLT♯, and black lines to GNA.

4.1. Mode trapping

Mode trapping is a phenomenon that occurs when a star has abrupt structural features (e.g., sharp chemical transitions) that cause the amplitudes of the pulsation modes to be significantly different from those that would occur if such features were not present. As a consequence, the period differences between consecutive periods are not constant; rather, they vary significantly, having pronounced minima which recur as the period increases (i.e., the “trapping cycle”), as described by Montgomery (2005), Córsico et al. (2019). In WD stars, the structure of the g-mode period spectrum is sensitive to the spatial profile of the Brunt-Väisälä frequency (Córsico et al. 2019). A clear signature of mode trapping is that the forward period spacing (defined as ΔΠ_k = Π_k + 1 − Π_k, as a function of the pulsation period, Π_k) displays strong deviations from a constant value. Moreover, this quantity has been shown to be sensitive to the chemical composition of the star, where chemical gradients in different regions of the star affect the period spacing at different ranges of periods (Córsico & Althaus 2006). In light of the clear differences in the squared Brunt-Väisälä frequency of the different models (shown in Fig. 4), we would expect models computed with the standard MLT to display larger mode trapping features than those computed under the assumption of either GNA or MLT♯. In Fig. 5, we show ΔΠ_k versus Π_k for PG1159 models where the MLT, MLT♯, and GNA prescriptions have been employed at the three different stages given in Figs. 3 and 4. In Fig. 5, we shows that for periods longer than ∼1700 s, MLT models (in blue) have much larger deviations from the mean period spacing than both GNA and MLT♯ models. This is true for the three evolutionary stages depicted in the figure.

Fig. 5. Forward period spacing for the PG1159 models with masses of M = 0.58 M⊙ (Min = 1.5 M⊙) and the three mixing prescriptions at the three different stages selected after the VLTP. As in the previous figures, blue dots correspond to MLT, orange dots to MLT♯, and black dots to GNA.

To ensure that the broad deviations in the constant period spacing observed in the MLT models are indeed coming from the O-peaks created by the thermal pulses (at 0.7 < q < 1.5, Fig. 4), we performed a similar analysis to the one described in Section 3.4 of Córsico & Althaus (2006). This consisted of artificially suppressing the Ledoux term B in different parts of the star. One of the advantages of using Eq. (15) is that the contribution to N² from changes in the chemical composition is almost completely contained in the Ledoux term B. This fact renders the method particularly useful for inferring the relative weight that each chemical transition region has in the mode-trapping properties of the models. With this in mind, we divided the stellar models into three regions, corresponding to the core (q < 0.7), O-peak region (0.7 < q < 1.45), and the envelope (q > 1.45).

In Fig. 6, we show the resulting period spacing for the model with log T_eff = 5 and the MLT prescription, when B is set to B = 0 outside each predefined region. Thus, in the first panel we show the case when the Ledoux term is only considered for the region of the peaks, the second panel considers only the Ledoux term in the core, and the third only in the envelope. Finally, for comparison, the fourth panel sets B = 0 in the whole stellar model, showing that the deviations from a constant forward period spacing seen in the previous panels and in Fig. 5 are consequences of the non-null Ledoux term. As we can see, Fig. 6 shows clearly that the core chemical transitions mostly affect mode trapping for Π < 1500 s while the O-peaks, created by the thermal pulses, lead to a distinctive trapping feature at longer periods. Envelope chemical transitions lead to mode trapping features over the whole range of plotted periods. In light of the similar chemical structure of the envelopes of MLT, MLT♯ and GNA models, Fig. 6 shows that the much larger trapping features observed in Fig. 5 for Π > 1700 s are indeed caused by the O-peaks left behind by the standard MLT treatment. This could be thought of as a useful seismic tool for deducing whether chemically induced mixing processes are present in the interiors of actual stars.

Fig. 6. Forward period spacing versus periods for the MLT case and log Teff = 5 when the Ledoux B term has been artificially suppressed in different regions and does not remain null in the regions specified by the different subscripts of B. The only exception is where B has been suppressed for the whole model (bottom), labeled as B = 0. We see that, when B is only considered not-null in the regions of peaks of the C-O profile, several periods are affected, with the most important contributions being for the periods from 1700 s and higher (top).

Finally, having identified the most significant contributions for the peaks in the period spacing, we can go back to Fig. 5 and test the changes in the amplitudes of ΔΠ statistically, for Π > 1700 s. Thus, we carried out an estimation for the mean value and the variance of the ΔΠ, which are shown in in Fig. 7. We obtained much higher variances for MLT and similar variances between MLT♯ and GNA.

Fig. 7. Mean value μ and variance σ2 of the forward period spacing using periods Π > 1700 s. The squares represent the stage of log Teff = 5.1, the circles are log Teff = 5.25, and the triangles are log Teff = 5. MLT shows significantly greater variances in all three stages.

4.2. Comparison with GW Vir stars

The computation of the forward period spacing in actual stars requires the observation of normal modes of consecutive radial order, which is not always possible. Therefore, when comparing with real pulsating stars, it is often more useful to rely on the residuals of the periods. Residuals are defined as the difference between the observed periods, Π_i^O, and the periods derived from a mean period spacing, Π_i^fit (Córsico et al. 2007). The mean period spacing is obtained making a linear least squares fit to the period sample. We know from both theory and observations (Córsico et al. 2006, 2021) that only the high luminosity GW Vir stars (i.e.m before the “knee” in the HR diagram, see Fig. 3) show pulsation periods longer than 1000 s. Therefore, given that the O-peaks created during the thermal pulses only affect the longer periods, for this exploratory work we will limited our comparisons to the model with log T_eff = 5.1 (see Fig. 8). We can see that the main differences between the residuals in MLT and MLT♯ and GNA arise for higher orders of k, having MLT significantly higher residuals in that regime.

Fig. 8. Linear period fit (upper panel) and residuals (lower panel) for the log Teff = 5.1 models computed under the three different convective theories.

Finally, we aim to test our models with actual observations of pulsating stars. The relevant period range for our work corresponds to periods longer than 1700 s, therefore, we focused our comparisons with GW Vir stars that exhibit most of their periods in that particular range. We selected the following stars: PG 1159-035 (Oliveira da Rosa et al. 2022), NGC 1501, HS 2324 (Córsico et al. 2021), and NGC 246 (Calcaferro et al. 2024a). We replicated the average period spacing found in the references and then proceeded to estimate and compare the residuals with the ones from our simulations. We show the comparison in Fig. 9. Given theoretical models are not asteroseismic models, the comparison is only qualitative; in other words, they are not the result of period-to-period fits to the observations, but instead they represent snapshots of theoretical models in a similar evolutionary stage to that of the observed stars.

Fig. 9. Residuals from the theoretical models (empty circles) and the observed pulsating stars (red triangles). As in the previous figures, blue corresponds to MLT, orange to MLT♯,  and black to GNA. There is a notable change in scale with respect to the comparison with PG 1159-035 (bottom) due to its larger residuals.

The most notable feature in Fig. 9 is the extremely large residuals found in PG 1159-035 (i.e., up to 10 s). Such a large extent of mode trapping is far greater than anything displayed by our models in any period range. Given that chemical transitions in the envelope would be able to produce such huge trapping features, we suspect that these high residuals values are connected to a different envelope structure than those corresponding to our models. This would be in line with the suggestion by Althaus et al. (2008), based on the measured rates of period changes, that PG 1159-035 harbors an envelope much thinner (in terms of mass) than predicted by standard stellar models. For the other stars, the comparison is inconclusive. On the regime of interest, k > 80, HS 2324 shows, on average, residuals that are in line with those observed in our MLT♯ model, notably smaller than those of our standard MLT computations. The residuals of NGC 1501 and NGC 246 are more in line with the predictions of the MLT models, yet the number of periods in the long period regime is rather small. Further studies, based on proper asteroseismic models of the relevant stars, are needed to reach any conclusion that would favor any of the mixing treatments considered in this work and to confirm, for example, whether thermohaline mixing processes or RT instabilities can take place during the TP-AGB phase.

5. Conclusions

In this work, we derived an extension to standard MLT that incorporates the existence of chemical gradients in the background stratification, which we named MLT♯ (see Appendix A). In our formulation, we included the prediction of the double diffusive convective theory GNA, developed by Grossman (1996). Next, we explored the importance of chemically driven convection and thermohaline processes during the TP-AGB phase and explored their consequences for pulsating GW Vir stars.

Using LPCODE, we computed the evolution of stars with initial masses M_i = 1, 1.5 and 3 M_⊙ during the TP-AGB phase where thermal pulses create layers of enhanced O-rich composition (O-peaks) and inversions of the chemical gradients may develop. For each initial mass, we computed three different TP-AGB evolutionary sequences, with the standard MLT, with MLT♯, and with the GNA. We find that when enough thermal pulses take place on the TP-AGB, negative chemical gradients establish that lead to the development of both chemically driven convection and thermohaline processes. In our simulations, this happens for the more massive sequences in our sample (M_i = 1.5 and 3 M_⊙). Both MLT♯ and GNA predict that the mixing timescales for these instabilities are shorter than the length of the TP-AGB phase itself. Consequently, when these processes are taken into account, the final chemical stratification in the outer regions of the CO cores of these models are smoother than those arising from stellar evolution computations that adopt the standard MLT prescription. On the contrary, our M_i = 1 M_⊙ sequences never develop thermohaline or Rayleigh-Taylor instabilities; in addition, all sequences retain a chemical profile with several O-peaks in the outer core, regardless of the adopted convective theory.

We also explored the possibility of using pulsating H-deficient pre-WD stars (GW Vir) to establish the prevalence of these O-peaks in real stars. In this connection, we computed the evolution of our three M_i = 1.5 M_⊙ sequences through a final post-AGB thermal pulse to construct models corresponding to GW Vir stars. Once they were constructed, we computed their adiabatic pulsational spectrum and studied their period spacing features.

We find that the O-peaks created during the TP-AGB produce mode-trapping features in the period spacing for modes with Π > 1700 s. Moreover, in such period range, the mode trapping features caused by the chemical transitions in the inner CO core are irrelevant. We do find, however, that chemical transitions in the envelope can lead to mode trapping features in the whole range of studied periods. Our models were computed with the standard MLT approach and, therefore, they include strong O-peaks, leading to much larger trapping features at Π > 1700 s than those computed with either MLT♯ or GNA. In light of these results we conclude that luminous GW Vir stars (i.e., those before the knee, as seen in Fig. 3 in the HR diagram offer a good testing ground for the occurrence of these mixing process. With this in mind, we chose GW Vir stars with the longest observed periods (PG 1159-035, NGC 1501, HS 2324, and NGC 246) and performed a residuals analysis of their periods in light of the predictions of our theoretical models. This qualitative comparison was inconclusive, with HS 2324 favoring models with smooth chemical profiles, while NGC 246 and NGC 1501 showed larger deviations from the mean period spacing, typical of our models with strong O-peaks in the outermost part of the core. Interestingly, we find that the prototype of the class, PG 1159-035, shows much larger mode trapping features than those observed in any of our models. We speculate that this could be related to the thinner envelopes (in comparison with canonical post-VLTP models) already hinted by Althaus et al. (2008) in light of the very large rates of period change displayed by this star. Besides being inconclusive, we believe our work shows that GW Vir stars offer a great opportunity to test chemically driven mixing processes on the TP-AGB phase.

Our results show that MLT♯ reproduces the predictions of the GNA theory with very good agreement, but through a significantly simpler and more numerically stable formulation. This makes it particularly suitable for stellar evolution codes where the implementation of GNA would be computationally expensive or unstable. Since MLT♯ modifies standard MLT only by incorporating the chemical gradient into the equations, it retains both the clarity of the physical interpretation and the flexibility of existing MLT-based frameworks. We conclude that MLT♯ is a practical and robust tool for modeling chemically driven convection and thermohaline mixing. This approach can be readily adopted in future evolutionary and pulsational studies where double-diffusive effects are relevant.

Acknowledgments

Part of this work was supported by PIP-2971 from CONICET (Argentina) and by PICT 2020-03316 from Agencia I+D+i (Argentina). We thank the Asociación Argentina de Astronomía for supporting the publication costs of this article.

References

Appendix A: Mathematical derivation of MLT♯ equations

In this section, we derive the equations that rule the convective motions in MLT♯. We will follow an almost identical approach as in Chapter 7 of Kippenhahn et al. (2013), indicating where the differences arise.

Considering the sum of radiative flux F_rad and convective flux, F_con, we obtain the radiative gradient ∇_rad needed to transport the whole flux by radiation,

$$\begin{array}{c}\hfill F_{\text{rad}}+F_{\text{con}}=\frac{4acG{T}^{4}m}{3\kappa P{r}^2}{\mathrm{\nabla }}_{\text{rad}}.\end{array}$$(A.1)

Considering that part of the flux is transported by convection, then F_rad is related to the actual temperature gradient ∇ by

$$\begin{array}{c}\hfill F_{\text{rad}}=\frac{4acG{T}^{4}m}{3\kappa P{r}^2}\mathrm{\nabla }.\end{array}$$(A.2)

On the other hand, the local convective flux of energy given by an element moving with a velocity v and in balance of pressure (DP = 0) and excess temperature DT⁶ is given by

$$\begin{array}{c}\hfill F_{\text{con}}=\rho v{c}_{P}DT.\end{array}$$(A.3)

The element will move over a distance l_m, the mixing length, until it mixes with the surrounding. Since the elements will start their motions at different distances relative to a given place (from 0 to l_m), they will have different velocities v and excess temperatures DT. We can consider an “average” element that has moved l_m/2 and then we have

$$\begin{array}{c}\hfill \frac{\mathit{DT}}{T}=\frac{1}{T}\frac{\partial DT}{\partial r}\frac{l_m}{2}=(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}})\frac{l_m}{2H_P}.\end{array}$$(A.4)

This is the point where the differences between the standard treatment of MLT and MLT♯ arise. Considering the density difference without neglecting the impact of the chemical instabilities, we have

$$\begin{array}{c}\hfill \frac{D\rho }{\rho }=-\delta \frac{\mathit{DT}}{T}+\phi \frac{D\mu }{\mu },\end{array}$$(A.5)

with

$$\begin{array}{c}\hfill \frac{D\mu }{\mu }={\mathrm{\nabla }}_{\mu }\frac{l_m}{2H_P}\end{array}$$(A.6)

and the radial buoyancy force per unit mass is

$$\begin{array}{c}\hfill k_r=-g\frac{D\rho }{\rho }.\end{array}$$(A.7)

Considering that, on average, half of this value acts on its motion, l_m/2, giving

$$\begin{array}{c}\hfill \frac{1}{2}{k}_r\frac{l_m}{2}=g\delta (\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}})\frac{l_m^2}{8H_P}-\phi {\mathrm{\nabla }}_{\mu }\frac{l_m}{4}.\end{array}$$(A.8)

Considering that half of the work goes into the kinetic energy, the average velocity of the elements will be given by

$$\begin{array}{c}\hfill v^2=g\delta (\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu })\frac{l_m^2}{8H_P}.\end{array}$$(A.9)

Assuming that the chemical gradient ∇_μ remains constant over the course of the timescale of the convective motion we need to obtain only the value for ∇ and ∇_e. We will use the following relation, from Kippenhahn et al. (2013):

$$\begin{array}{c}\hfill \frac{{\mathrm{\nabla }}_{\text{e}}-{\mathrm{\nabla }}_{\text{ad}}}{\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}}=\frac{6ac{T}^3}{\kappa {\rho }^2{c}_P{l}_{m}v},\end{array}$$(A.10)

where we assume the convective element to be a sphere of diameter l_m. This means that the mixing length is assumed to be the same as the size of the blob of matter.

Replacing A.9 in the expression above, we can obtain the first of the equations that of MLT♯:

$$\begin{array}{c}\hfill ({\mathrm{\nabla }}_{\text{e}}-{\mathrm{\nabla }}_{\text{ad}})(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }{)}^{1/2}=2U(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}),\end{array}$$(A.11)

with U equally defined as in Section 2.

Next, in order to obtain the second equation, we can replace A.4 and A.9 in A.3, obtaining the following expression for the convective flux:

$$\begin{array}{c}\hfill F_{\text{con}}=\frac{\rho c_P\sqrt{g\delta }}{4\sqrt{2}}{l}_m^{2}T{H}_P^{-3/2}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }{)}^{1/2}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}).\end{array}$$(A.12)

Also, from A.1 and A.2, we have

$$\begin{array}{c}\hfill F_{\text{con}}=\frac{4acG{T}^{4}m}{3\kappa P{r}^2}({\mathrm{\nabla }}_{\text{rad}}-\mathrm{\nabla }).\end{array}$$(A.13)

Finally, using dP/dr = −Gmρ/r² and the definition for H_P, we can combine A.12 and A.13 in order to obtain the second equation of MLT♯, expressed as

$$\begin{array}{c}\hfill (\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}}-\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu }{)}^{1/2}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\text{e}})=\frac{8}{9}U({\mathrm{\nabla }}_{\text{rad}}-\mathrm{\nabla }).\end{array}$$(A.14)

Eqs. A.11 and A.14 can be solved numerically for obtaining ∇ and ∇_e. Replacing the gradients in A.9 will give the convective velocity of the element.

All Figures

	Fig. 1. Chemical abundances of He (dotted lines), C (continuous lines) and O (dashed lines) of stellar models of initial mass Min = 1 M⊙ (left panel), 1.5 M⊙ (middle), and 3 M⊙ (right) at the end of the AGB phase. The blue lines correspond to the case where MLT was used, the orange ones are for MLT♯, and the black ones are for GNA.
In the text

Fig. 2. Top two panels: Surface luminosity (top) and helium luminosity (bottom) during the ninth thermal pulse of the sequence with Min = 1.5 M⊙. Four stages, labeled as a, b, c, and d, are highlighted, which correspond to the chemical stratifications shown in the four bottom panels. Bottom panels: Formation of the O peak during the thermal pulses. Panels b and c show that the O peak is formed in the same way using both prescriptions, while panel d shows how the additional instabilities considered in MLT♯ lead to the mixing of the O peak formed in the previous TP. We also highlight the position of the eighth and ninth O peaks.

In the text

	Fig. 3. Evolutionary track for the Min = 1.5 M⊙ star. The VLTP is highlighted as well as the three points (marked with crosses), where the pulsational calculations were made.
In the text

	Fig. 4. Chemical profiles (top), Ledoux term B (middle), and the logarithm of the squared Bränt-Väisälä frequency (bottom) for a PG1159 model with a mass of M = 0.58 M⊙ (Min = 1.5 M⊙) at three different stages after the VLTP. As in Fig. 1, blue lines correspond to MLT, orange lines to MLT♯, and black lines to GNA.
In the text

	Fig. 5. Forward period spacing for the PG1159 models with masses of M = 0.58 M⊙ (Min = 1.5 M⊙) and the three mixing prescriptions at the three different stages selected after the VLTP. As in the previous figures, blue dots correspond to MLT, orange dots to MLT♯, and black dots to GNA.
In the text

Fig. 6. Forward period spacing versus periods for the MLT case and log Teff = 5 when the Ledoux B term has been artificially suppressed in different regions and does not remain null in the regions specified by the different subscripts of B. The only exception is where B has been suppressed for the whole model (bottom), labeled as B = 0. We see that, when B is only considered not-null in the regions of peaks of the C-O profile, several periods are affected, with the most important contributions being for the periods from 1700 s and higher (top).

In the text

	Fig. 7. Mean value μ and variance σ2 of the forward period spacing using periods Π > 1700 s. The squares represent the stage of log Teff = 5.1, the circles are log Teff = 5.25, and the triangles are log Teff = 5. MLT shows significantly greater variances in all three stages.
In the text

	Fig. 8. Linear period fit (upper panel) and residuals (lower panel) for the log Teff = 5.1 models computed under the three different convective theories.
In the text

	Fig. 9. Residuals from the theoretical models (empty circles) and the observed pulsating stars (red triangles). As in the previous figures, blue corresponds to MLT, orange to MLT♯,  and black to GNA. There is a notable change in scale with respect to the comparison with PG 1159-035 (bottom) due to its larger residuals.
In the text