# [

## Abstract

We report on modelling in stellar astrophysics with the ANTARES code. First, we describe properties of turbulence in solar granulation as seen in high-resolution calculations. Then, we turn to the first 2D model of pulsation-convection interaction in a cepheid. We discuss properties of the outer and the HeII ionization zone. Thirdly, we report on our work regarding models of semiconvection in the context of stellar physics.

Stellar convection, pulsation and semiconvection] Simulations of stellar convection, pulsation and semiconvection Herbert J. Muthsam et al.] Herbert J. Muthsam, Friedrich Kupka, Eva Mundprecht, Florian Zaussinger , Hannes Grimm-Strele & Natalie Happenhofer 2010 \volumeIAU Symposium No. 271 \pagerange000–007 \jnameAstrophysical Dynamics: From Stars to Galaxies \editorsA.C. Editor, B.D. Editor & C.E. Editor, eds.

## 1 Introduction

In this paper we present applications of the ANTARES code to various problems in stellar physics. ANTARES is a radiation-hydrodynamics code. Its core functionality has recently been described by [Muthsam et al. (2010)]. Using this core functionality, we first discuss properties of turbulence in solar granulation, studied in extremely high resolution in 3D.

We have in the meanwhile extended the core functionality of ANTARES, and the applications discussed in the next two chapters make use of these extensions. So, we have included a spherical and optionally radially moving grid for the study of the interaction of radial stellar pulsation with convection. We present a first description of the pulsation-convection interaction of a 2D model of a cepheid with realistic equation of state and grey radiative transfer. For the importance of studying the pulsation-convection interactions consult, for example, the reviews [Buchler (1997)] and [Buchler (2009)]. Consider that since the pioneering work due to Christy, Cox and Kippenhahn in the 1960’s improvement in nonlinear stellar pulsation modelling has basically been due to better microphysics and to some extent to improvement in numerics. Convection has, however, (practically) always been included by a mixing-length type approach where already the basic form of the equations may be prone to some doubt and where additionally a number of parameters needs to be set in a way which can barely considered to be really convincing. Problems probably stemming from such shortcomings are discussed, e.g., in [Buchler (1997)]. For a recent discussion of a specific uncertainty which arises from traditional convection modelling in pulsating stars (namely, excitation of double mode pulsations or lack thereof) see [Smolec & Moskalik (2009)].

Subsequently, we turn our attention towards modelling of semiconvection in the context of stellar physics. Although it is textbook knowledge that a faithful description of semiconvection is highly important for properly dealing with the late stages of stellar evolution, but little work has been undertaken with regard to multidimensional modelling in the stellar context. The only major numerical studies are due to [Merryfield (1995)] and [Biello (2001)]. See also the results due to [Bascoul (2007)]. We present here our first 2D results discussing, in particular, the properties of diffusive interfaces. These results mainly refer to Boussinesq models, although we report also on compressible models and discuss their connection with the Boussinesq case. For enabling a more extensive investigation of compressible models it is necessary to implement a method which allows large time-steps, restricted by the macroscopic rather than the sound velocity. In this connection, an implementation of the method of [Kwatra et al.(2009)] is in an advanced stage of development.

## 2 High-resolution solar granulation simulations

In order to investigate the turbulent state of solar granulation we have performed hydrodynamic simulations with extremely high resolution. While they are basically similar to the ones described in [Muthsam et al. (2010)], they differ in that the highest resolution is achieved by making use of two grid refinement zones (stacked within each other) instead of one, and the resolution in the smaller one is considerably finer than in the old run. In addition, the old run referred to an exploding granule, whereas we concentrate here, in the highest resolution, on a strong downdraft surrounded by normal granules in order to figure out to what extent our previous results apply here as well.

More specifically, the basic computational domain spans now (the first coordinate is the vertical) with a cell size of . The first grid refinement zone has an extension of and a grid spacing of . The second and therefore finest grid refinement zone contains the downdraft mentioned above and parts of the surrounding granules. It has an extent of with a cell size of .

Using the high-resolution calculations, we want to discuss the role of resolution in the representation of the turbulent state of the solar granulation layer. In earlier simulations of an exploding granule [Muthsam et al. (2010)] a host of vortex tubes of small diameter has shown up, located preferentially near the granular lanes and the downdrafts. The questions popping up are now whether a similar large number of vortex tube shows up even in more normal granulation and whether these previous calculations were fully resolved; they used a grid spacing of , i.e. worse than even our present first refinement.

In order to answer the question at what resolution the turbulent field of solar granulation is (possibly) fully resolved, we compare in figure 1 two snapshots of our highest resolution subdomain. The first figure represents the state shortly after we have turned on the second grid refinement and thus still have the effective resolution due to the first refinement, whereas the next snapshot is taken seconds later. Due to increased resolution, the number and intensity of the vortex tubes has drastically increased. A closer study of the properties of this special sort of turbulence encountered mainly in granular downdrafts is under way.

## 3 Cepheid pulsation and convection: a 2D model

