Table Of Content

AstrophysicsandSpaceScience DOI10.1007/s•••••-•••-••••-• Numerical simulations of the κ-mechanism with convection T. Gastine1 • B. Dintrans1 0 1 0 2 n a J (cid:13)c Springer-Verlag•••• 5 Abstract A strong coupling between convection and These driving zones should nevertheless be located at ] pulsations is known to play a major role in the dis- a precise radius inside the star, neither too close to R appearance of unstable modes close to the red edge of the surface nor to deep, in order to balance the damp- S the classical Cepheid instability strip. As mean-field ing that mainly occurs at the surface. It defines the . h models of time-dependent convection rely on weakly- so-called transition region that shapes the limit be- p constrainedparameters(e.g.Baker1987),wetacklethis tweenthequasi-adiabaticinteriorandthestronglynon- - o problem by the means of 2-D Direct Numerical Simu- adiabatic surface. For classical Cepheids that pul- tr lations (DNS) of κ-mechanism with convection. sate in the fundamental acoustic mode, this transi- s Using a linear stability analysis, we first determine tion region is located at a temperature T ≃ 4 × a [ the physical conditions favourable to the κ-mechanism 104 K corresponding to the second helium ionisation to occur inside a purely-radiative layer. Both the in- (Baker and Kippenhahn 1965). However, the bump lo- 2 stabilitystripsandthenonlinearsaturationofunstable v cationisnotsolelyresponsiblefortheacousticinstabil- modes are then confirmed by the corresponding DNS. 7 ity. A careful treatment of the κ-mechanism would in- 8 We next present the new simulations with convection, volve dynamical couplings with convection, metallicity 4 where a convective zone and the driving region over- effects and both realistic equations of state and opac- 0 lap. The coupling between the convective motions and . itytables(Bono, Marconi, and Stellingwerf1999). The 1 acoustic modes is then addressed by using projections purpose of our model is to simplify the hydrodynamic 0 onto an acoustic subspace. 0 approachwhileretainingtheleading-order phenomenon 1 –e.g. thelocationoftheopacitybumporitsamplitude– Keywords Hydrodynamics - Instabilities - Stars: os- v: cillations - Convection - Methods: numerical suchthatfeasibledirectnumericalsimulationsoftheκ- i mechanism with convection can be achieved. X r a 1 Introduction 2 The first step: radiative models in 1-D Eddington (1917) discovered an excitation mechanism of stellar oscillations that is related to the opacity be- Wefirstfocusonradialmodespropagatinginapurely- haviour in ionisation zones: the κ-mechanism, where radiative and partially-ionised layer and then restrict κ denotes the opacity. This mechanism occurs when our study to the 1-D case. Our model consists of a the opacity varies during compression phases so as local zoom about an ionisation region and is composed to block the emerging radiative flux (Zhevakin 1963; byamonatomicandperfectgas(withγ =c /c =5/3) p v Cox and Whitney1958). Ionisationregionscorrespond underbothaconstantgravity~gandkinematicviscosity to a strong increase in opacity that leads to the opac- ν. All quantities in our model are dimensionless, e.g. ity bumps responsibleforthelocalexcitationofmodes. temperatureisgiveninunitofthesurfacetemperature ofthestar,lengthinunitofafractionofthestarradius T.Gastine and velocity in unit of the sound speed at the surface. B.Dintrans In order to easily investigate the influence of the 1LATT,UniversitédeToulouse, CNRS,14 avenueEdouard Be- opacity bump on the stability, we adopt the following lin,F-31400Toulouse,France 2 Fig. 1 Influence of the hollow parameters on the conduc- −2 tivity profile for Kmax =10 and Tbump/Tsurf =3.5. conductivity hollow to mimic this bump (as κ ∝ 1/K, see Gastine and Dintrans 2008a, hereafter Paper I): Fig. 2 a) Vertical profile of the Ψ coefficient (Eq. 3). Thesuperimposedgreenzonerepresentsthelocationwhere Ψ=1±0.5. b) The real part of the work integral plotted −π/2+arctan(σT+T−) for the three different equilibrium models. The dot-dashed K (T )=K 1+A , (1) 0 0 max(cid:20) π/2+arctan(σe2) (cid:21) black line is for the constant radiative conductivity case. c) Corresponding radiative conductivity profiles. d) Corre- withA=(Kmax−Kmin)/Kmax andT± =T0−Tbump± sponding equilibrium fields dKT/dz (Eq. 2). e. HereT istheequilibriumtemperatureprofile,T 0 bump is the hollow central temperature, and σ, e and A de- luminosity. Ψ is the ratio between the thermal energy note its slope, width and relative amplitude, respec- embeddedbetweenthegivenradiusandthesurfaceand tively. Examplesofcommonvaluesoftheseparameters the energy radiated during an oscillation period. are provided in Fig. 1. 2.2 The linear stability analysis 2.1 Conditions for the instability Inordertocheckthesetwoconditions,weperformalin- Using the well-known work integral formalism (e.g. earstabilityanalysisbyusingtheLinearSolverBuilder Unno et al. 1989), it is easy to show that the follow- code(Valdettaro et al.2007). Bystartingwiththecon- ing condition is necessary to get unstable modes by κ- ductivity profiles shown in Fig. 1, we first compute the mechanism in variable stars: equilibrium setup that is solution of both the hydro- static and radiative equilibria given by dK ∂lnK T 0 <0 where K ≡ . (2) T dz (cid:18)∂lnT (cid:19) 0 ρ0 ∇~p =ρ ~g and div[K (T )∇~T ]=0, (4) 0 0 0 0 0 However, this condition is not sufficient as the ionisation region should also be located at a given location where p0 and ρ0 denote the equilibrium pressure and in the star, neither too deep nor too close to its sur- density. Then, we solve the linear oscillation equations face. This area is the so-called transition region that for the perturbations by seeking normal modes of the separatesthequasi-adiabaticinteriorfromthestrongly form exp(λt) with λ = τ +iω, where ω is the mode non-adiabatic surface and is defined by (e.g. Cox 1980) frequency and τ its growth or damping rate (unstable modes corresponding to τ > 0). By varying the hollow parameters, we then determine the physical con- hc T i∆m Ψ= v 0 ≃1, (3) ditions that lead to unstable modes excited by the κ- ΠL mechanism in our layer (cf. Paper I). where ∆m is the integrated mass between the consid- The two criteria (2) and (3) predict that modes are eredpointandthesurface,ΠthemodeperiodandLthe unstable for a “sufficient” hollow (adequate amplitude Numericalsimulationsoftheκ-mechanismwithconvection 3 and width) that is located in the transition region. Using the work integral computed from the obtained eigenvectors,weindeedrecoverinstabilitystripsforthe κ-mechanismbyconsideringthethreedifferentequilib- rium setups (Fig. 2): • The first case (dotted blue line, T /T = 1.7) bump surf corresponds to a hot star with an ionisation region close to the surface. As the surrounding density is small, Ψ≪1 and the conductivity hollow hardly in- fluencestheworkintegral: drivingisunabletoprevail over damping, leading to τ <0 or stable modes. • In the second case (solid green line, T /T = bump surf 2.1), the radiative conductivity begins to decrease significantly at the location of the transition region whereΨ≃1. Asaconsequence,drivingbecomesim- portant. Furthermore, the radiative conductivity in- creaseoccurswherenon-adiabatic effectsarealready significant, i.e. Ψ<1there. Itmeansthatnodamp- ing occurs between the hollow position and the surface astheradiative fluxperturbations arefrozen-in. Drivingisovercomingdampinginthiscase,therefore τ >0 and modes are unstable. Fig. 3 a) Temporal power spectrum for the momentum • The third case (dashed red line, T /T = 2.8) bump surf in the (z, ω)-plane. b) The resulting spectrum after in- corresponds to a cold star of which ionisation re- tegrating over depth. c) Comparison between normalised gionislocateddeeperinthestellaratmospherewhere momentum profiles for n = (0, 2) modes according to the Ψ ≫ 1. Ionisation then occurs in a quasi-adiabatic DNS power spectrum (solid black line) and to the linear location. As a consequence, the excitation provided stability analysis (dotted blue line). bytheconductivityhollowcannotbalancethedamp- ing arising at the top of the layer, thus τ < 0 and To determine which modes are present in the DNS modes are stable. in the nonlinear-saturation regime, we first perform In conclusion, our first radiative model of the κ- a temporal Fourier transform of the momentum field mechanism is well able to reproduce the main physics ρu(z,t) and plot the resulting power spectrum in the involved in the oscillations of Cepheid stars: (i) one (z,ω)-plane (Fig. 3a). The mean spectrum follows by should have dK /dz < 0 to drive the oscillations; (ii) integrating ρu(z,ω) over depth (Fig. 3b). With this T the thermal engine underlying the conductivity hollow method, acoustic modes are extracted because they c should also be located in the transition region where emerge as “shark-fin profiles” about definite eigenfre- Ψ≃1 (Fig. 2). quencies (Dintrans and Brandenburg 2004). Several discrete peaks corresponding to normal 2.3 Study of the nonlinear saturation modes well appear but the fundamental mode close to ω = 5.439 clearly dominates. Moreover, the lin- 0 To confirm both the results obtained previously in the ear eigenfunctions are compared to the mean profiles linear stability analysis and to investigate the result- computed from a zoom taken in the DNS power spec- ing nonlinear saturation, we perform direct numerical trumabouteigenfrequenciesω =5.439andω =11.06 0 2 simulations of our problem. The idea is to start from (Fig.3c). Theagreementbetweenthelineareigenfunc- the most favourable initial conditions found in the sta- tions (dotted blue lines) and the DNS profiles (solid bility analysis and then advance the general hydro- black lines) is remarkable. In summary, Fig. 3 shows dynamic equations in time. We use a modified ver- that several overtones are present in the DNS, even for sionofthepublic-domainfinite-differencePencilCode2 long times. These overtones are however linearly stable that includes an implicit solver of our own for the suggestingthatsomeunderlyingenergytransfersoccur radiative diffusion term in the temperature equation between modes through nonlinear couplings. (Gastine and Dintrans 2008b, hereafter Paper II). These nonlinear interactions are investigated thanks toatoolthatmeasuresthesoundgenerationbyturbu- 2Seehttp://pencil-code.googlecode.com/. lent flows (e.g. Bogdan, Cattaneo, and Malagoli 1993). 4 It involves projections of the DNS physical fields onto 3 Convective models in 2-D the regular and adjoint eigenvectors that are solutions to the linear oscillation equations. The projection co- 3.1 From radiative to convective setups efficients then give the time evolution of each acoustic After the previous results obtained on purely radia- mode that is present in the DNS separately. Further- tive models in 1-D, we next address the convection- more, as the evolution of the kinetic energy of each pulsation coupling that is suspected to quench the ra- mode is also accessible with this method, one can fol- dial oscillations of Cepheids close to the red edge. The low the energy transfer between modes. idea is to slightly modify the unstable radiative setup to obtain a convective zone superimposed to the ionisation region in 2-D simulations. The convective instability obeys Schwarzschild’s criterion given by (e.g. Chandrasekhar 1961) dT g gK (T ) 0 0 0 > =⇒ F > . (5) bot (cid:12) dz (cid:12) c c (cid:12) (cid:12) p p (cid:12) (cid:12) (cid:12)where(cid:12) F = K |dT /dz| is the imposed bottom flux. bot 0 0 For a given gravity g and conductivity hollow K (T ), 0 0 convectionwillthereforedevelopforalargeenoughbot- tom flux satisfying Eq.(5). Fig. 4 Time evolution of the kinetic energy ratio for the n=0 and n=2 acoustic modes. Figure 4 shows an example of the time evolution of the kinetic energy ratio E /Etot for the two modes n kin n = 0 (i.e. the fundamental one) and n = 2 (i.e. the second overtone). Here E denotes the kinetic energy n embeddedintheacousticmodeofordern,whileEtot is kin the total kinetic energy in the simulation. After its linear transient growth, agiven fraction of the fundamental mode energy is progressively transferred to the second overtone until the nonlinear saturation is achieved abovetimet≃150,withfinally85%ofthetotalkinetic Fig.5 Snapshotofthevorticityfield∇~ ×~vina2-Dmodel energyinthen=0modeandtheremaining15%inthe with convection and κ-mechanism. n=2 one. We recall that it is at a first glance surpris- ingtodetectthissecondovertoneatlongtimesbecause Figure5displaysasnapshotofthevorticityfieldob- it is linearly stable and should normally be damped by tainedinsuchaconvectivesimulationattimet=2500. diffusive effects. In order to ensure that the thermal relaxation is well Thereasonforitspresenceliesinafavored coupling achieved,weranthissimulationforamuchlongertime in the period ratio between these two modes. Indeed, than in the purely-radiative case and the resolution is the fundamental period is P0 = 2π/ω0 ≃ 1.155 in this 256×256 gridpoints. Asexpected, convection develops simulation while the n = 2 one is P2 ≃ 0.568. As a in the middle of the box where the conductivity hollow consequence, the corresponding period ratio is close to is located1. Moreover, because of the density contrast one half (P2/P0 ≃ 0.491) such that the second over- across the convection zone, the vorticity is trapped in tone is involved in the nonlinear saturation through a strong convective plume downdrafts that easily over- 2:1 resonance with the fundamental mode. The lin- shoot in the stable zone below due to their low Péclet ear growth of the fundamental mode is then balanced number (Dintrans 2009). by the pumping of energy to the stable second overtone that behaves as an energy sink, leading to the full 1An animation of this simulation is provided at limit-cycle stability. http://www.ast.obs-mip.fr/users/tgastine/conv_kappa.avi. Numericalsimulationsoftheκ-mechanismwithconvection 5 3.2 Evolution of acoustic modes with convection (dashed green line). Indeed, after a long transient during which it reaches a non-trifling value (∼ 30%), it To address the time evolution of the acoustic modes progressively tends to zero at long times while the con- thatarepresentintheDNSwithconvection,weproject vective motions are responsible for the whole kinetic as before the physical fields onto an appropriate acous- energy content. The radial oscillations excited by κ- tic subspace to get the projection coefficient c (t). mechanism are thus quenched by convection and this ℓn Here ℓ denotes the mode degree while n is its order situation is relevant to the red edge where the unstable (with k =(2π/L )ℓ and L is the width of the box). acoustic modes are known to be damped by convective x x x motions occuring at the star’s surface. We thus have shown that the Rayleigh number is a key parameter for this problem. But further inves- tigations of the convection-pulsation coupling are in progress, by especially checking time-dependent theo- riesofconvection forwhichseveralfreeparameterscan be constrained by such kind of direct numerical simulations (e.g. Yecko, Kolla´th, and Buchler 1998). Acknowledgements Calculations were carried out on the CALMIP machine of the University of Toulouse (http://www.calmip.cict.fr/). References Fig. 6 Evolution of the acoustic kinetic energy ratio for Baker,N.,Kippenhahn,R.: Astrophys.J.142,868(1965) twosimulationswithaweakconvection(solidblueline)and Baker, N.H.: In: Hillebrandt, W., Meyer-Hofmeister, E., a stronger one (dashed green line). Thomas,H.C.(eds.)PhysicalProcessesinComets,Stars and Active Galaxies, p. 105 (1987) Bogdan, T.J., Cattaneo, F., Malagoli, A.: Astrophys. J. Fig. 6 shows the total energy embedded in acoustic 407, 316 (1993) modes for two simulations with convection that differ Bono, G., Marconi, M., Stellingwerf, R.F.: Astrophys. J. by their Rayleigh number given by Suppl. Ser. 122, 167 (1999) Chandrasekhar, S.: Hydrodynamic and hydromagnetic stability. Oxford University Press: Clarendon, Oxford Ra=−gd4conv ds = (∇−∇adia)gd4conv, (6) (1961) νχc dz νχH Cox,J.P.: Theoryofstellarpulsation.PrincetonUniversity p p Press, Princeton (1980) where H = −(dlnp /dz)−1 is the pressure scale Cox, J.P., Whitney, C.: Astrophys. J. 127, 561 (1958) p 0 Dintrans,B.: CommunicationsinAsteroseismology158,45 height, d the vertical extent of the convection zone, conv (2009) χ=K /ρ c theradiativediffusivityandstheentropy. 0 0 p Dintrans, B., Brandenburg, A.: A&A 421, 775 (2004) Thisnumberallowstoquantifythestrengthofthecon- Eddington, A.S.: The Observatory 40, 290 (1917) vective motions, that is, the critical Rayleigh number Gastine, T., Dintrans, B.: Astron. Astrophys. 484, 29 above which convection develops is about 103 for poly- (2008a) Gastine, T., Dintrans, B.: Astron. Astrophys. 490, 743 tropic stratifications (e.g. Gough et al. 1976). (2008b) Inthefirstsimulation(Ra=4×104,solidblueline), Gough, D.O., Moore, D.R., Spiegel, E.A., Weiss, N.O.: As- the Rayleigh number is slightly super-critical and the trophys. J. 206, 536 (1976) convection intensity is weak. In this case, we see that Unno, W., Osaki, Y., Ando, H., Saio, H., Shibahashi, H.: the acoustic energy ratio remains large (i.e. & 80%, Nonradialoscillationsofstars.UniversityofTokyoPress, Tokyo (1989) the remaining 20% being in the convection) and the Valdettaro, L., Rieutord, M., Braconnier, T., Fraysse, V.: nonlinear saturation is reached in the same way than JCoAM 205(1), 382 (2007) in the previous radiative models. In other words, the Yecko, P.A., Kolla´th, Z., Buchler, J.R.: Astron. Astrophys. modepropagationisnotaffectedbyconvectivemotions 336, 553 (1998) and the convection can be seen as frozen-in. Zhevakin, S.A.: Annu. Rev. Astron. Astrophys. 1, 367 (1963) On the contrary, increasing the Rayleigh number leads to a more vigorous convection and the acoustic energyratioismuchmoreaffectedbyconvectiveplumes ThismanuscriptwaspreparedwiththeAASLATEXmacrosv5.2.