EDP Sciences
Free Access
Issue
A&A
Volume 501, Number 2, July II 2009
Page(s) 659 - 677
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/200811414
Published online 13 May 2009

The core helium flash revisited

II. Two and three-dimensional hydrodynamic simulations

M. Mocák1,2 - E. Müller1 - A. Weiss1 - K. Kifonidis1

1 - Max-Planck-Institut für Astrophysik, Postfach 1312, 85741 Garching, Germany
2 - Institut d'Astronomie et d'Astrophysique, Université Libre de Bruxelles, CP 226, 1050 Brussels, Belgium

Received 25 November 2008 / Accepted 3 April 2009

Abstract
Context. We study turbulent convection during the core helium flash close to its peak by comparing the results of two and three-dimensional hydrodynamic simulations.
Aims. In a previous study we found that the temporal evolution and the properties of the convection inferred from two-dimensional hydrodynamic studies are similar to those predicted by quasi-hydrostatic stellar evolutionary calculations. However, as vorticity is conserved in axisymmetric flows, two-dimensional simulations of convection are characterized by incorrect dominant spatial scales and exaggerated velocities. Here, we present three-dimensional simulations that eliminate the restrictions and flaws of two-dimensional models and that provide a geometrically unbiased insight into the hydrodynamics of the core helium flash. In particular, we study whether the assumptions and predictions of stellar evolutionary calculations based on the mixing-length theory can be confirmed by hydrodynamic simulations.
Methods. We used a multidimensional Eulerian hydrodynamics code based on state-of-the-art numerical techniques to simulate the evolution of the helium core of a 1.25 $M_{\odot}$ Pop I star.
Results. Our three-dimensional hydrodynamic simulations of the evolution of a star during the peak of the core helium flash do not show any explosive behavior. The convective flow patterns developing in the three-dimensional models are structurally different from those of the corresponding two-dimensional models, and the typical convective velocities are lower than those found in their two-dimensional counterparts. Three-dimensional models also tend to agree more closely with the predictions of mixing length theory. Our hydrodynamic simulations show the turbulent entrainment that leads to a growth of the convection zone on a dynamic time scale. In contrast to mixing length theory, the outer part of the convection zone is characterized by a subadiabatic temperature gradient.

Key words: stars: evolution - hydrodynamics - convection - stars: interiors

1 Introduction

The core helium flash is the most violent event in the life of a star with an initial mass between approximately 0.7 $M_{\odot}$ and 2.2 $M_{\odot}$ (Sweigart & Gross 1978). Electron degeneracy in the helium core at the time of the flash leads to a thermonuclear runaway that at its peak produces large amounts of energy ( ${\sim}10^{10}~L_{\odot}$) within the stellar core over a very short period of time ($\sim $days). The temperature rises up to several 108 Kelvins until the degeneracy of the electron gas is eventually lifted. Nevertheless, the event does not seem to be catastrophic for the star. It only results in a slow expansion of the helium core (typically with a few $\mbox{~m~s${}^{-1}$ }$), since energy transport by turbulent convection, heat conduction, and radiation seems to be able to deliver most of the flash energy quiescently from the stellar interior to the outer stellar layers. However, computations of the flash have a confusing history of predicting either severe explosions (Deupree 1984a; Wickett 1977; Deupree 1984b; Cole et al. 1985; Deupree & Cole 1983; Cole & Deupree 1980; Zimmermann 1970; Edwards 1969; Cole & Deupree 1981; Villere 1976) or quiescent behavior (Lattanzio et al. 2006; Dearborn et al. 2006; Deupree 1996).

Previous two-dimensional (2D) hydrodynamic simulations of turbulent convection within the helium core during the peak of the core helium flash (Mocák et al. 2008) support a quiescent scenario in agreement with stellar evolutionary calculations. However, they also show that the convection zone, powered by helium burning, is characterized by convective velocities that are roughly four times higher than those predicted by mixing-length theory, and the width of the convection zone grows on a dynamical timescale. This may lead to mixing of hydrogen into the helium core as known from one-dimensional simulations of extremely metal-poor stars (Campbell & Lattanzio 2008; Cassisi et al. 2003; Fujimoto et al. 1990; Schlattl et al. 2001; Weiss et al. 2004a; Hollowell et al. 1990).

Table 1:   Some properties of the initial model.

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg21.eps}\hspace*{1cm}
\includegraphics[width=8.5cm,clip]{1414fg22.eps}
\end{figure} Figure 1:

Left panel: temperature (in units of $10^{7}\mbox{~K}$) distribution of the mapped (dashed) and stabilized (solid) initial model displayed with the density (in units of  $10^{5}\mbox{~g~cm${}^{-3}$ }$) and entropy (in $k_{\rm B}$ per baryon) stratification of the stellar evolution model (dotted and dash-dotted), respectively. The numbers indicate the pressure scale height Hp (1/Hp = -|dln p/dr| in 108 cm) and the density scale height $H_\rho $ in brackets (1/$H_\rho $ = -|dln $\rho $/dr| in 108 cm) at positions given by the closest black dots, respectively. Right panel: chemical composition of the initial model. The dotted vertical lines mark the edges of the convection zone.

Open with DEXTER

It is well known that two-dimensional hydrodynamic simulations of turbulence are seriously biased due to the imposed symmetry restrictions. Opposite to three-dimensional (3D) flows, the turbulent kinetic energy increases from small to large scales in 2D simulations, i.e., the energy cascade to smaller length scales characteristic of turbulent flows is not reproduced (Canuto 2000; Bazan & Arnett 1998). Hence, the mean convective velocities, the amount of overshooting, and the size of turbulent structures are too large. Thus, as pointed out already by e.g., Muthsam et al. (1995) and Bazan & Arnett (1998), three-dimensional simulations are required to validate the predictions of two-dimensional simulations.

With increasing computational capabilities, the importance of multidimensional hydrodynamic simulations in stellar evolution studies grows rapidly because they are essentially ``parameter free''. In contrast, canonical one-dimensional stellar evolutionary calculations involve free parameters like the mixing length or the overshooting distance (Böhm-Vitense 1958; Canuto & Mazzitelli 1991), which are chosen in an appropriate manner to fit observational data (Montalbán et al. 2004). Comparison of the results obtained with both approaches is crucial, because it allows one to constrain, validate or disprove the free parameters used in stellar evolutionary calculations, and because it can provide a clue to the applicability limits of the canonical (1D) treatment (Herwig et al. 2006; Meakin & Arnett 2007b; Arnett et al. 2007; Meakin & Arnett 2006; Dearborn et al. 2006; Asida & Arnett 2000; Kercek et al. 1999,1998; Bazan & Arnett 1998).

In the following we present an investigation of the core helium flash by means of two and three-dimensional hydrodynamic simulation using state-of-the-art numerical techniques, a detailed equation of state, and a time-dependent gravitational potential. Our hydrodynamic models are characterized by a convectively unstable layer (the convection zone) embedded in between two stable layers composed of several chemical nuclear species and of a partially degenerate electron gas. Similar systems were studied in the past by many authors e.g., Hurlburt et al. (1994,1986), Muthsam et al. (1995), Singh et al. (1995,1998) and Brummell et al. (2002) assuming, however, that the stellar matter is composed of a single ideal Boltzmann gas. This gives extra relevance to our simulations because they allow us to study, e.g., the dependence of turbulent entrainment and the structure of convective boundary layers on the composition of the stellar gas, and on the composition gradients present in the stellar model.

We introduce the stellar model used as input for our multidimensional simulations in the next section. We then briefly describe our hydrodynamics code in Sect. 3, and discuss and compare the results of the two and three-dimensional simulations in Sect. 4. In particular, we compare the temporal evolution, the kinetic energy density, the power spectra, the thermodynamic fluctuation amplitudes, the stability of the flow structures, and the turbulent entrainment at the convective boundaries of our models. We discuss our results in the context of the predictions of mixing-length theory in Sect. 5. Finally, we present a high resolution 2D simulation covering almost 1.5 days of core helium flash evolution in Sect. 6, partially described by Mocák et al. (2008). We end with a summary of our findings and some conclusions in Sect. 7.

2 Initial model

The initial model was obtained with the stellar evolution code Garstec (Weiss & Schlattl 2000,2007). Some of its properties are listed in Table 1, and the distributions of temperature, density and composition are displayed in Fig. 1. The model possesses a white dwarf like degenerate structure with an off-center temperature maximum resulting from plasma and photo-neutrino cooling. Its central density is ${\sim}7$ $\times $ 105 g cm-3. At the outer edge of the isothermal region in the center of the helium core the temperature jumps up to $T_{\rm max} \sim 1.7$ $\times $ 108 K, and the adjacent convection zone has a superadiabatic temperature gradient. The core is predominantly composed of 4He (mass fraction X(4He) > 0.98), and some small amounts of 1H, 3He, 12C, 13C, 14N, 15N and 16O, respectively. For our hydrodynamic simulations we adopt the abundances of 4He, 12C and 16O as the triple-$\alpha$ reaction dominates the nuclear energy production during the flash. The remaining composition is assumed to be adequately represented by a gas with a mean molecular weight equal to that of 20Ne.

As our multidimensional hydrodynamics code Herakles (see next section) utilizes an Eulerian computational grid, our initial data for the hydrodynamic simulations are obtained by polynomial interpolation from the stellar model which was computed on a Lagrangian grid. The initial hydrodynamics model obtained in this way is not in hydrostatic equilibrium, because the equation of state of Timmes & Swesty (2000) implemented in our hydrodynamics code differs from that of Rogers et al. (1996) implemented in the Garstec code (for a given density and temperature in the initial stellar model, the pressure differs by about 1%). Stabilization of this initial hydrodynamics model resulted in a small decrease of the temperature compared to the stellar model (Fig. 1).

3 Hydrodynamics code

The hydrodynamic simulations were performed with an enhanced version of the grid-based code Herakles (Mocák et al. 2008). The adopted mathematical model implemented in the code consists of the Euler equations in spherical coordinates ($r, \theta$, and $\phi $; see Appendix A) coupled to the source terms describing thermal transport, self-gravity, and nuclear burning. The hydrodynamic equations are integrated using the PPM reconstruction scheme (Colella & Woodward 1984), and a Riemann solver for real gases according to Colella & Glaz (1984). The chemical species are evolved by a set of additional coupled continuity equations (Plewa & Müller 1999). Source terms in the momentum and energy equations due to self-gravity and nuclear burning are treated by means of operator splitting at the end of every integration step. The gravitational potential is approximated by a one-dimensional Newtonian potential which is obtained from the spherically averaged mass distribution. The stiff nuclear reaction network is integrated using the semi-implicit Bader-Deuflhard method, (Press et al. 1992; Bader & Deuflhard 1983) which allows for large time integration steps.

In Herakles, a program cycle consists of two hydrodynamic time steps and proceeds as follows (in 2D simulations step (3) is omitted):

1.
The hydrodynamic equations are integrated in r-direction (r-sweep) for one time step including the effects of heat conduction. The time averaged gravitational forces are computed, and the momentum and the total energy are updated to account for the respective source terms. Subsequently, the equation of state is called to update the thermodynamic state due to the change of the total energy.

