Letter to the Editor
A test of timedependent theories of stellar convection^{⋆}
^{1}
MaxPlanckInstitut für Sonnensystemforschung, MaxPlanckStrasse 2, 37191 KatlenburgLindau, Germany
email: gastine@mps.mpg.de
^{2}
Institut de Recherche en Astrophysique et Planétologie, CNRS/Université de Toulouse, 14 Av. Edouard Belin, 31400 Toulouse, France
Received: 22 February 2011
Accepted: 16 March 2011
Context. In Cepheids close to the red edge of the classical instability strip, a coupling occurs between the acoustic oscillations and the convective motions close to the surface. The best topical models that account for this coupling rely on 1D timedependent convection (TDC) formulations. However, their intrinsic weakness comes from the large number of unconstrained free parameters entering into the description of turbulent convection.
Aims. We compare two widely used TDC models with the first 2D nonlinear direct numerical simulations (DNS) of the convectionpulsation coupling in which the acoustic oscillations are selfsustained by the κmechanism.
Methods. The free parameters appearing in the Stellingwerf and Kuhfuß TDC recipes are constrained using a χ^{2}test with the timedependent convective flux that evolves in nonlinear simulations of highlycompressible convection with the κmechanism.
Results. This work emphasises some inherent limits to TDC models, that is, the temporal variability and nonuniversality of their free parameters. More importantly, within these limits, Stellingwerf’s formalism is found to give better spatial and temporal agreements with the nonlinear simulation than Kuhfuß’s one. It may therefore be preferred in 1D TDC hydrocodes or stellar evolution codes.
Key words: hydrodynamics / convection / stars: variables: Cepheids / stars: oscillations / methods: numerical
Appendix A is available in electronic form at http://www.aanda.org
© ESO, 2011
1. Introduction
The first theoretical calculations of the Cepheids instability strip done in the 60’s assumed that convection was steady with respect to oscillations. Unfortunately, this “frozenin convection” approximation led to a cooler red edge than the observed one, because the strong coupling between convection and pulsations occurring in cool Cepheids was ignored (e.g. Baker & Kippenhahn 1965). Then, following the pioneering works of Unno (1967) and Gough (1977), several timedependent convection (TDC) models have been developed to investigate the influence of the convection on the pulsational stability (e.g. Stellingwerf 1982; Kuhfuß 1986; Gehmeyr & Winkler 1992a). The last uptodate TDC models actually succeed in predicting both a red edge close to the observed one and realistic luminosity curves (e.g. Yecko et al. 1998; Bono et al. 1999; Feuchtinger 1999; Kolláth et al. 2002).
However, all of these TDC models suffer from a common weakness due to the numerous free parameters, usually known as α coefficients, which describe the turbulent convection (e.g., the eight dimensionless parameters in Smolec & Moskalik 2010). These parameters are either obtained from a fit to the observations or hardly constrained by theory when taking the asymptotic limit of stationary convection (Vitense 1953; BöhmVitense 1958). A parametric study carried out by Yecko et al. (1998) emphasises the intrinsic degeneracy of TDC models, because similar instability strips have been obtained with different sets of parameters.
Direct numerical simulations (hereafter DNS) are able to constrain these TDC models because they fully account for the nonlinearities involved in the convectionpulsation coupling. These nonlinear simulations are challenging as they require both largeamplitude oscillations and convective motions. In our last 2D simulations, we improved the way acoustic waves are generated by reproducing the selfconsistent excitation operating in variable stars, known as the κmechanism (Gastine & Dintrans 2011, hereafter GD2011). The resulting coupling of acoustic modes with convection is therefore more consistent because the mode amplitude is not imposed artificially. We showed in GD2011 that the convective plumes may either quench the radial oscillations or coexist with the acoustic modes, depending mainly on the density contrast of the equilibrium model.
The purpose of this letter is to compare the fully nonlinear results with two main TDC models: (i) the first one refers to an initial formulation of Stellingwerf (1982) that was used and improved by Bono & Stellingwerf (1994) and Bono et al. (1999); (ii) the other one was developed by Kuhfuß (1986) and Gehmeyr & Winkler (1992a) and implemented both in the Vienna hydrocode (e.g. Wuchterl & Feuchtinger 1998; Feuchtinger 1999) and the FloridaBudapest one (e.g. Yecko et al. 1998; Kolláth et al. 2002). In these models, a single equation for the turbulent kinetic energy ℰ_{t} is added to the classical meanfield equations and the main secondorder correlations, such as the convective flux, are expressed as a function of ℰ_{t} alone (see e.g. Baker 1987; Gehmeyr & Winkler 1992b; Buchler 2000, 2009).
We investigate here a particular simulation where the oscillations strongly modulate the convective flux in more detail. The nonlinear results are first compared at each snapshot with the TDC recipes by computing χ^{2}statistics of the relevant α coefficients. Second, the mean values of α are used to compare the optimal TDC fluxes with the DNS ones. The formulation of Stellingwerf is found to be closer to the nonlinear result than the Kuhfuß one: (i) the temporal variation of the α coefficient is weaker; (ii) the mean convective flux is closer to the DNS one, especially in the description of the overshooting area.
2. The hydrodynamic model
We consider a local 2D box of size L_{x} × L_{z}, filled by a perfect monatomic gas and centred on both sides of an ionisation region. Its associated opacity bump is shaped by a hollow in the temperaturedependent radiative conductivity profile K(T). Such a configuration can lead to unstable acoustic modes due to the driving term ∇·K(T)∇T in the energy equation; i.e., this is the κmechanism (Gastine & Dintrans 2008a,b). Furthermore, this hollow in K(T) is deep enough that the equilibrium temperature gradient is locally superadiabatic and convective motions develop here according to Schwarzschild’s criterion.
The governing hydrodynamic equations are written in nondimensional form by choosing L_{z} as the length scale and as the velocity one, hence the time scale (with c_{p} the specific heat and T_{top} the surface temperature). The resulting nonlinear set of equations is advanced in time with the highorder, finitedifference pencil code^{1}, which is fully explicit except for the radiative diffusion term that is solved implicitly thanks to a parallel alternate direction implicit (ADI) solver. With the chosen units, the simulation box spans about 10% of the star radius on both sides of the ionisation region, while the timestep is about 1 minute, such that a simulation typically spans over 4500 days (see GD2011).
3. Results
The time evolution of the convective flux obtained in fully nonlinear 2D simulations is compared with the following TDC expressions developed by Stellingwerf (1982) and Kuhfuß (1986):
(1)where ∇ = dlnT/dlnp, ∇_{ad} = 1 − c_{v}/c_{p} and
(2)with p and ρ the pressure and density, respectively, and the brackets denote a horizontal average. Two free, dimensionless parameters α_{St} and α_{Ku} are introduced that the simulations allow to constrain. We first focus on a 2D DNS that is similar to the G8 simulation in GD2011, which is a simulation in which the total kinetic energy is almost entirely contained in acoustic modes excited by the κmechanism (80%), with the rest in convective plumes (20%). The convective motions do not affect the growth and the nonlinear saturation of the unstable radial acoustic mode much. We can therefore expect that the heat transport is strongly modulated by wave motions, which is an ideal frame for testing the accuracy of timedependent convection models.
3.1. Temporal modulation of the convective flux in the DNS
Fig. 1 a) Mean vertical profiles of radiative ℱ_{rad} (solid blue line), turbulent enthalpy ℱ_{conv} (dashed green line), kinetic ℱ_{kin} (dotted black line), modes ℱ_{mod} (dotdashed magenta line), and total ℱ_{tot} (dotted red line) fluxes, normalised to the bottom flux F_{bot}. b) Temporal power spectrum for the convective flux alone. 