Initial condition and computational domain. We model a cepheid with an effective temperature of , luminosity , mass , hydrogen content and metallicity . Our computational domain reaches from to , thus the outer of the star are computed.

This domain is equipped with a polar grid. In radial direction 510 grid points plus 4 ghost cells at each boundary are used. The r-range covers where is fixed at and varies with time from about to . The grid is stretched in radial direction by a factor . The mesh sizes are = varying from at the top to at the bottom. In angular direction there are 800 grid points, the distance between two adjacent gridpoints in the H-ionisation zone is and in the He-ionisation zone. The corresponding aspect ratios are and .

Fluxes and pulsation in the HeII-ionisation zone. The work integral can be decomposed into an average part and a perturbational part defined as . Here, and . is the momentum vector and the velocity. denotes the horizontal average of the radial momentum component.

While varies locally from to in the HeII convection zone and is greater than , the horizontal average is considerably smaller, , compared to .

The averaged convective flux varies from to of the input energy flux at the lower boundary, the mean value during one period being . Locally values of up to corresponding to can be reached. The averaged kinetic flux varies from to the mean value is . The data are from 12 consecutive periods and the same flux pattern can be observed in each though the extent may vary. In the expanded state one observes two convection centers: a new one forming at the top and an older one which is the remainder of the downwards moving and slowly disapearing plume of the previous period. At contraction the two centers are very close to each other. These centers of convective activity lead also to the two (occasionally more) stripes situated atop of each other and visible, for example, in the quantities depicted in figure 3. The centers of activity can also directly be seen in figure 2.

Remarks on the H-ionisation zone.