2.
Step (1) is repeated in $\theta $-direction ($\theta $-sweep).

3.
Step (1) is repeated in $\phi $-direction ($\phi $-sweep).

4.
The nuclear reaction network is solved in all zones with significant nuclear burning (i.e., where T > 108 K), and then the equation of state is called to update the pressure and the temperature.

5.
In the subsequent time step, the order of Step (1)-(3) is reversed to guarantee second-order accuracy of the time integration, and Step (4) is repeated with the updated state quantities.

6.
The size of the time step for the next cycle is determined.
The velocities in the convection zone, even at the peak of the core helium flash, are subsonic corresponding to Mach numbers $M \sim 0.01$ (Mocák et al. 2008). In this regime, the applicability of Riemann solver based methods, like PPM, has been questioned (Almgren et al. 2006; Schneider et al. 1999; Turkel 1999). However, a recent study by Meakin & Arnett (2007a) based on a direct comparison of anelastic (Kuhlen et al. 2003) and fully compressible simulations computed with PPM shows that at Mach numbers around 10-2, PPM can capture the properties of convective flows well.

Table 2:   Some properties of the two and three-dimensional simulations: number of grid points in r (Nr), $\theta $ ( $N_{\theta }$) and, $\phi $ ($N_{\phi }$) direction, spatial resolution in r ($\Delta r$ in 106 cm), $\theta $ ( $\Delta \theta $), and $\phi $ ( $\Delta \phi $) direction, characteristic velocity $v_{\rm c}$ (in 10 $^{6}\mbox{~cm~s${}^{-1}$ }$) of the flow, expansion velocity at temperature maximum  $v_{\rm exp}$ (in  $\mbox{~cm~s${}^{-1}$ }$), typical convective turnover timescale  $\tau _{\rm cnv}$ (in s) and maximum evolutionary time  $t_{\rm max}$ (in s), respectively.

4 Hydrodynamic simulations

We performed two two-dimensional and one three-dimensional simulation whose properties are summarized in Table 2[*]. The simulations hefl.2d.a and hefl.3d cover 9500 s and 6000 s of evolutionary time during the peak of the core helium flash, respectively, and were computed on a computational grid spanning a wedge of  $120\mbox{$^\circ$ }$ in angular directions centered at the equator. The grid had a radial resolution of $\Delta r = 5.55$ $\times $ 106 cm, and an angular resolution of $\Delta \theta = \Delta \phi = 1.33\mbox{$^\circ$ }$. The rather wide angular grid appeared to be necessary for the three-dimensional simulations, due to the size of the largest vortices ($\sim $40$^\circ$) found in previous two-dimensional simulations (Mocák et al. 2008).

The 2D model hefl.2d.b was evolved for roughly 36 h (130 000 s) on a grid covering the full  $180\mbox{$^\circ$ }$ angle in $\theta $ direction. It was simulated at almost twice the resolution as the other models. The properties of model hefl.2d.b and its evolution are discussed in Sect. 6, and to some extent also in Sect. 4. We used reflective boundary conditions in every coordinate direction and for all models, except model hefl.2d.a. For this model, it turned out to be necessary to impose periodic boundary conditions in angular direction, because reflective boundaries together with the 120 $^\circ$ wedge size affect the large scale convective flow adversely leading to higher convection velocities.

We are limited with our explicit hydrodynamics code by the CFL condition which is most restrictive near the center of the spherical grid. Therefore, we cut off the inner part of the grid at a radius of 2 $\times $ 108 cm that is still sufficiently far away from the radius of the temperature inversion at $r \sim 5$ $\times $ 108 cm. To trigger convection, we imposed a random flow field with a maximum (absolute) velocity of  $10 \mbox{~cm~s${}^{-1}$ }$, and a random density perturbation $\Delta \rho / \rho \le 10^{-2}$. Imposing some explicit non-radial perturbations is necessary, because a spherically symmetric model evolved with Herakles on a grid in spherical coordinates will remain spherically symmetric forever. This is different from the study of Asida & Arnett (2000), who did not perturb the initial model by an artificial random flow, as instabilities were growing from round-off errors. The different perturbation techniques seem not to influence the final thermally relaxed steady state (Meakin & Arnett 2007b).

Because thermal transport of energy by conduction and radiation is roughly seven orders of magnitude smaller than the convective energy flux, it has been neglected in our simulations. Most of the liberated nuclear energy is carried away by convection. All computed models are non-rotating, because rotation seems not to play an important role during the core helium flash (Lattanzio et al. 2006).

4.1 Important timescales

A major disadvantage of multidimensional hydrodynamic simulations of stars today is related to the fact, that stars evolve on a nuclear timescale which is typically of the order of 109 years and even with the best supercomputers nowadays, we are not able to cover more than a day of hydrodynamic stellar evolution (Mocák et al. 2008; Herwig et al. 2006; Meakin & Arnett 2007b). Nevertheless, hydrodynamic simulations are an important tool for studying the evolution of stars. Let us first define three major timescales, which one can associate with nuclear burning and hydrodynamic flow.

  • Nuclear timescale $\tau_{\rm brn}$: an e-folding time, $\tau_{\rm brn} = T / \dot{T} \approx c_V T/ \dot{\epsilon}_{\rm nuc}$, where  $\dot{\epsilon}_{\rm nuc}$ is the energy release rate of the nuclear processes, T is the temperature, and cV is the specific heat.

  • Sound crossing timescale $\tau_{\rm hyd}$: a time, during which a region of a radius $\delta r$ reacts to any pressure imbalance. Such a reaction occurs with the local sound speed $c_{\rm s}$ and therefore the sound crossing time scale is $\tau_{\rm hyd} = \delta r / c_{\rm s}$ ( $c_{\rm S}\sim 1.7$ $\times $ $10^{8}\mbox{~cm~s${}^{-1}$ }$ in center of our models).

  • Convective turnover timescale $\tau _{\rm cnv}$: a time, which a convective element needs to cross the convection zone of width  $R_{\rm cnv}$ with a velocity  $v_{\rm cnv}$, i.e., $\tau_{\rm cnv} = 2~R_{\rm cnv} / v_{\rm cnv}$.
When the convective turnover timescale $\tau _{\rm cnv}$ in a stellar convection zone becomes comparable or smaller than the nuclear burning timescale  $\tau_{\rm brn}$, the border between quiescent quasi-hydrostatic evolution and an explosion of a star is narrow. At this stage convection is losing its ability to transport as much energy as it is produced by nuclear reactions, and the convective layers begin to expand rapidly. Whenever such a situation occurs in stellar evolution, the convective flow should be studied by hydrodynamic simulations, as they are able to correctly follow the reaction of the flow to any pressure imbalance occurring on the sound crossing timescales  $\tau_{\rm hyd}$. According to stellar evolution calculations, nuclear timescale of our model M (Table 1) is $\tau_{\rm brn} \sim 20$ days which is much longer compared to the convective turnover timescale $\tau _{\rm cnv}$ $\sim $ 2000 s. This implies that convection is able to transport the released nuclear energy away from burning regions, and the star remains in quasi-hydrostatic equilibrium. Our hydrodynamic simulations (Table 2) give timescales that are similar roughly within a factor of 2-4. However, one has to keep in mind that the quantitative conclusions based on the two-dimensional simulations which have the highest resolutions we performed are anyway seriously biased due to the assumption of axial symmetry. On the other hand, the 3D simulation hefl.3d is not fully resolved yet and the timescales mentioned above will change for the same simulation with higher resolution.

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg41.eps}
\end{figure} Figure 2:

Temporal evolution of the total nuclear energy production rate E of models hefl.3d (solid) and hefl.2d.a (dotted), respectively. $L_{\odot }$ is the solar luminosity.

Open with DEXTER

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg42.eps}
\end{figure} Figure 3:

Temporal evolution of the total kinetic energy $E_{\rm K}$ of models hefl.3d (solid) and hefl.2d.a (dotted), respectively.

Open with DEXTER

4.2 Temporal evolution

The 3D model hefl.3d and the corresponding 2D model hefl.2d.a evolve initially (t < 1200 s) similarly. Convection sets in after roughly 1000 s, and hot rising bubbles appear in the region where helium burns in a thin shell ($r \sim 5$ $\times $ 108 cm). After another ${\sim}200$ s, the bubbles cover the complete height of the convection zone (R $\sim $ 4.7 $\times $ 108 cm). The flow eventually approaches a quasi-steady state after roughly 2000 s consisting of several upstreams (or plumes). The total nuclear energy production (Fig. 2) during the quasi-steady state is about 10% higher in the 3D model due to the strong dependence of the triple-$\alpha$ reaction rate on temperature which is slightly higher in the 3D model. Its evolution is correlated with the total kinetic energy (Fig. 3) mainly during the onset of convection when the value of the total kinetic energy rises from 1039 erg to 1045 erg, i.e., hot gas near the position of temperature maximum begins to move and carries the liberated nuclear energy off the burning region, thereby inhibiting a thermonuclear runaway.

 \begin{figure}
\par\includegraphics[width=8.5cm,clip]{1414fg43.eps}\par\includegraphics[width=8.5cm,clip]{1414fg44.eps}
\end{figure} Figure 4:

Temporal evolution of a log10 of the angular averaged kinetic energy density (in $\mbox{~erg }\mbox{~g}^{-1}$) of models hefl.3d ( upper) and hefl.2d.a ( lower), respectively.

Open with DEXTER

The kinetic energy density (Fig. 4) shows small fluctuations in the fully evolved convection zone, and is by an order of magnitude larger in model hefl.2d.a than in model hefl.3d. This is in agreement with other studies, as it is well know that two-dimensional turbulence is more intensive (Muthsam et al. 1995). The convection zone in the Fig. 4 is characterized by the strip of large kinetic energy density and the adjacent convectively stable layers by waves induced by convection. The region of high kinetic energy density appearing at ${\sim}3000$ s in the bottom layer in model hefl.3d (see top panel of Fig. 4) is an artifact caused by the proximity of the reflective inner radial boundary (at 2 $\times $ 108 cm) of the computational grid.

We observe buoyancy driven gravity waves in the convectively stable layers (Figs. 411) (Hurlburt et al. 1994,1986; Meakin & Arnett 2007b; Zahn 1991), but we do not discuss these waves any further here, because their properties are likely biased by the reflective boundaries (Asida & Arnett 2000). We only point out that the differences in amplitude seen in Fig. 4 are physical, i.e., gravity waves have a lower amplitude in 3D than in 2D, although the energy carried by them is likely similar (Kiraga et al. 1999).

Table 3:   Root mean square fluctuation amplitudes within the convection zone.

4.3 Structure of the convective flow

Fully evolved convection (t > 2000 s) in the 3D model hefl.3d differs significantly from that in the corresponding 2D model hefl.2d.a. The convective flow is dominated in the 2D model by vortices having (angular) diameters ranging from 30$^\circ$ to 50$^\circ$ (Fig. 5), and an aspect ratio of close to one. The vortices are qualitatively similar to those found in other two-dimensional simulations (Hurlburt et al. 1994,1986; Bazan & Arnett 1998; Porter & Woodward 1994). This vortex structure of 2D turbulence is quite typical, and arises from the self-organization of the flow (Fornberg 1977; McWilliams 1984).