Open with DEXTER 
Fig. 2 Evolution of the convective flux in a (t, z) plane: a) in the DNS; b) with Stellingwerf’s formalism ℱ_{St}; c) with Kuhfuß’s formalism ℱ_{Ku} (for α_{St} = α_{Ku} = 1). The vertical extent is centred on both sides of the convection zone to emphasise its oscillations. These snapshots were computed after 1800 periods of oscillations, i.e. well after the nonlinear saturation was achieved. 

Open with DEXTER 
The imposed bottom flux F_{bot} is transported mainly through the computational domain by the radiative ℱ_{rad}, enthalpy, and kinetic ℱ_{kin} fluxes. Following GD2011, the enthalpy flux is divided into the classical convective flux ℱ_{conv} and the contribution coming from the (unstable) fundamental mode (hereafter ℱ_{mod}) with
(3)where the primed quantities denote the fluctuations about the horizontal average, and θ is the temperature eigenfunction of the fundamental mode. The resulting timeaveraged and normalised fluxes are given in Fig. 1a.
The bulk of the total flux is transported by the radiative flux, except in the convective zone where ℱ_{conv}/F_{bot} ≃ 20%, while the kinetic flux is negligible (ℱ_{kin}/F_{bot} ≤ 1%). The value of ℱ_{mod} is hardly as large as ℱ_{conv} in the convective zone. This quantity is a good signature of the amplitude of the acoustic modes – the higher ℱ_{mod}, the larger the radial oscillations (GD2011) – therefore confirming the efficiency of the κmechanism in this simulation.
The signature of the temporal modulation of the convective flux is extracted from its Fourier spectrum in time; that is, we first compute , with ω the frequency and then integrate over the vertical direction z to get the power spectrum (Fig. 1b). Several discrete peaks appear about given frequencies, of which the physical nature is emphasised after superimposing the linear acoustic eigenfrequencies (the vertical dashed blue lines). The match with the fundamental mode frequency ω_{00} = 3.85 is perfect, while the weakamplitude secondary peaks instead correspond to harmonics of ω_{00} (i.e. 2ω_{00},3ω_{00},... , shaped by downwarddirected vertical grey arrows in Fig. 1b) than to the acoustic overtones ω_{01}, ω_{02},... It means that the amplitude of the fundamental acoustic mode is large enough to generate several harmonics through a nonlinear cascade, and this again indicates both the robustness of the κdriving and the relevance of this kind of convectionpulsation simulation to check the TDC recipes.
3.2. DNS vs. TDC models: convective patterns
We first compare the time evolution of the horizontally averaged convective flux obtained in the DNS (Eq. (3)) with its TDC counterparts ℱ_{St} and ℱ_{Ku} (Eq. (1)). As we are interested here in the qualitative agreement between the convective patterns in the plane (t, z), we simply assume α_{St} = α_{Ku} = 1.
The three resulting patterns are displayed in Fig. 2 for a time interval spanning about four periods of the fundamental acoustic mode. The black areas denote positive values for the convective flux and therefore delimit the convective zone. An oscillatory behaviour is clearly visible in each panel, with a period that looks similar to the one of the fundamental acoustic mode; that is, almost 4 oscillation cycles are depicted. This is consistent with the large peak shown around ω_{00} in the power spectrum in Fig. 1b.
The two TDC patterns also display good agreement with the nonlinear simulation regarding the overshooting phenomenon. We indeed recover the same dark red structures below the convective zone, and they correspond to the downdrafts entering the radiative zone (where this penetration is associated with a negative convective flux). Nevertheless, we note that the overshooting seems to be more vigorous in Kuhfuß’s model than in Stellingwerf’s, because these dark red structures fill a larger surface in Fig. 2c in the bottom radiative zone than in Fig. 2b.
3.3. DNS vs. TDC models: statistics of coefficients α
Fig. 3 Time evolution around the mean value of coefficients α_{St} (upper panel) and α_{Ku} (bottom panel). Horizontal grey spans mark the limits of the associated relative standard deviation. 