The upper boundary of our computational domain is situated at an optical deph of a few thousandths. Figure 4 shows a part (in direction) of the upper zone. Color coding is for momentum densities (bright resp. red is downwards). The convective motions are very vigorous. We have observed Mach numbers for inflow of up to 4 and for outflow of up to 3. Remarkably, the lower boundary of the downdrafts often coincides the line of maximum temperature gradient (just above the bright line in figure 4. These results could only be achieved with grid refinement. The refined grid spacing is one third of the coarse grid spacing in the radial and one fourth in the angular direction.

## 4 Numerical simulation of semiconvection

When the convective transport does not depend on the temperature gradient only but also on an additional scalar field like salt (in the ocean) or helium (in stars), then double-diffusive convection can occur. This name implies that two diffusion coefficients ( and ) influence the fluids motion. In the oceans double-diffusive convection can form salt-fingers, while the well-known Latte Macchiato layers are a result of another double-diffusive phenomenon, which in astrophysics is called semi-convection.

In the evolution of stars semi-convection plays an important role. When cold hydrogen-rich matter is stratified above hot helium-rich matter, double-diffusive layering can occur under certain circumstances (Ledoux stable, ). The main question, which arises is the influence of the double-diffusive mixing processes on the evolution of high mass stars () around main sequence turnoff.

A stable layered structure in semi-convection zones (cf. Latte Macchiato), separated by diffusive boundary layers, is assumed ([Huppert & Moore (1976)], [Spruit(1992)]). In particular, the mass and heat transport through these diffusive boundaries are of major interest for explaining mixing time scales and the life span of the entire semi-convective region. There are several unknowns when describing single layers. Especially, the thickness of a single layer as function of thermal and saline Rayleigh number (, ) and the superadiabatic gradient is not well understood. A comprehensive hydrodynamical model, based on [Spruit(1992)] and [Muthsam et al. (1999)] (see also [Muthsam et al.(1995)]) was derived for the double diffusive convection case to tackle these uncertainties.

In [Zaussinger & Spruit(2010)] numerical simulations in 2D of semi-convection, based on compressible and incompressible formulations, have been performed for a wide initial parameter range of the Prandtl number , the Lewis number , the modified Rayleigh number and the stability parameter . It is proven there that simulations done in the Boussinesq approximation are directly comparable with fully compressible fluids as long as the layer heights are small enough (). The aim of this study was the verification and extent of validity of existing power laws in the form , where is the thermal Nusselt number. An extrapolation to the stellar parameter space was done subsequently.

In [Zaussinger(2010)] the parameter dependence of the fluxes is analysed. Single and double layer simulations were at the focus of that study. While the diffusion associated dimensionless numbers (, ) can be compared directly between compressible and incompressible fluids, other physical quantities like the fluxes (, ) or the Rayleigh numbers (, ) need a special treatment. The thermal Nusselt number, which compares the diffusive heat flux to the total heat flux, has to be corrected by the adiabatic stratification flux , which is by definition not present in incompressible fluids modelled by the Boussinesq approximation. Neglecting this difference would lead to wrong comparisons in simulations with flat temperature gradients.

Linear stability analysis shows that semi-convection in an oscillatory instability. This could be observed in all our simulations. Figure 5 illustrates the development of a semiconvection cell in terms of the Nusselt number when starting from a linear stratification. Nusselt numbers in the range indicate very diffusive regimes, where the total heat flux is dominated by conduction. High values correspond to more convective regimes. A diffusive phase (1) is followed by an oscillatory phase (2). The fully evolved semi-convective role (3) leads to the expected ([Spruit(1992)]) step like structure in helium and temperature at the upper and lower boundaries. With a suitable choice of parameters, a preassigned step in the middle of the domain is found to be stable for many eddy turnover times.

While adoption of linear stratification as initial condition leads to very long simulation times, a ‘direct initial step stratification’ is considered. It has been shown that both approaches lead to the same results. Strong damping mechanisms in the parameter space and lead to very long simulation times. By starting from a step like stratification these could be reduced significantly and a formerly inaccessible parameter range came into reach.

The verification of theoretical and experimental power laws under consideration for the stellar parameter space was the main focus of this study. It turned out that several existing ‘fitting formulas’ ([Castaing et al. (1989)], [Spruit(1992)], [Niemela et al. (2000)]) describe an upper limit. Even the flux relation,

(1) |

could be verified for semi-convective stable layers (see figure 6). Equipped with these results an extension to the stellar parameter regime was considered. The main problem of existing 1D stellar evolution codes is the lack of information about to the semi-convective mixing efficiency. The superadiabatic gradient is not known in advance and has to be estimated. The semi-convective zone in a stellar model is limited in extent and consequently the evolution of the star on small time-scales does not depend much on the way it is calculated. But the additional diffusion in helium could be important for later evolutionary stages. This leads to the assumption that the heat flux could be considered as known rather than the real temperature gradient . Based on [Spruit(1992)] an implicit function was derived, relating the temperature gradient, the Rayleigh number, the fluxes and the the vertical extent of the semi-convective zone to each other.

Acknowledgements This research has been supported by grants P18224, P20762, P20973 and P21742 of the Austrian Science Foundation and by grant KU 1954/3 within SPP 1276 of the Deutsche Forschungsgemeinschaft. Calculations were performed within DEISA (project SOLEX), at the Leibniz Computer Center, Munich (project SOLARSURF), at computers of the Max Planck Society and at the VSC cluster, Vienna. We are grateful to G. Houdek, Vienna, for providing us with starting models for cepheids.

Trampedach Your simulation of semi-convection looked 2-D. Is that correct? Will it be extended to 3-D?

MuthsamYes, the present models are 2D. We will turn to 3D in the forseeable future.

Gough The stated purpose of the final solution you presented was to determine whether the interface between turbulent convection layer is purely diffusive, but you did not tell us whether that is the case. Please would you tell us whether the transport of heat and helium across the interface which you showed us, which maintains its integrity despite being modulatory is purely diffusive or whether there is significant (small scale) convective transport?

Muthsam You obviously refer to the movie with the pre-assigned step in the middle of the domain. The interface itself looks quite stable. Only on isolated spots (upstreams or downstreams) these streams can penetrate it somewhat.

PouquetDo you see any advantage to use a hybrid (MPI/OpenMP) parallelization method?

Muthsam For really large runs you have to use MPI anyway for lack of shared memory. If you have, on one node, multiple cores you save messages when using them in the OpenMP rather than in the MPI mode.

### References

- Bascoul, G. P. 2007, IAU Symposium, 239, 317
- Biello, J. A. 2001, Ph.D. Thesis
- Buchler, J. R. 1997, Variables Stars and the Astrophysical Returns of the Microlensing Surveys, 181
- Buchler, J. R. 2009, American Institute of Physics Conference Series, 1170, 51
- Castaing, B., Gunaratne, G., Kadanoff, L., Libchaber, A., & Heslot, F. 1989, Journal of Fluid Mechanics, 204, 1
- Huppert, H. E., & Moore, D. R. 1976, Journal of Fluid Mechanics, 78, 82
- Kupka, F., Ballot, J., & Muthsam, H. J. 2009, Communications in Asteroseismology, 160, 30
- Kwatra, N., Su, J., Grétarsson, J. T., & Fedkiw, R. 2009, Journal of Computational Physics, 228, 4146
- Merryfield, W. J. 1995, ApJ, 444, 318
- Mundprecht, E. 2010, PhD thesis, University of Vienna
- Muthsam, H. J., Göb, W., Kupka, F., Liebich, W., & Zoechling, J. 1995, A&A, 293, 127
- Muthsam, H. J., Göb, W., Kupka, F., & Liebich, W. 1999, NewAst, 4, 405
- Muthsam, H. J., Kupka, F., Löw-Baselli, B., Obertscheider, C., Langer, M., & Lenz, P. 2010, NewAst, 15, 460
- Niemela, J. J., Skrbek, L., Sreenivasan, K. R., & Donnelly, R. J. 2000, APS Meeting Abstracts, 2
- Smolec, R., & Moskalik, P. 2009, American Institute of Physics Conference Series, 1170, 73
- Spruit, H. C. 1992, A&A, 253, 131
- Zaussinger, F. 2010, PhD thesis, University of Vienna
- Zaussinger, F., & Spruit, H. 2010 in preparation for A&A