The convective flow in the 3D model hefl.3d consists of column-shaped plumes (Figs. 511), and contrary to the 2D model hefl.2d.a, it does not show any dominant angular mode. The typical angular size of turbulent features ranges from 10$^\circ$ to 30$^\circ$ (Fig. 5). The power spectra of angular velocity fluctuations show that turbulent elements have an almost time-independent characteristic angular size of 30$^\circ$-50$^\circ$ in case of the 2D model, while the spectra computed for the 3D model change with time and exhibit no dominant angular mode.

We find turbulent flow features across the whole convection zone resulting from the interaction of convective up and downflows. Close to the edges of the convective zone we observe the smallest turbulent flow features that form when compact turbulent plumes are decelerated and break-up (Brummell et al. 2002). Shear instabilities likely play an important role in the development of the turbulent flow as well (Cattaneo et al. 1991).

 \begin{figure}
\par\includegraphics[width=5.6cm,clip]{1414f45a.eps}\hspace*{0.5c...
...ps}\hspace*{1.5cm}
\includegraphics[width=7.9cm,clip]{1414f45e.eps}
\end{figure} Figure 5:

Upper panels: snapshots taken at ${\sim }4754$ s showing contour plots of the absolute value of the velocity (in units of $10^{6} \mbox{~cm~s${}^{-1}$ }$) for the 2D model hefl.2d.a ( left), and for the 3D model hefl.3d in a meridional plane of azimuth angle $\phi = 50\mbox{$^\circ$ }$( middle) and $\phi = 70\mbox{$^\circ$ }$ ( right), respectively. Lower panels: normalized power spectra of angular fluctuations in the absolute velocity as a function of radius for the 2D model helf.2d.a ( left) averaged over time from 2100 s to 9500 s, and for the 3D model hefl.3d ( right) averaged over time from 2250 s to 6000 s, and azimuthal angle. The dashed vertical lines mark the edges of the convective zone of the initial model according to the Schwarzschild criterion.

Open with DEXTER

Tables 3 and 4 provide time-averaged (over a time period of about 2500 s) root mean square fluctuation amplitudes of various variables of the convective flow inside the convection zone[*] and near its edges, respectively for models hefl.3d and hefl.2d.a. The temperature and density fluctuation amplitudes are 30-40% larger in the 2D model than in the 3D one. This is expected, because vortices are stable in 2D flows but decay in 3D ones (Fig. 5). The fluctuation amplitudes in the composition (4He and 12C) are larger by 10-30% in the 3D model which is a result of more ``broken'' and hence more non-uniform mixing of chemical elements than in the 2D one. The temperature and density fluctuation amplitudes are a factor of 2-3 larger near the inner edge than near the outer edge[*] of the convective layer. At both edges are the fluctuation amplitudes by a factor of 2-4 larger in 2D than in 3D models. The fluctuation amplitudes in the composition (4He and 12C) in both models differ close to the convective boundaries as well, by up to a factor of 2.

Table 4:   Root mean square fluctuation amplitudes at the edges of the convection zone (see Table 3).

 \begin{figure}
\par\includegraphics[width=5.9cm,clip]{1414f46a.eps}\hspace*{0.2c...
...s}\hspace*{0.2cm}
\includegraphics[width=5.9cm,clip]{1414f46i.eps}
\end{figure} Figure 6:

Cuts through the 3D model hefl.3d at t = 4815 s showing the angular variation of temperature (in units of 108 K; left panels), 12C mass fraction (in units of 10-3; middle panels), and radial velocity (in units of $10^{5} \mbox{~cm~s${}^{-1}$ }$; right panels), respectively, at three different radii: r1 = 4.8 $\times $ 108 cm (temperature maximum; top), r2 = 6.5 $\times $ 108 cm (center of the convection zone; middle), and r3 = 9.3 $\times $ 108 cm ( top of convection zone, bottom). The black lines in the right panels mark the boundaries between positive and negative radial velocities i.e., upflows and downflows.

Open with DEXTER

Cuts through the 3D model hefl.3d at t = 4815 s showing the angular variation of temperature, 12C mass fraction, and radial velocity at three different radii (Fig. 6) demonstrate that the helium core at the peak of the core helium flash is a very turbulent environment at all heights of the convection zone. The bottom of the convection zone contains hot filaments of gas where the temperature exceeds that of the environment by about 1%. The filaments contain ashes from helium burning, i.e., 12C and 16O, and they move across the whole bottom of the convection zone in a random way. The filaments are correlated with upflows, as the hot gas of burned matter is forced by buoyancy to rise towards the top of the convection zone.

The apparent turbulent nature of the convective flow indicated by our simulations implies that the treatment of mixing in stars as a diffusive process may lead to inaccurate or even incorrect results. Convective flows are rather advective, as pointed out by Woodward et al. (2008).

4.4 Stability of flow structures

To analyze the size and the stability of the vortices we introduce an auto-correlation function of the radial velocity

\begin{displaymath}%
A_{1}(r_{0};r) = \frac{\langle v_{r}(r_{0})~ v_{r}(r)\rangl...
...2}_{\Omega,t}~
\langle v_{r}(r)^{2}\rangle^{1/2}_{\Omega,t}}
\end{displaymath} (1)

that measures the radial extent of flow patterns at a given radius r0, or equivalently the radial size of vortices. The notation  $\langle \rangle_{\Omega,t}$ indicates averaging over angles and time. A second auto-correlation function

\begin{displaymath}%
A_{2}(r_{0};t-t_{0}) = \frac{\langle v_{r}(r_{0},t_{0})~
v...
...\Omega} ~
\langle v_{r}(r_{0},t)^{2} \rangle^{1/2}_{\Omega}}
\end{displaymath} (2)

provides a measure of the lifetime of flow patterns at radius r0, beyond time t0. Here, we average only over angles. Both correlation functions have the properties $-1 \le A_{1,2} \le 1$, and A1(r0;r0) = 1 and A2(r0;t0-t0) = 1, respectively. They are similar to the autocorrelation function used by Chan et al. (1982), and Chan & Sofia (1986).

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg47.eps}
\end{figure} Figure 7:

Auto-correlation function A1(r0;r) measuring the radial extent of flow patterns (see Eq. (1)) at different radii r0 (4.8 $\times $ 108 cm, 5.9 $\times $ 108 cm, 7 $\times $ 108 cm, 8.1 $\times $ 108 cm, and 9.25 $\times $ 108 cm) and $t \sim 4000$ s measuring the radial extent of flow patterns for the 2D model hefl.2d.a (dotted lines) and the 3D hefl.3d (solid lines), respectively.

Open with DEXTER

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg48.eps}
\end{figure} Figure 8:

Auto-correlation function A2(r0;t-t0) measuring the lifetime of flow patterns (see Eq. (2)) at two different epochs for the 3D model hefl.3d (solid: $t_{0} \sim 2260$ s, dashed: $t_{0} \sim 2900$ s), and for the 2D model hefl.2d.a (crosses, diamonds), respectively. The radius r0 is 7.6 $\times $ 108 cm.

Open with DEXTER

Figure 7 displays the radial auto-correlation for models hefl.2d.a and hefl.3d and confirms the extension of the convective flow across the whole convective region as determined by the Schwarzschild criterion in the initial stellar model. The broad plateaus with $A_{1} \approx 1$ corresponding to the 2D model hefl.2d.a bear eveidence of the axial symetry imposed in the two-dimensional case which leads to pronouced circular vortices. In the three-dimensional case, the distributions of A1 tend to differ from unity at nearly all radii.

Figure 8 shows the temporal auto-correlations with two typical results for different t0. The three-dimensional model always shows the typical behaviour of a decrease of the function value to 0.5 within 200-250 s. From this we conclude that the flow pattern fluctuates always in the same way. This is different in the two-dimensional model, where A2 can keep high values for a much longer time implying rather persistent structure (the vortices) of the convective flow.

4.5 Turbulent entrainment and the width of the convection zone

Convection may induce mixing in convectively stable layers adjacent to convectively unstable regions. Following Meakin & Arnett (2007b), we prefer to call this process turbulent entrainment (or mixing), a term also known in oceanography, see e.g., Fernando (1991). The commonly used term (convective) overshooting accounts only for localized ascending or descending plumes crossing the edge of the convection zone. If the filling factor of these plumes or their crossing frequency is high, they can change the entropy in convectively stable regions surrounding convection zones, a process that is known as penetration (Brummell et al. 2002). Turbulent entrainment accounts for both overshooting and penetration.

 \begin{figure}
\par\includegraphics[width=7.4cm,clip]{1414f49a.eps}\vspace*{5mm}...
....eps}\vspace*{5mm}
\includegraphics[width=7.4cm,clip]{1414f49c.eps}
\end{figure} Figure 9:

Carbon mass fraction X(12C) ( top), entropy S ( middle), and entropy gradient $\nabla S$ ( bottom) as a function of radius near the outer edge of the convection zone of model hefl.3d at three different epochs: t1 = 2000 s (dashed), t2 = 4000 s (dash-dotted), and t3 = 6000 (solid). In addition, the initial X(12C) profile is shown in the top panel (dotted).

Open with DEXTER

Contrary to Hurlburt et al. (1994) we determine the depth of the entrainment neither by the radius where the kinetic energy carried into the stable layers is zero, nor by the radius where the kinetic energy has dropped to a certain fraction of its maximum value (Brummell et al. 2002). We find both conditions insufficient, because the kinetic energy flux becomes ``zero'' much faster than the convective flux (Fig. 12), which is another possible indicator of the depth of the entrainment (see next Sect. 4.6). Instead, we rather use the 12C mass fraction, as it is low outside the convection zone during the flash (X(12C) < 2 $\times $ 10-3), and as it can rise there only due to turbulent entrainment.

In this study, we use the condition X(12C) = 2 $\times $ 10-3 to define the edges of the convection zone. Due to the turbulent entrainment, these edges are pushed towards the stellar center and surface (Fig. 9). Hence, the width of the convection zone increases on the dynamic timescale, which is in contradiction with the predictions of one-dimensional hydrostatic stellar modeling, where the width of the convection zone is determined by the local Schwarzschild or Ledoux criterion. However, the criterion for the width of the convection zone cannot be a local one due to turbulent entrainment caused by convection.

The speed, at which the radius of the outer edge of the convection zone increases with time due to entrainment, is estimated for models hefl.2d.a and hefl.3d to be at most 14  $\mbox{~m~s${}^{-1}$ }$. The radius of the inner edge of the convection zone changes at a much smaller rate (Meakin & Arnett 2007b; Bazan & Arnett 1998), as the region interior to the convection zone is more stable against convection and has a higher density than the region exterior to the convection zone (Singh et al. 1995). The entrainment at the bottom of the convection zone also leads to a heating of the cool interior layers (Cole et al. 1985; Deupree & Cole 1983). This seems to be a robust feature of convection zones driven by nuclear burning, and is observed in other studies too (e.g., Asida & Arnett 2000).