Open with DEXTER 
A onetoone comparison between the convective flux in the DNS and its TDC predictions requires finding the optimal values of α_{St} and α_{Ku}. This is done by performing several χ^{2}tests at different snapshots in the simulation to track the variations in coefficients α. The resulting fluctuations of α_{St} and α_{Ku} around their mean values are shown in Fig. 3.
The dispersion of the Stellingwerf coefficient α_{St} (upper panel) is weaker than the Kuhfuß one (lower panel), since its values are almost within a 5% range around the mean . In contrast, several outliers with values greater than 10% are found in the evolution of α_{Ku} which then appears more chaotic around the mean value . As a consequence, the relative standard deviation (depicted in gray in Fig. 3) is weaker in Stellingwerf’s case than in Kuhfuß’s.
3.4. DNS vs. TDC models: convective fluxes with optimal α’s
We recall that TDC models assume that the coefficients α entering Eq. (1) are constant. The optimal value of these coefficients has been deduced by adjusting the TDC recipes with the instantaneous convective flux throughout the nonlinear simulation. The final test, given in Fig. 4, then consists in comparing the mean nonlinear convective flux taken over the entire simulation and its best TDC approximations built from these optimal α values.
This figure emphasises that Stellingwerf’s formulation gives a better agreement than Kuhfuß’s one. Indeed, the Kuhfuß model overestimates the overshooting as the (negative) convective flux remains nonnegligible until the bottom of the radiative zone. In contrast, the Stellingwerf profile accounts for the local penetration of convective plumes better and shows the same exponentiallike decay in the negative convective flux when sinking in the radiative zone (e.g. Dintrans 2009). However, the two models are fairly similar in the bulk of the convective zone where convection is fully developed. One also notes that they both predict a negative flux at the top of the convective zone, which is an upper overshooting of convection motions near the surface that is not observed in the nonlinear simulation.
4. Conclusion
Fig. 4 Mean convective flux in the DNS (solid black line), compared with the best TDC predictions based on the optimal α values that came out of the χ^{2}test in Sect. 3.3 (Stellingwerf’s model is a dashed blue line and Kuhfuß’s a dotted green line). 