The region just interior to the convection zone shows less entrainment (Meakin & Arnett 2007b; Bazan & Arnett 1998), as the square of the Brunt-Väisälä frequency is almost ten times larger there than that in the region just outside the outer edge of the convection zone. The Brunt-Väisälä frequency is a good stability indicator since it is related to the behavior of convective elements within a convection zone, a fact also pointed out by Hurlburt et al. (1994). The Brunt-Väisälä frequency can be written as (Meakin & Arnett 2007b):

\begin{displaymath}%
N^2 = -g \left( \frac{\partial \ln \rho}{\partial ~ r} -
\frac{\partial \ln \rho}{\partial ~ r}
\bigg\vert _{s} \right),
\end{displaymath} (3)

where g, $\rho $, and r are the gravitational acceleration, the density and the radius, respectively. A layer is convectively stable, if N2 > 0, and unstable otherwise[*] (Kippenhahn & Weigert 1990).

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414f410.eps}
\end{figure} Figure 10:

Square of the Brunt-Väisälä buoyancy frequency as a function of radius for the 3D model hefl.3d at t = 4638 s.

Open with DEXTER

The Brunt-Väisälä frequency differs in the 2D and 3D simulations slightly, and has on average a very small negative value inside the convection zone: N2 = -4.4 $\times $ 10-6 rad $^2\mbox{~s}^{-2}$ in the 2D model hefl.2d.a, and N2 = -1.2 $\times $ 10-6 rad $^2\mbox{~s}^{-2}$ in the 3D model hefl.3d. Just outside (inside) the inner (outer) edge of the convection zone of the 2D model hefl.2d.a we find N2 = 0.580 rad $^2\mbox{~s}^{-2}$ (0.053 rad $^2\mbox{~s}^{-2}$), and N2 = 0.583 rad $^2\mbox{~s}^{-2}$ (0.052 rad $^2\mbox{~s}^{-2}$) for the 3D model hefl.3d, respectively (Fig. 10). In the 2D model hefl.2d.b, these frequencies are higher by about a factor of two.

 \begin{figure}
\par\mbox{\includegraphics[width=8.5cm,clip]{1414iso1.eps}\includ...
...ps}\hspace{2.5cm}
\includegraphics[width=6cm,clip]{1414iso4.eps} }\end{figure} Figure 11:

Different views of isosurfaces of the velocity field for the 3D model hefl.3d at t = 3000 s. The blue isosurface corresponds to radial downflows with vr = -6 $\times $ $10^{5} \mbox{~cm~s${}^{-1}$ }$, and the yellow and brown isosurfaces show radial upflows with vr = 6 $\times $ $10^{5} \mbox{~cm~s${}^{-1}$ }$, and 1 $\times $ $10^{4}\mbox{~cm~s${}^{-1}$ }$ (gravity waves), respectively. The edge sizes of the box are 1.2 $\times $ 109 cm and 2.4 $\times $ 109 cm, respectively. The yellow-greenish sphere in the bottom right panel marks the top of the convection zone according to the Schwarzschild criterion.

Open with DEXTER

Due to its high stability, the radius of the bottom edge of the convection zone did not change during the time covered by our simulations hefl.2d.a and hefl.3d, except for an initial jump inwards over one radial grid zone, when it was touched by convective downflows for the first time. However, entrainment may move the edge further towards the stellar center later in the evolution. The entrainment rate (i.e., the velocity at which the convective boundary moves) is lower by a factor of $\sim $6 at the bottom of the convection zone in our 2D simulation hefl.2d.b (see Sect. 6) as compared to that at the outer boundary ( ${\sim}14\mbox{~m~s${}^{-1}$ }$; see above). Similar behavior was also observed in 3D simulations of oxygen burning shell (Meakin & Arnett 2007b). It implies that the entrainment rate at the inner edge of the convection zone is ${\sim}2.3\mbox{~m~s${}^{-1}$ }$ in our core helium flash simulations. The corresponding change in radius is only $\sim $1.4 (2.2) $\times $ 106 cm or about a third of the width of a radial zone during the time covered by the simulations hefl.3d (hefl.2d.a), and hence too small to be seen. As these estimates are resolution dependent, the values presented should be considered as order of magnitude estimates, only.

The entrainment is more efficient in the 2D model hefl.2d.a than in the 3D model hefl.3d. In 2D, the observed convective flow structures are large and fast rotating vortices that due to the imposed axisymmetry are actually tori (Bazan & Arnett 1998). They have a high filling factor near the edge of the convection zone where they overshoot or penetrate (Brummell et al. 2002). 3D structures crossing the edge of the convection zones are smaller (localized) plumes with a lower filling factor and smaller velocities.

-1mm

We studied convective stability in detail only for the layer above the convection zone, where the gas is only weakly degenerate, and we are thus able to compare our results with those of similar systems simulated in 3D at high Reynolds numbers ( ${\sim}10^{4}$) by Brummell et al. (2002). According to these authors, a stable layer allows for more ``overshooting'', if its entropy gradient is smaller. We found that the turbulent entrainment lowers the entropy gradient at the outer edge of the convection zone in our simulations (Fig. 9). Hence, the stability of the exterior stable layer decreases with time, allowing for more turbulent entrainment.

Contrary to Brummell et al. (2002), who studied only single fluid flow, our simulations involve a mixture of different fluids of different composition. This complicates the above argumentation, as in a multi-fluid flow a shallow entropy gradient does not necessarily imply that the layer is less stable against turbulent entrainment. We plan to address this issue in more detail elsewhere. We have not analyzed the stability of the region below the inner edge of the convection zone, because it is highly degenerate and appears to be very stable. No significant entrainment of 12C was observed there during the entire simulation of models hefl.2d.a and hefl.3d, respectively. Simulations covering longer periods of time are therefore needed (see Sect. 6). The analysis of the energy fluxes inside the convection zone and near the outer edge of the convection zone, which will be discussed in the next subsection, will provide more information about the phenomenon of the turbulent entrainment.

4.6 Energy fluxes

Energy fluxes are a useful tool for understanding convective flows (Hurlburt et al. 1994; Brummell et al. 2002; Hurlburt et al. 1986; Muthsam et al. 1995; Meakin & Arnett 2007b; Chan & Sofia 1986). They allow one to discriminate the energy transport due to different processes and mechanisms (e.g., due to convection, heat conduction, etc.). We thus analyzed various energy fluxes and source terms that are defined in Appendix B. All energy fluxes and source terms (see, Fig. 12) are averaged over several convective turnover timescales, as they change considerably with time due to the appearance of plumes (Muthsam et al. 1995).

The convective energy flux is mainly positive as heat is mostly transported upwards by convection. It is larger in the 2D model hefl.2d.a, as in 2D the convective flow structures are more laminar and ordered, and thus experience less dissipation. The smaller value of $F_{\rm C}$ in the 3D model hefl.3d is a result of a less ordered flow throughout the whole convection zone. This is in agreement with Table 3 which shows that the fluctuations in temperature and density are smaller in the 3D model.

A kinetic energy flux arises from deviations from the mean convective flow (i.e., mainly from the upflow-downflow asymmetry). Typically it is largest in the most turbulent regions of the flow, where on top of the kinetic flux due to the up- and down-flow asymmetry, there is also a significant contribution due to the asymmetry of localized turbulent (convective) elements. Directly connected to this is the offset between the maxima of the kinetic and the convective flux (see Fig. 12), which reflects the fact that the convective flow decays more efficiently when the flow becomes strongly turbulent.

 \begin{figure}
\par\mbox{\includegraphics[width=8.8cm]{1414flx1.eps}\includegrap...
...=8.8cm]{1414flx3.eps}\includegraphics[width=8.8cm]{1414flx4.eps} }\end{figure} Figure 12:

Various energy fluxes and source terms as a function of radius averaged (from 2000 s to 6000 s) over about four convective turnover timescales. Panel  a) shows the convective flux $F_{\rm C}$ of the 2D model hefl.2d.a (solid-red) and the 3D model hefl.3d (dashed-black) together with the kinetic energy flux $F_{\rm K}$ in the 2D (dash-dotted-red) and 3D (dash-dot-dotted-black) model, respectively. The dotted vertical lines mark the edges of the convection zone in the initial model according to the Schwarzschild criterion. Panel  b) gives the source terms due to the work done by buoyancy forces PA (dash-dot-dotted black) and due to volume changes PP (dashed black) in the 3D model hefl.3d, and in the 2D model hefl.2d.a (dash-dotted-red, solid red) respectively. The vertical lines enclose the nuclear burning zone (T > 108 K). Panels  c) and  d) show an enlarged view of the energy fluxes and source terms displayed in panels  a) and  b) near the outer edge of the convection zone.

Open with DEXTER

The work done by buoyancy forces is positive in the whole region of dominant nuclear burning (see Fig. 12) indicating either less dense and hot gas moving upwards or more dense and cooler gas moving downwards. A negative value of the work done by buoyancy forces implies the opposite situation, i.e., less dense and hot gas flowing downward or denser and cooler gas flowing upward. The latter is known as buoyancy braking leading to a deceleration of the flow motion (Brummell et al. 2002), and to the unusual situation that hot matter tends to sink, and cool matter is likely to rise.

The gas should expand while rising up through the convection zone, i.e., the work done by buoyancy (PA) and by volume changes (PP) should always be anti-correlated. The anti-correlation is clearly seen only in the 2D model hefl.2d.a, whereas in the 3D model hefl.3d the quantities are on average anti-correlated only in the central region of the convection zone. At the inner edge of the convection zone, buoyancy drives gas upwards which is on average simultaneously compressed, probably by the broad downflows. At the outer edge where buoyancy braking occurs, the gas on average expands. Hence, it must expand faster than the upflows cool and are being compressed (a situation observed for the 2D model hefl.2d.a).

All the fluxes discussed here agree qualitatively well with those of our previous high-resolution 2D simulations (Mocák et al. 2008). They are also qualitatively very similar to those of the high-resolution 3D simulations of Brummell et al. (2002) who investigated a stratified model with a convectively stable region located on top of a convectively unstable region both consisting of an ideal gas with a very high Reynolds number ( ${\sim}10^{4}$). The angular and time averaged radial distributions of the kinetic and convective fluxes seem to be robust in the convectively unstable region, as our 3D results are qualitatively similar to those obtained in several other 3D studies (Hurlburt et al. 1994; Brummell et al. 2002; Meakin & Arnett 2007b).

The outer part of the convection zone, where buoyancy braking occurs, resembles the overshooting region due to active penetration of plumes described by Hurlburt et al. (1994,1986) and Brummell et al. (2002). Note, however, that this region is convectively stable at the beginning of their simulations, i.e., buoyancy braking takes already place inside the convection zone in our models.

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414ve2d.eps}\hspace*{1cm}
\includegraphics[width=8cm,clip]{1414ve3d.eps}
\end{figure} Figure 13:

Radial distributions of the time (from 2000 s to 4000 s) and angle-averaged velocity components (vr, $v_\theta $, $v_\phi $) and velocity modulus ( $v_{\rm abs}$) for the 2D model hefl.2d.a ( left), and the 3D model hefl.3d ( right), respectively. The panels also show the velocity predicted by the mixing-length theory ( $v_{\rm mlt}$).

Open with DEXTER

The distribution of the kinetic flux (Fig. 12) exhibits structural differences in the convection zone between 2D and 3D flows. The typical evolved 2D flow contains well defined vortices (Fig. 5) whose central parts never interact with the region of dominant nuclear burning. This results in a reduced kinetic energy flux at $r \sim 5.5$ $\times $ 108 cm, as this region corresponds to the central region of the vortices, which do not experience any strong radial motion. On the other hand, the distribution of the kinetic energy flux in the 3D model hefl.3da is rather smooth as a result of the column-shaped flow structures (Fig. 11).

The convective flux changes sign in stable layers since the downflows or upflows penetrating into the stable zones are suddenly too hot or cold compared to the surrounding gas. The penetration continues until the momentum of the convective elements is used up or diffusion smooths out the perturbations, and the convective flux approaches zero (Brummell et al. 2002). This fingerprint of turbulent entrainment is clearly present at both convective boundaries of our models (Fig. 12). It is less pronounced at the inner convective boundary in the 3D model hefl.3d, because the simulation cover only a short period of time and the boundary was not weakened enough for the entrainment to be sufficiently visible by a negative convective flux.

The convective energy flux is relatively strong even in regions where the kinetic energy flux is almost zero. Therefore, the ``zero'' kinetic flux criterion (Brummell et al. 2002; Hurlburt et al. 1994) seems to be a bad indicator of entrainment which may extend well beyond the location where the kinetic flux becomes small.

In fact, what is happening at the convective boundaries is an exchange between the potential energy of the stratification given by the buoyancy jump ${\rm d}b = N^2 ~{\rm d}r$ and the kinetic energy of the turbulence. Turbulence looses its kinetic energy by doing work against gravity, which leads to a reduction of the buoyancy jump, and hence to stability weakening of the stable boundary layer to the effects of the turbulent entrainment. The buoyancy jump ${\rm d}b$ is a direct measure of the stability of the boundary layer. To mix gas into the boundary layer the buoyancy must be reduced. This is accomplished through the buoyancy flux $q = P_A / \rho$, which is related to the temporal variation of the buoyancy jump by ${\rm d}b/{\rm d}t = - {\rm div}(q)$, where PA and $\rho $ are the sink/source term of the kinetic energy due to buoyancy forces and the density, respectively.

The convective flux can directly be related to the buoyancy flux that is a function of PA by a linear relation $F_{\rm C} = F_{\rm C} (q)$ described in more detail by Meakin & Arnett (2007b). Figure 12 shows that this agrees well with what we observe in our simulations. It also supports our previous conclusion that entrainment is well indicated by the convective flux and the related buoyancy flux which via the equation ${\rm d}b/{\rm d}t = - {\rm div}(q)$ leads to the decrease of the buoyancy jump in the stable layer. This in turn reduces the convective stability of that layer.

The properties of the entrainment at the outer convective boundary differs in models hefl.2d.a and hefl.3d (see Fig. 12). In the 2D model entrainment is reaching deeper into the stable layer due to a more active convection zone with higher typical velocities (Fig. 13) than in the 3D model. The radial distribution of the work done by buoyancy PA is qualitatively similar in both the 2D and 3D models, i.e., it is negative indicating buoyancy braking.

Nevertheless, the work done by gas compression or expansion PP is different. In the 2D model hefl.2d.a the gas on average is compressed (PP is positive), while in the 3D model hefl.3d the gas is expanding at the boundary. This again confirms that 2D and 3D convective flows are qualitatively different.

4.7 The flow within the convection zone

The velocities in the 3D model hefl.3d match those predicted by mixing-length theory better than in the 2D model hefl.2d.a, where they are about a factor of $\sim $3 larger. This behavior was also observed in other hydrodynamic simulations of convective flows; see e.g., Muthsam et al. (1995); Meakin & Arnett (2007b). The radial velocities in the regions above and below the convection zone are smaller than the angular velocities in both models, which is a typical feature of gravity waves (Asida & Arnett 2000).

In the 3D model hefl.3d, the flow in the convectively stable layer beneath the convection zone exhibits some numerical artefacts due to the proximity of the inner grid boundary (see Fig. 13). The radial distributions of the time and angle-averaged components of the velocity field and of the density fluctuations show pronounced maxima at ${\sim}3$ $\times $ 108 cm and a sharp cut-off at r = 2.5 $\times $ 108 cm (Figs. 1314). The sharp cut-off is caused by the artificial damping we had to apply to the velocity field in the innermost grid zones to prevent numerical instabilities from spreading limitless to the convection zone.

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414f414.eps}
\end{figure} Figure 14:

Radial distributions of the time (from 2800 s to 4800 s) and angle-averaged density fluctuations in the 3D model hefl.3d (solid red line). The panels also show the angular variation of the respective quantity at a given radius (gray shaded region). The dotted vertical lines mark the edges of the convection zone as determined by the Schwarzschild criterion.

Open with DEXTER

Although, the flow velocities in the 3D model match those predicted by mixing-length theory very well, one should keep in mind that with increasing resolution the flow velocities will likely increase due to the reduced numerical viscosity. This trend is confirmed by the velocities obtained for the high-resolution 2D model hefl.2d.b that are higher than in the low-resolution model hefl.2d.a (Table 2).

Near both edges of the convective zone there are large narrow peaks visible in the radial distributions of the time and angle-averaged density fluctuations (Fig. 14). These peaks are not caused by compression or expansion, but they are a result of the density discontinuities at the edges of the convection zone (Arnett et al. 2007; Meakin & Arnett 2007a), because any angle-dependent radial perturbation will cause large angular variations of density at these discontinuities. The same holds also for the temperature discontinuities at the convective borders.

The temperature fluctuations within the convection zone are rather uniformly distributed, but they are more intense near the outer edge of the convection zone, where they are only weakly correlated with the radial velocity (Fig. 6). At the top of the convection zone the emerging rising plumes are embedded in an environment which sinks down (Fig. 6, bottom panels). This situation is similar to the sinking down-drafts with upwelling centers found in simulations by Nordlund & Dravins (1990) and Cattaneo et al. (1991).

4.8 Upflow-downflow asymmetry

The 2D and 3D simulations share the common property of an upflow-downflow asymmetry (Brummell et al. 2002; Hurlburt et al. 1994; Muthsam et al. 1995). The downflows cover a much larger volume fraction of the convection zone than the upflows (Fig. 15). The filling factor of the downflows increases with decreasing depth (Rast et al. 1993; Meakin & Arnett 2007b) across most of the convection zone, and the downflows are more dominant in the 3D model hefl.3d than in the 2D models. In the model hefl.3d (hefl.2d.a, hefl.2d.b) we find that the absolute velocities are about 40% (15%, 4%) higher in the upflows than in the downflows. Hence, the downflows in the convection zone of a star at the peak of the core helium flash are slower and broader than the faster and narrower upflows.

5 Mixing length theory and simulations

Mixing length theory (or MLT) commonly used for treating convection in stellar evolutionary calculations relies on assumptions and parameters that are often chosen based on convenient ad-hoc arguments about the convective flow, like e.g., the value of the mixing length, the amount of upflow-downflow symmetry or the position where, within the convection zone, convective elements start to rise (Kippenhahn & Weigert 1990; Weiss et al. 2004b).

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414f415.eps}
\end{figure} Figure 15:

Fractional volume occupied by upflow streams $A_{{\rm up}}$ averaged over $\sim $2000 s as a function of radius for the 3D model hefl.3d (solid), the 2D low-resolution model hefl.2d.a (dashed-dotted), and the 2D high-resolution model hefl.2d.b (dotted), respectively.

Open with DEXTER

 \begin{figure}
\par\includegraphics[width=7.8cm,clip]{1414adg1.eps}\vspace*{0.5cm}
\includegraphics[width=7.8cm,clip]{1414adg2.eps}
\end{figure} Figure 16:

Upper panel: radial distributions of the adiabatic temperature gradient  $\nabla _{\rm ad}$ (dotted) and of the temperature gradients of model hefl.3d (solid), respectively. The latter distribution is a linear fit to the gradients averaged over angle and over the first 200 s of the evolution of the model. The gray shaded region marks the convection zone CVZ. Lower panel: same as above, but showing the radial distributions in the evolved convection zone averaged over roughly 3000 s of evolutionary time.

Open with DEXTER

MLT assumes that the temperature of a convective element (blob) is the same as that of the ambient medium surrounding it when it starts to rise. However, as a blob will not rise until it is hotter than the surroundings, this MLT assumption is contradictory. MLT further assumes that once the blobs begin to rise they carry their surplus of heat lossless over a distance given by the mixing length before they release it to the surrounding gas instantaneously at the end of their path. These assumptions are also not fulfilled in general. Our simulations show that convective elements typically start their rise deep inside the star from the region of dominant nuclear burning where they are accelerated by buoyant forces. The assumptions of MLT that convective blobs form and begin their motion at different depths of the convection zone, and that the average convective blob propagates a distance equal to half of the assumed mixing length before dissolving with the surrounding gas (Kippenhahn & Weigert 1990), therefore do not hold. MLT finally also assumes a perfect correlation between the thermodynamic variables and the velocity of the flow in a convection zone. However, the results of our simulations falsify this assumption (Fig. 6).

According to MLT, the temperature fluctuations in a convection zone are directly proportional to the mixing length and to the deviation of the temperature gradient of the model $\nabla_{{\rm sim}} =
({\rm d}~\ln T / {\rm d}~\ln p)_{{\rm sim}}$ from the adiabatic one $\nabla_{{\rm ad}} = ({\rm d}~\ln T/{\rm d}~\ln p)_{{\rm ad}}$:

\begin{displaymath}%
\frac{T^{'}}{T} = (\nabla_{{\rm sim}} -
\nabla_{{\rm ad}}) \frac{1}{H_{p}} \frac{\Lambda}{2}
\end{displaymath} (4)

where $\Lambda$ is the mixing length, Hp the pressure scale height, T' the absolute value of the temperature deviation from the mean (horizontally averaged) temperature T, and p the pressure, respectively.

 \begin{figure}
\par\includegraphics[width=5.7cm,clip]{1414vl1.eps}\hspace*{5mm}
...
...s}\hspace*{5mm}
\includegraphics[width=5.7cm,clip]{1414vl3.eps}
\par\end{figure} Figure 17:

Snapshots of the spatial distribution of the velocity modulus |v| (in units of $10^{6} \mbox{~cm~s${}^{-1}$ }$) for the 2D model hefl.2d.b at 24 000 s ( left), 60 000 s ( middle), and 120 000 s ( right), respectively.

Open with DEXTER

Since T'/T and $\nabla_{{\rm sim}} - \nabla_{{\rm ad}}$ can directly be obtained from our simulations, we attempted to test MLT in a qualitative manner. Our simulations show that in the outer part of the convection zone, i.e., in the region where the buoyancy force is getting smaller (Fig. 12, panel b), the temperature gradient of the models  $\nabla_{{\rm sim}}$ becomes subadiabatic (see lower panel of Fig. 16). Equation (4) which was derived for the adiabatic rise of convective bubbles would then imply that the temperature of convective elements should be lower than that of the surrounding gas (hence no convection) in the outer part of the convection zone, or that the value of $\Lambda$ should be negative. However, convective elements do not rise adiabatically in our hydrodynamic simulations and the subadiabatic gradient means only that the convective elements start to cool faster than their surroundings. It does not imply necessarily that the elements are already cooler than the surrounding gas which would prevent the gas from being convectively active. Note that initially, the temperature gradient is superadiabatic in the whole convection zone (see upper panel of Fig. 16), because the stellar evolutionary model used as initial input for our simulations is computed under the assumptions of the MLT.

6 Long-term evolution

Two-dimensional simulations are biased due to the imposed symmetry restriction which tend to overestimate the activity in the convection zone, but qualitative similarities with geometrically unconstrained 3D models exists. Moreover, many phenomena that we observe in 2D models happen also in the 3D ones, just slower. Hence, we think it is justified to explore the long-term evolution (i.e., covering a few tens of hours instead of a few hours) of our initial core helium flash models by performing cost-effective 2D instead of very costly 3D simulations.

In the following we describe the long-term evolution of the 2D model hefl.2d.b whose early evolution, covering 8 h, was discussed in detail by Mocák et al. (2008). The model is characterized by a very dynamic flow involving typical convective velocities of ${\sim}1.8$ $\times $ $10^{6} \mbox{~cm~s${}^{-1}$ }$. Our long-term hydrodynamic simulation of this model covering 36 h (see Figs. 1718) has revealed that the global and angle-averaged maximum temperatures continue to rise at the initial rate of  $40\mbox{~K~s${}^{-1}$ }$ that is 60% lower than the rate predicted by stellar evolutionary calculations. As a consequence, the typical convective velocities increase by about 50% and reach a level of ${\sim}2.8$ $\times $ $10^{6} \mbox{~cm~s${}^{-1}$ }$ at the end of the simulation (Fig. 17).

Hydrodynamic simulations of convection driven by nuclear burning covering several convective turnover timescales show a rapid growth of the convection zone due to the turbulent entrainment (Meakin & Arnett 2007b). An analysis of our simulations based on a tracing of the radial position of the convective boundaries (defined by the condition X(12C) = 2 $\times $ 10-3; see Sect. 4.5), shows a similar behavior.

 \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg62.eps}
\end{figure} Figure 18:

Temporal evolution of the horizontally averaged temperature maximum  $\langle T \rangle_{\rm max}$ (solid), and of the global temperature maximum  $T_{\rm max}$ (solid thin) in the long-term 2D model hefl.2d.b. The dashed line correspond to the temporal evolution of the maximum temperature in the stellar evolutionary calculation.

Open with DEXTER

 \begin{figure}
\par\mbox{\includegraphics[width=8.5cm]{1414f63a.eps}\includegrap...
...8.5cm]{1414f63c.eps}\includegraphics[width=8.5cm]{1414f63d.eps} }\end{figure} Figure 19:

Upper panels: angular averaged 12C mass fraction as a function of radius near the inner ( left) and outer edge ( right) of the convection zone in the long-term 2D model hefl.2d.b at $t = 60~000\mbox{~s}$ (dashed) and $t = 120~000\mbox{~s}$ (dash-dotted), respectively, with temperature stratification (thick) at $t = 120~000\mbox{~s}$. The vertical dotted lines mark the boundaries of the convection zone at t = 0 s. Lower panels: temporal evolution of the position of the inner ( left) and outer ( right) edge of the convection zone in model hefl.2d.b, respectively.

Open with DEXTER

Turbulent motion near the upper edge of the convection zone pumps material into the convectively stable layer at an entrainment rate of roughly 14 $\mbox{~m~s${}^{-1}$ }$ without any significant slowdown over the whole duration of the simulation (Fig. 19) covering $\sim $130 000 s (or more than 250 convective turnover timescales). The entrainment rate at the inner convective boundary is about a factor of six smaller ($\sim $2.3 $\mbox{~m~s${}^{-1}$ }$) and increasing during the second half of the simulation ( t > 60 000; see Fig. 19). These entrainment rates have to be considered as upper limits because of the imposed axisymmetry which leads to exaggerated convective velocities and large filling factors for the penetrating plumes. The turbulent entrainment causes a growth of the convection zone on a dynamic timescale, in agreement with the hydrodynamic models of oxygen shell burning (Meakin & Arnett 2007b).

Both interfaces at the edges of the convection zone remain sharp during the whole length of the simulation (Fig. 19). The entrainment is correlated with a decrease of the temperature at the outer edge of the convection zone due to the decreasing entropy (Fig. 9). At the inner edge of the convection zone the entrainment leads to heating of the cold region at r < 4.72 $\times $ 108 cm that is cooled by neutrinos (Fig. 19). Contrary to the finding of Asida & Arnett (2000), the heating does not penetrate deeper into the star than the mixing of 12C and the other nuclear ashes.

Due to the growth of the convection zone and due to nuclear burning the mean 4He, 12C and 16O mass fractions change in the convection zone at rates listed in Table 5. The mean value of the 12C mass fraction in the convection zone decreases at a rate of $R_i \sim -7.4$ $\times $ $10^{-9}~{\rm s}^{-1}$ until $t \sim 40~000$ s to a value of X(12C) = 5.2 $\times $ 10-3. Then it begins to increase again at roughly the same (absolute) rate $R_f \sim +7.71$ $\times $ $10^{-9}~{\rm s}^{-1}$. The 12C mass fraction decreases at the beginning, because the volume of the convection zone grows initially almost discontinuously due to the sudden start of the entrainment. Hence, nuclear reactions are for a start unable to produce enough carbon to compensate for the volume increase. At $t \sim 110~000$ s the 12C mass fraction has risen again to its initial value of Xi (12C) $\sim $ 5.5 $\times $ 10-3, and at the end of the simulation at t = 130 000 s the Xf (12C) $\sim $ 5.9 $\times $ 10-3, a value that is 7% higher than the initial one.

Table 5:   Rate of change of the mean mass fractions of 4He, 12C, and 16O in the 2D model hefl.2d.b within the first 40 000 s (Ri; in units of $10^{-9}~{\rm s}^{-1}$), and within the time interval $40~000~{\rm s} \le t \le 130~000~{\rm s}$ (Rf; in units of $10^{-9}~{\rm s}^{-1}$) with the initial Xi and final Xf mass fractions (in units of 10-3), respectively.

The mean 16O mass fraction shows a similar trend as that of 12C as its production depends directly on the 12C mass fraction. The mass fraction of 4He rises within the first 40 000 s because convection is dredging up fresh 4He from the convectively stable layers. Later in the evolution the mass fraction of 4He decreases, as it is being constantly burned.

7 Summary

We performed, analyzed and compared 2D axisymmetric and 3D hydrodynamic simulations of the core helium flash. In agreement with our previous study of 2D hydrodynamic models of the core helium flash we find that the core helium flash in three dimensions neither rips the star apart, nor that it significantly alters its structure.

The evolved convection of the 3D models differs qualitatively from that of the axisymmetric ones. The typical convective structure in the 2D simulations is a vortex with a diameter roughly equal to the width of convection zone, whereas the 3D structures are smaller in extent and have a plume-like shape. The typical convective velocities are much higher in the 2D models than in the 3D ones. In the latter models the convective velocities tend to fit those predicted by the mixing length theory better. Both 2D and 3D models are characterized by an upflow-downflow asymmetry, where downflows occupy larger volume fraction but the upflows are faster.

The upper part of the evolved convection zone is characterized by a subadiabatic temperature gradient, where buoyancy braking takes places, i.e., rising convective plumes start to slow down in this region and eventually descend back into deeper layers of the star. Convection should not exist in that mass layer according to mixing length theory.

Our hydrodynamic simulations show the presence of turbulent entrainment at both the inner and outer edge of the convection zone, which results in a growth of the convection zone on a dynamic timescale. While entrainment occurs at an almost constant speed at the outer boundary of the convection zone, it tends to accelerate at the inner one. The entrainment rates are higher at the outer edge than at the inner edge of the convection zone, as the latter is more stable against entrainment. In essence, these findings show that the extent of the inner convection zone powered by helium burning during the flash predicted by the canonical stellar evolution theory is not correct.

The fast growth of the convection zone due to the entrainment has some potentially interesting implications. As the entrainment is not considered in canonical stellar evolutionary calculations, stars evolving towards the core helium flash may never reach a state as the one given by our initial stellar model. Hence, this may influence the growth of the convection zone observed in our hydrodynamic simulations, as the thermodynamic conditions at the edges of the convection zone may differ.

If a rapid growth of the convection zone indeed occurs, the main core helium flash studied here will never be followed by subsequent mini-flashes, as convection will lift the electron degeneracy in the helium core in less than one month. In addition, the helium core will likely experience an injection of hydrogen from the stellar envelope within 10 days and undergo a violent nuclear burning phase powered by the CNO cycle. However, the growth of the convection zone within the core that is simulated in our models does not have to continue until it will reach the outer convection zone extending up to the surface of the star. Hence, mixing of nuclear ashes to the stellar atmosphere (or a dredge-up) does not necessarily take place. But a fast dynamic growth of the inner convection zone will lead to a change of the composition of the stellar core, and consequently of the luminosities of low-mass stars on the horizontal branch.

The possibility of a dredge-up during the core helium flash explains perhaps early R-stars (Izzard et al. 2007) that are carbon enriched and show evidence of s-process elements. Their production could be triggered by hydrogen injection to the helium core. Mixing of these s-process elements together with carbon from the helium and CNO burning layers to the stellar atmosphere may produce signatures typical for the early R-stars.

We find significant differences between the properties and the evolution of 2D and 3D models having a convection zone powered by (semi-)degenerate helium burning. However, as our 3D models are likely not yet fully converged and as they cover only a relatively short period of evolutionary time, long-term 3D simulations using a higher grid resolution are needed to obtain a better and more reliable understanding of the hydrodynamics of the core helium flash. Long-term 3D simulations are especially important, as it is unclear, whether the upper boundary of the helium powered convection zone will be able to reach or even penetrate to the hydrogen envelope, hence a secondary CNO flash and the possible mixing of ash from the helium layer to the upper atmosphere may not occur.

Acknowledgements
The simulations were performed at the Leibniz-Rechenzentrum of the Bavarian Academy of Sciences & Humanities on the SGI Altix 4700 system. The authors want to thank Frank Timmes for some of his publicly available Fortran subroutines which we used in the Herakles code, and Kurt Achatz whose unpublished hydrodynamic simulations of the core helium flash performed during his diploma work in 1995 inspired this work. We are indebted to Casey Meakin for several enlightening discussions and helpful comments. The authors are grateful to an anonymous referee for his/her comments which helped to improve this manuscript.