Open with DEXTER 
The main weakness of all theories of timedependent convection (TDC) lies in the large number of free parameters. This is particularly awkward when convection is strongly coupled with pulsations, such as near the red border of the Cepheid instability strip where similar results are obtained with different sets of parameters (Yecko et al. 1998). Additional constraints must be found to reduce the intrinsic degeneracy of models. But the modelling of the interplay between the turbulent convection and oscillations is really a difficult task owing to the strong nonlinearities involved in this coupling.
One solution to tackle this problem and to bring new constraints consists in performing fully nonlinear simulations, and this is the path we chose in this work following our first study in GD2011. Two widely used TDC theories, namely the Stellingwerf (1982) and Kuhfuß (1986) ones, are compared with results coming from 2D nonlinear simulations of compressible convection in which strong acoustic oscillations are selfsustained by the κmechanism. The heat transport is then modulated by the fundamental acoustic mode such that this kind of simulation is relevant for investigating the convectionpulsation coupling (Figs. 1, 2).
Focusing on the two TDC formulae for the convective flux, we computed the evolution of free coefficients “α” from a χ^{2}test applied to the fully nonlinear results (Fig. 3). A large temporal variability is found in both cases, and it weakens the robustness of the TDC assumption α = const. Moreover, the mean values are not universal. By applying the same method to other simulations performed in GD2011 (Table 1), we indeed do not recover the same with a relative standard deviation of about 12% (Stellingwerf) and 18% (Kuhfuß).
Values of the optimal Stellingwerf and Kuhfuß α coefficients pulled out of the nonlinear simulations in GD2011.
Within these limits, Stellingwerf’s formulation is found to agree with the nonlinear simulations better than does Kuhfuß’s: (i) the final mean convective flux ℱ_{St} is closer to its DNS counterpart with a much better estimation of the bottom overshooting; (ii) the temporal dispersion of the α_{St} coefficient is weaker; that is, the stability is enhanced (Fig. 4). This result means that the timedependent convective flux scales better with a law than , and the Stellingwerf formulation may probably be preferred in the 1D hydrocode used in, e.g., the topical Cepheids models. However, this study emphasises that both formalisms lead to an artificial overshooting at the top of the convection zone. We also noted that this TDC test involves the exact value of the turbulent kinetic energy ℰ_{t} provided by the nonlinear simulation, and not by the dedicated 1D TDC equation, which is inherently less accurate. As a consequence, the obtained profiles are certainly the best we can expect from these TDC recipes.
In this work, the temporal modulation of the convective flux was ensured by the acoustic modes excited by κmechanism. An interesting prospect could be the case of a modulation based on the internal gravity waves excited by convection itself. Indeed, convection can excite gravity waves in variable stars, either by convective elements penetrating stably stratified regions, such as in solartype stars (e.g. Dintrans et al. 2005); or through the socalled “convective blocking” mechanism, as in white dwarfs or Gamma Doradus stars (e.g. Pesnell 1987). In both cases, the convective flux is ipso facto modulated by gravity waves, and it may be interesting in that respect to also check the accuracy of timedependent convection models.
Acknowledgments
This work was granted access to the HPC resources of CALMIP under the allocation 2010P1021 (http://www.calmip.cict.fr).
References
 Baker, N., & Kippenhahn, R. 1965, ApJ, 142, 868 [NASA ADS] [CrossRef] (In the text)
 Baker, N. H. 1987, in Physical Processes in Comets, Stars and Active Galaxies, ed. W. Hillebrandt, E. MeyerHofmeister, & H.C. Thomas (Berlin & New York: SpringerVerlag), 105 (In the text)
 BöhmVitense, E. 1958, Z. Astrophys., 46, 108 [NASA ADS] (In the text)
 Bono, G., & Stellingwerf, R. F. 1994, ApJS, 93, 233 [NASA ADS] [CrossRef] (In the text)
 Bono, G., Marconi, M., & Stellingwerf, R. F. 1999, ApJS, 122, 167 [NASA ADS] [CrossRef] (In the text)
 Brandenburg, A., & Dobler, W. 2002, CoPhC, 147, 471 [NASA ADS] [CrossRef] (In the text)
 Buchler, J. R. 2000, ASPCS, 203, 343 [NASA ADS] (In the text)
 Buchler, J. R. 2009, in AIPCS, ed. J. A. Guzik, & P. A. Bradley (Melville: American Institute of Physics), 1170, 51 (In the text)
 Dintrans, B. 2009, CoAst, 158, 45 [NASA ADS] (In the text)
 Dintrans, B., Brandenburg, A., Nordlund, Å., & Stein, R. F. 2005, A&A, 438, 365 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Feuchtinger, M. U. 1999, A&A, 136, 217 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gastine, T., & Dintrans, B. 2008a, A&A, 484, 29 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gastine, T., & Dintrans, B. 2008b, A&A, 490, 743 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gastine, T., & Dintrans, B. 2011, A&A, 528, A6 [GD2011] [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gehmeyr, M., & Winkler, K.H. A. 1992a, A&A, 253, 92 [NASA ADS] (In the text)
 Gehmeyr, M., & Winkler, K.H. A. 1992b, A&A, 253, 101 [NASA ADS] (In the text)
 Gough, D. O. 1977, ApJ, 214, 196 [NASA ADS] [CrossRef] (In the text)
 Kolláth, Z., Buchler, J. R., Szabó, R., & Csubry, Z. 2002, A&A, 385, 932 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Kuhfuß, R. 1986, A&A, 160, 116 [NASA ADS] (In the text)
 Pesnell, W. D. 1987, ApJ, 314, 598 [NASA ADS] [CrossRef] (In the text)
 Smolec, R., & Moskalik, P. 2010, A&A, 524, A40 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Stellingwerf, R. F. 1982, ApJ, 262, 330 [NASA ADS] [CrossRef] (In the text)
 Unno, W. 1967, PASJ, 19, 140 [NASA ADS] (In the text)
 Vitense, E. 1953, Z. Astrophys., 32, 135 [NASA ADS] (In the text)
 Wuchterl, G., & Feuchtinger, M. U. 1998, A&A, 340, 419 [NASA ADS] (In the text)
 Yecko, P. A., Kolláth, Z., & Buchler, J. R. 1998, A&A, 336, 553 [NASA ADS] (In the text)
Online material
Appendix A: Setup of the simulations
A.1. The equilibrium model
As in GD2011, our system represents a zoom around an ionisation region. Since we are computing local simulations, the vertical gravity g = −ge_{z} and the kinematic viscosity ν are assumed to be constant. Following our purely radiative model of the κmechanism (Gastine & Dintrans 2008a,b), the ionisation region is represented by a temperaturedependent radiative conductivity profile that mimics an opacity bump:
with
(A.2)where T_{bump} is the position of the hollow in temperature and σ, e, and denote its slope, width, and relative amplitude, respectively. We assume both radiative and hydrostatic equilibria; that is,
where F_{bot} is the imposed bottom flux. Following GD2011, we chose L_{z} as the length scale, i.e. [x] = L_{z}, top density ρ_{top} and top temperature T_{top} as density and temperature scales, respectively. The velocity scale is then , while time is given in units of .
Table A.1 then summarises the parameters of the numerical simulations presented in this study in these dimensionless units. The penultimate column of this table contains the value of the frequency ω_{00} of the fundamental unstable radial mode excited by the κmechanism, which lies between 3 and 4 for every DNS. The last column gives the value of the Rayleigh number, which quantifies the strength of the convective motions. It is given by
where L_{conv} is the width of the convective zone, χ = K_{0}/ρ_{0}c_{p} the radiative diffusivity, and s the entropy.
Dimensionless parameters of the numerical simulations.
A.2. The nonlinear equations
With the parallel version of the alternate direction implicit (ADI) solver for the radiative diffusion implemented in the pencil code (see GD2011), we advance the following hydrodynamic equations in time:
where ρ, u, and T denote density, velocity, and temperature, respectively, while K(T) is given by Eq. (A.1). The operator D/Dt = ∂/∂t + u·∇ is the usual total derivative, while S is the (traceless) rateofstrain tensor given by
Finally, we impose the condition that all fields are periodic in the horizontal direction, while stressfree boundary conditions (i.e., u_{z} = 0 and du_{x}/dz = 0) are assumed for the velocity in the vertical one. Concerning the temperature, a perfect conductor at the bottom (i.e., flux imposed) and a perfect insulator at the top (i.e., temperature imposed) are applied.
To ensure that both the nonlinear saturation and thermal relaxation are achieved, the simulations were computed over very long times, typically t ≳ 3000. As the eigenfrequency of the unstable acoustic mode ω_{00} ∈ [3 − 4] (see Table A.1), this corresponds approximately to 1800 periods of oscillation.
All Tables
Values of the optimal Stellingwerf and Kuhfuß α coefficients pulled out of the nonlinear simulations in GD2011.
All Figures
Fig. 1 a) Mean vertical profiles of radiative ℱ_{rad} (solid blue line), turbulent enthalpy ℱ_{conv} (dashed green line), kinetic ℱ_{kin} (dotted black line), modes ℱ_{mod} (dotdashed magenta line), and total ℱ_{tot} (dotted red line) fluxes, normalised to the bottom flux F_{bot}. b) Temporal power spectrum for the convective flux alone. 

Open with DEXTER  
In the text 
Fig. 2 Evolution of the convective flux in a (t, z) plane: a) in the DNS; b) with Stellingwerf’s formalism ℱ_{St}; c) with Kuhfuß’s formalism ℱ_{Ku} (for α_{St} = α_{Ku} = 1). The vertical extent is centred on both sides of the convection zone to emphasise its oscillations. These snapshots were computed after 1800 periods of oscillations, i.e. well after the nonlinear saturation was achieved. 

Open with DEXTER  
In the text 
Fig. 3 Time evolution around the mean value of coefficients α_{St} (upper panel) and α_{Ku} (bottom panel). Horizontal grey spans mark the limits of the associated relative standard deviation. 

Open with DEXTER  
In the text 
Fig. 4 Mean convective flux in the DNS (solid black line), compared with the best TDC predictions based on the optimal α values that came out of the χ^{2}test in Sect. 3.3 (Stellingwerf’s model is a dashed blue line and Kuhfuß’s a dotted green line). 

Open with DEXTER  
In the text 