Appendix A: Hydrodynamic equations in spherical polar coordinates

The hydrodynamic equations of a non-viscous multi-component reactive gas subject to a gravitational potential $\Phi$ and having a heat conductivity K are given in spherical polar coordinates ( $r, \theta, \phi$) by

\begin{displaymath}%
\partial_t \rho
\!+\! \frac{1}{r^2} \partial_{r} (r^2 \rh...
...frac{1}{r\sin\theta} \partial_{\phi} (\rho v_{\phi})
\!=\! 0
\end{displaymath} (A.1)


    $\displaystyle \partial_t (\rho v_{r})
+ \frac{1}{r^2} \partial_{r} (r^2 \rho v_...
...2)
+ \frac{1}{r\sin\theta} \partial_{\theta}(\sin\theta~ \rho v_{r} v_{\theta})$  
    $\displaystyle \qquad + \frac{1}{r\sin\theta} \partial_{\phi} (\rho v_{r} v_{\ph...
...}^2}{r}
- \frac{\rho v_{\phi}^2}{r}
+ \partial_{r} p
=
- \rho \partial_{r} \Phi$ (A.2a)
    $\displaystyle \partial_t (\rho v_{\theta})
+ \frac{1}{r^2} \partial_{r} (r^2\rh...
..._{\theta}^2)
+ \frac{1}{r\sin\theta} \partial_{\phi} (\rho v_{\theta} v_{\phi})$  
    $\displaystyle \qquad + \frac{\rho v_{\theta} v_{r}}{r}
- \frac{\rho v_{\phi}^2 ...
...ta}
+ \frac{1}{r} \partial_{\theta} p =
- \frac{\rho}{r} \partial_{\theta} \Phi$ (A.2b)
    $\displaystyle \partial_t (\rho v_{\phi})
+ \frac{1}{r^2} \partial_{r} (r^2\rho ...
...
v_{\theta} v_{\phi})
+ \frac{1}{r\sin\theta} \partial_{\phi} (\rho v_{\phi}^2)$  
    $\displaystyle \qquad + \frac{\rho v_{\phi} v_{r}}{r}
+ \frac{\rho v_{\phi} v_{\...
...sin\theta} \partial_{\phi} p
=
- \frac{\rho}{r \sin\theta} \partial_{\phi} \Phi$ (A.2c)


                                  $\displaystyle \partial_t (\rho e)
+ \frac{1}{r^2} \partial_{r} \left\{ r^2 [v_{r}
(\rho e + p) - K \partial_{r} T ]
\right\}$  
    $\displaystyle \qquad + \frac{1}{r \sin\theta} \partial_{\theta}
\left\{ \sin\theta~ [v_{\theta} (\rho e + p) -
\frac{K}{r} \partial_{\theta} T]
\right\}$  
    $\displaystyle \qquad + \frac{1}{r \sin\theta} \partial_{\phi}
\left\{ v_{\phi} (\rho e + p) -
\frac{K}{r \sin\theta} \partial_{\phi} \Phi
\right\}=$  
    $\displaystyle - \rho \left( v_{r}\partial_{r} \Phi
+ \frac{v_{\theta}}{r} \part...
...c{v_{\phi}}{r \sin\theta} \partial_{\phi} \Phi
\right)
+ \rho \dot{\varepsilon}$ (A.3)


    $\displaystyle \partial_{t} (\rho X_{k})
+ \frac{1}{r^{2}} \partial_{r}(r^{2} \r...
...
+ \frac{1}{r \sin\theta} \partial_{\theta}
(\sin\theta~ \rho X_{k} v_{\theta})$  
    $\displaystyle \qquad + \frac{1}{r \sin\theta} \partial_{\phi} (\rho X_{k} v_{\phi})
= \rho \dot{X}_{k}, \quad k = 1 \dots N_{\rm nuc}$ (A.4)

where $\rho $, vr, $v_\theta $, $v_\phi $, p, e, T, $\dot{\varepsilon}$, Xk, and $\dot{X}_k$ are the density, the radial velocity, the $\theta $-velocity, the rotation velocity, the pressure, the total specific energy, the temperature, the energy generation rate per mass due to reactions, the mass fraction of species k, and the change of this mass fraction due to reactions, respectively. $N_{\rm nuc}$ is the number of species the gas is composed of.

Appendix B: Energy fluxes

The various contributions to the total energy flux (Hurlburt et al. 1986; Achatz 1995) can be obtained by first integrating the hydrodynamic energy equation given in Appendix A over the angular coordinates $\theta $ and $\phi $. Then, one decomposes both the specific enthalpy $\varepsilon + p/\rho$ (where  $\varepsilon$ is the specific thermal energy) and the specific kinetic energy vivi/2 into a horizontal mean and a perturbation, $f \equiv \overline f+f'$, and obtains[*]

\begin{displaymath}%
\partial_tE + \partial_{r}(F_C + F_K + F_R + F_E)=0
\end{displaymath} (B.1)

where


 
                                 $\displaystyle E = \int_V \rho e~{\rm d}V$ (B.2)
    $\displaystyle F_{\rm C} = \oint v_{r}\rho\cdot \left(\varepsilon+\frac{p}{\rho}\right)'
~r^2\mbox{d}\Omega$ (B.3)
    $\displaystyle F_{\rm K} = \oint v_{r}\rho\cdot \left(\frac12 v_iv_i\right)' ~r^2\mbox{d}\Omega, \quad i=1,2,3$ (B.4)
    $\displaystyle F_{\rm R} = -\oint K\partial_{r}T ~r^2\mbox{d}\Omega$ (B.5)
    $\displaystyle F_{\rm E} = 4\pi r^2\overline{v_{r}\rho}\cdot
\left(~\overline{\varepsilon+\frac p\rho}+\overline{\frac12v_iv_i}
+\Phi~\right).$ (B.6)

Here, the gravitational potential $\Phi$ is assumed to be constant for simplicity. The sum of the various flux terms Fi give the total energy transported per unit time across a sphere of radius r by different physical processes. One has the convective (or enthalpy) flux, $F_{\rm C}$, the flux of kinetic energy, $F_{\rm K}$, and the flux due to heat conduction and radiation, $F_{\rm R}$. Finally, $F_{\rm E}$, includes all terms causing a spherical mass flow, i.e., the model's expansion or contraction, while $F_{\rm C}$ and $F_{\rm K}$ rest on deviations from this mean energy flow (vortices). The latter are the major contributors to the heat transport by convection.

In a similar way one can also formulate a conservation equation for the mean horizontal kinetic energy that provides further insight into the effects of convective motions. Using the other hydrodynamic equations (Eqs. (A.1) to (A.2c)), and the relation $\partial_t (\rho v_iv_i/2) = v_i \partial_t(\rho v_i) - v_i v_i \partial_t\rho/2$, one finds

\begin{displaymath}%
\partial_tE_{\rm K} + \partial_{r}(F_{\rm K} + F_P + F_{\rm E,K})
= P_A + P_P + P_{\rm E,K}
\end{displaymath} (B.7)

With $F_{\rm K}$ as introduced above, one obtains


                                                 $\displaystyle E_{\rm K} = \int_V \frac\rho2 v_i v_i ~{\rm d}V$ (B.8)
    $\displaystyle F_P = -\oint v_{r} p' ~r^2\mbox{d}\Omega$ (B.9)
    $\displaystyle F_{\rm E,K} = 4\pi r^2\overline{v_{r}\rho}\cdot
\left(~\overline{\frac p\rho+\frac{v_i v_i}2}~\right)$ (B.10)
    $\displaystyle P_A = -\oint v_{r}\rho'\partial_{r}\Phi ~r^2\mbox{d}\Omega$ (B.11)
    $\displaystyle P_P = \oint p'\partial_i v_i ~r^2\mbox{d}\Omega$ (B.12)
    $\displaystyle P_{\rm E,K} = 4\pi r^2\cdot
\left(~\overline p~\overline{\partial_i v_i}-
\overline v_{r}\overline\rho~\partial_{r}\Phi~\right), \quad i=1,2,3$ (B.13)

where the Pi are source or sink terms of the kinetic energy. They are separated into the effect of buoyancy forces (PA), and the work due to density fluctuations (PP, volume changes). By analyzing the various Pi one can determine what causes the braking or acceleration of the convective flow. The acoustic flux, FP, describes the vertical transport of density fluctuations. $F_{\rm E,K}$ and  $P_{\rm E,K}$ describe the effect of expansion (volume work, and work against the gravitational potential), similar to FE in Eq. (B.6).

References

 

Footnotes

...[*]
We simulated another 3D model at lower resolution than model hefl.3d (using a wedge of 90 $\times $ 80 $\times $ 80 zones centered at the equator with $\Delta r = 11.1$ $\times $ 106 cm and $\Delta \theta = \Delta \phi = 1.5\mbox{$^\circ$ }$), but we do not discuss this model here any further, because of its larger numerical viscosity and its almost 50% smaller convective velocities.
... zone[*]
Extension of the convection zone in radius r is according to the Schwarzschild criterion: 4.72 $\times $ $10^{8}~{\rm cm} \le r \le 9.2$ $\times $ 108 cm.
... edge[*]
Inner and outer edge of the convection zone is located at r = 4.72 $\times $ 108 cm and r = 9.2 $\times $ 108 cm, respectively.
... otherwise[*]
The Brunt-Väisälä frequency is related to the bulk Richardson number known from oceanography, which is a more direct measure of the stability of the edges of a convection zone in the presence of a turbulent flow (Meakin & Arnett 2007b).
... obtains[*]
Note that the equations given in Mocák et al. (2008) that correspond to Eqs. (B.2), (B.6) and (B.8) contain some small typographical errors.

All Tables

Table 1:   Some properties of the initial model.

Table 2:   Some properties of the two and three-dimensional simulations: number of grid points in r (Nr), $\theta $ ( $N_{\theta }$) and, $\phi $ ($N_{\phi }$) direction, spatial resolution in r ($\Delta r$ in 106 cm), $\theta $ ( $\Delta \theta $), and $\phi $ ( $\Delta \phi $) direction, characteristic velocity $v_{\rm c}$ (in 10 $^{6}\mbox{~cm~s${}^{-1}$ }$) of the flow, expansion velocity at temperature maximum  $v_{\rm exp}$ (in  $\mbox{~cm~s${}^{-1}$ }$), typical convective turnover timescale  $\tau _{\rm cnv}$ (in s) and maximum evolutionary time  $t_{\rm max}$ (in s), respectively.

Table 3:   Root mean square fluctuation amplitudes within the convection zone.

Table 4:   Root mean square fluctuation amplitudes at the edges of the convection zone (see Table 3).

Table 5:   Rate of change of the mean mass fractions of 4He, 12C, and 16O in the 2D model hefl.2d.b within the first 40 000 s (Ri; in units of $10^{-9}~{\rm s}^{-1}$), and within the time interval $40~000~{\rm s} \le t \le 130~000~{\rm s}$ (Rf; in units of $10^{-9}~{\rm s}^{-1}$) with the initial Xi and final Xf mass fractions (in units of 10-3), respectively.

All Figures

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg21.eps}\hspace*{1cm}
\includegraphics[width=8.5cm,clip]{1414fg22.eps}
\end{figure} Figure 1:

Left panel: temperature (in units of $10^{7}\mbox{~K}$) distribution of the mapped (dashed) and stabilized (solid) initial model displayed with the density (in units of  $10^{5}\mbox{~g~cm${}^{-3}$ }$) and entropy (in $k_{\rm B}$ per baryon) stratification of the stellar evolution model (dotted and dash-dotted), respectively. The numbers indicate the pressure scale height Hp (1/Hp = -|dln p/dr| in 108 cm) and the density scale height $H_\rho $ in brackets (1/$H_\rho $ = -|dln $\rho $/dr| in 108 cm) at positions given by the closest black dots, respectively. Right panel: chemical composition of the initial model. The dotted vertical lines mark the edges of the convection zone.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg41.eps}
\end{figure} Figure 2:

Temporal evolution of the total nuclear energy production rate E of models hefl.3d (solid) and hefl.2d.a (dotted), respectively. $L_{\odot }$ is the solar luminosity.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg42.eps}
\end{figure} Figure 3:

Temporal evolution of the total kinetic energy $E_{\rm K}$ of models hefl.3d (solid) and hefl.2d.a (dotted), respectively.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8.5cm,clip]{1414fg43.eps}\par\includegraphics[width=8.5cm,clip]{1414fg44.eps}
\end{figure} Figure 4:

Temporal evolution of a log10 of the angular averaged kinetic energy density (in $\mbox{~erg }\mbox{~g}^{-1}$) of models hefl.3d ( upper) and hefl.2d.a ( lower), respectively.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=5.6cm,clip]{1414f45a.eps}\hspace*{0.5c...
...ps}\hspace*{1.5cm}
\includegraphics[width=7.9cm,clip]{1414f45e.eps}
\end{figure} Figure 5:

Upper panels: snapshots taken at ${\sim }4754$ s showing contour plots of the absolute value of the velocity (in units of $10^{6} \mbox{~cm~s${}^{-1}$ }$) for the 2D model hefl.2d.a ( left), and for the 3D model hefl.3d in a meridional plane of azimuth angle $\phi = 50\mbox{$^\circ$ }$( middle) and $\phi = 70\mbox{$^\circ$ }$ ( right), respectively. Lower panels: normalized power spectra of angular fluctuations in the absolute velocity as a function of radius for the 2D model helf.2d.a ( left) averaged over time from 2100 s to 9500 s, and for the 3D model hefl.3d ( right) averaged over time from 2250 s to 6000 s, and azimuthal angle. The dashed vertical lines mark the edges of the convective zone of the initial model according to the Schwarzschild criterion.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=5.9cm,clip]{1414f46a.eps}\hspace*{0.2c...
...s}\hspace*{0.2cm}
\includegraphics[width=5.9cm,clip]{1414f46i.eps}
\end{figure} Figure 6:

Cuts through the 3D model hefl.3d at t = 4815 s showing the angular variation of temperature (in units of 108 K; left panels), 12C mass fraction (in units of 10-3; middle panels), and radial velocity (in units of $10^{5} \mbox{~cm~s${}^{-1}$ }$; right panels), respectively, at three different radii: r1 = 4.8 $\times $ 108 cm (temperature maximum; top), r2 = 6.5 $\times $ 108 cm (center of the convection zone; middle), and r3 = 9.3 $\times $ 108 cm ( top of convection zone, bottom). The black lines in the right panels mark the boundaries between positive and negative radial velocities i.e., upflows and downflows.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg47.eps}
\end{figure} Figure 7:

Auto-correlation function A1(r0;r) measuring the radial extent of flow patterns (see Eq. (1)) at different radii r0 (4.8 $\times $ 108 cm, 5.9 $\times $ 108 cm, 7 $\times $ 108 cm, 8.1 $\times $ 108 cm, and 9.25 $\times $ 108 cm) and $t \sim 4000$ s measuring the radial extent of flow patterns for the 2D model hefl.2d.a (dotted lines) and the 3D hefl.3d (solid lines), respectively.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg48.eps}
\end{figure} Figure 8:

Auto-correlation function A2(r0;t-t0) measuring the lifetime of flow patterns (see Eq. (2)) at two different epochs for the 3D model hefl.3d (solid: $t_{0} \sim 2260$ s, dashed: $t_{0} \sim 2900$ s), and for the 2D model hefl.2d.a (crosses, diamonds), respectively. The radius r0 is 7.6 $\times $ 108 cm.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=7.4cm,clip]{1414f49a.eps}\vspace*{5mm}...
....eps}\vspace*{5mm}
\includegraphics[width=7.4cm,clip]{1414f49c.eps}
\end{figure} Figure 9:

Carbon mass fraction X(12C) ( top), entropy S ( middle), and entropy gradient $\nabla S$ ( bottom) as a function of radius near the outer edge of the convection zone of model hefl.3d at three different epochs: t1 = 2000 s (dashed), t2 = 4000 s (dash-dotted), and t3 = 6000 (solid). In addition, the initial X(12C) profile is shown in the top panel (dotted).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414f410.eps}
\end{figure} Figure 10:

Square of the Brunt-Väisälä buoyancy frequency as a function of radius for the 3D model hefl.3d at t = 4638 s.

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\includegraphics[width=8.5cm,clip]{1414iso1.eps}\includ...
...ps}\hspace{2.5cm}
\includegraphics[width=6cm,clip]{1414iso4.eps} }\end{figure} Figure 11:

Different views of isosurfaces of the velocity field for the 3D model hefl.3d at t = 3000 s. The blue isosurface corresponds to radial downflows with vr = -6 $\times $ $10^{5} \mbox{~cm~s${}^{-1}$ }$, and the yellow and brown isosurfaces show radial upflows with vr = 6 $\times $ $10^{5} \mbox{~cm~s${}^{-1}$ }$, and 1 $\times $ $10^{4}\mbox{~cm~s${}^{-1}$ }$ (gravity waves), respectively. The edge sizes of the box are 1.2 $\times $ 109 cm and 2.4 $\times $ 109 cm, respectively. The yellow-greenish sphere in the bottom right panel marks the top of the convection zone according to the Schwarzschild criterion.

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\includegraphics[width=8.8cm]{1414flx1.eps}\includegrap...
...=8.8cm]{1414flx3.eps}\includegraphics[width=8.8cm]{1414flx4.eps} }\end{figure} Figure 12:

Various energy fluxes and source terms as a function of radius averaged (from 2000 s to 6000 s) over about four convective turnover timescales. Panel  a) shows the convective flux $F_{\rm C}$ of the 2D model hefl.2d.a (solid-red) and the 3D model hefl.3d (dashed-black) together with the kinetic energy flux $F_{\rm K}$ in the 2D (dash-dotted-red) and 3D (dash-dot-dotted-black) model, respectively. The dotted vertical lines mark the edges of the convection zone in the initial model according to the Schwarzschild criterion. Panel  b) gives the source terms due to the work done by buoyancy forces PA (dash-dot-dotted black) and due to volume changes PP (dashed black) in the 3D model hefl.3d, and in the 2D model hefl.2d.a (dash-dotted-red, solid red) respectively. The vertical lines enclose the nuclear burning zone (T > 108 K). Panels  c) and  d) show an enlarged view of the energy fluxes and source terms displayed in panels  a) and  b) near the outer edge of the convection zone.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414ve2d.eps}\hspace*{1cm}
\includegraphics[width=8cm,clip]{1414ve3d.eps}
\end{figure} Figure 13:

Radial distributions of the time (from 2000 s to 4000 s) and angle-averaged velocity components (vr, $v_\theta $, $v_\phi $) and velocity modulus ( $v_{\rm abs}$) for the 2D model hefl.2d.a ( left), and the 3D model hefl.3d ( right), respectively. The panels also show the velocity predicted by the mixing-length theory ( $v_{\rm mlt}$).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414f414.eps}
\end{figure} Figure 14:

Radial distributions of the time (from 2800 s to 4800 s) and angle-averaged density fluctuations in the 3D model hefl.3d (solid red line). The panels also show the angular variation of the respective quantity at a given radius (gray shaded region). The dotted vertical lines mark the edges of the convection zone as determined by the Schwarzschild criterion.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414f415.eps}
\end{figure} Figure 15:

Fractional volume occupied by upflow streams $A_{{\rm up}}$ averaged over $\sim $2000 s as a function of radius for the 3D model hefl.3d (solid), the 2D low-resolution model hefl.2d.a (dashed-dotted), and the 2D high-resolution model hefl.2d.b (dotted), respectively.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=7.8cm,clip]{1414adg1.eps}\vspace*{0.5cm}
\includegraphics[width=7.8cm,clip]{1414adg2.eps}
\end{figure} Figure 16:

Upper panel: radial distributions of the adiabatic temperature gradient  $\nabla _{\rm ad}$ (dotted) and of the temperature gradients of model hefl.3d (solid), respectively. The latter distribution is a linear fit to the gradients averaged over angle and over the first 200 s of the evolution of the model. The gray shaded region marks the convection zone CVZ. Lower panel: same as above, but showing the radial distributions in the evolved convection zone averaged over roughly 3000 s of evolutionary time.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=5.7cm,clip]{1414vl1.eps}\hspace*{5mm}
...
...s}\hspace*{5mm}
\includegraphics[width=5.7cm,clip]{1414vl3.eps}
\par\end{figure} Figure 17:

Snapshots of the spatial distribution of the velocity modulus |v| (in units of $10^{6} \mbox{~cm~s${}^{-1}$ }$) for the 2D model hefl.2d.b at 24 000 s ( left), 60 000 s ( middle), and 120 000 s ( right), respectively.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm,clip]{1414fg62.eps}
\end{figure} Figure 18:

Temporal evolution of the horizontally averaged temperature maximum  $\langle T \rangle_{\rm max}$ (solid), and of the global temperature maximum  $T_{\rm max}$ (solid thin) in the long-term 2D model hefl.2d.b. The dashed line correspond to the temporal evolution of the maximum temperature in the stellar evolutionary calculation.

Open with DEXTER
In the text

  \begin{figure}
\par\mbox{\includegraphics[width=8.5cm]{1414f63a.eps}\includegrap...
...8.5cm]{1414f63c.eps}\includegraphics[width=8.5cm]{1414f63d.eps} }\end{figure} Figure 19:

Upper panels: angular averaged 12C mass fraction as a function of radius near the inner ( left) and outer edge ( right) of the convection zone in the long-term 2D model hefl.2d.b at $t = 60~000\mbox{~s}$ (dashed) and $t = 120~000\mbox{~s}$ (dash-dotted), respectively, with temperature stratification (thick) at $t = 120~000\mbox{~s}$. The vertical dotted lines mark the boundaries of the convection zone at t = 0 s. Lower panels: temporal evolution of the position of the inner ( left) and outer ( right) edge of the convection zone in model hefl.2d.b, respectively.

Open with DEXTER
In the text


Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.