{'fraction_non_alphanumeric': 0.025608732157850547, 'fraction_numerical': 0.014903442485306465, 'mean_word_length': 4.779619894864537, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
3D Anelastic Simulations of Convection in Massive Stars

After briefly describing the anelastic approximation and Glatzmaier's code, we present results from our preliminary studies of core convection during the hydrogen burning phase of a 15 solar mass star, as well as our most recent results concerning convection in the oxygen shell of a 25 solar mass star.


Introduction
Two-dimensional studies of turbulent convection in the interiors of massive stars have been quite common (Deupree 2000;Asida & Arnett 2000;Bazan & Arnett 1998), but the problem is inherently three-dimensional and many aspects of the convective flow will not be reproduced correctly in two-dimensional models (Canuto 2000). Only recently has it become possible to begin tackling the problem of simulating convection in all its three-dimensional g(l)ory, and several major projects are under way to do just this (Dearborn et al. 2001;Brun 2002;Porter & Woodward 2000). The price of this additional dimension is additional computational complexity, the need for faster computers and more data storage capabilities, not to mention difficulties in visualizing and analyzing the resulting data. In this paper we describe three-dimensional simulations based on anelastic hydrodynamics, which we hope capture the essential phenomena of the problems at hand, while allowing more rapid progress at a lower computational cost than comparable studies using fully compressible codes. This paper is divided into three parts. After very briefly describing the anelastic approximation and Glatzmaier's implementation of it in Section 2, we present our hydrogen core convection model and results of this simulation in Section 3, and the oxygen shell convection model and results in Section 4.
In addition to the static plots shown here, please look at http://www.supersci.org for movies of these models.

The Numerical Method
The simulations presented in this paper were performed using anelastic hydrodynamics. This approximation allows a background density stratification, but filters out short-period, acoustic oscillations, thereby allowing larger timesteps 2 Kuhlen, Woosley, & Glatzmaier than fully compressible codes. We believe that these attributes make anelastic hydrodynamics an excellent method for studying stable convective regions in massive stars.

The Anelastic Approximation
For completeness we have included the anelastic equations below. Any interested reader can find a much more detailed exposition of the anelastic approximation in Glatzmaier (1984) and references therein.
Barred symbols denote reference state quantities, all others perturbations.
Here E is the rate of strain and I the identity tensor. A vanishing coefficient of bulk viscosity is assumed for the viscous force calculation, and viscous heating is neglected.
Of particular interest is the mass conservation equation in which the Eulerian time derivative of the density perturbation drops out. This in effect filters out sound waves and greatly ameliorates the Courant condition, which in turn allows larger timesteps. Note that by design anelastic codes cannot model any shocks, or (super)sonic flow in general. The anelastic approximation attempts to solve for an average solution of the flow, ignoring the propagation of any small scale pressure fluctuations.

Glatzmaier's Implementation
The numerical implementation of these equations has been performed by Glatzmaier, and a detailed description of the code can be found in Glatzmaier (1984). The code is fully spectral, with an expansion in spherical harmonics to cover the angular variations and in Chebyshev polynomials for the radial dependence.
The models described here have no radiative energy transport, so there are only two ways to transport energy: convective and diffusive heat flux. The nature of the diffusive heat flux is explained below in Section 2.5.

The Reference State
One of the most important steps in producing these simulations is the generation of the reference state. We have incorporated a general equation of state and have added temperature and density dependent energy generation and neutrino losses.

3D Anelastic Simulations of Convection in 3D
3 Radial Structure As the reference state is one-dimensional, we take advantage of our 1D stellar evolution code KEPLER (Weaver, Zimmerman, & Woosley 1978) and fit its density and pressure to a polytrope (P = Kρ 1+1/n ). For the radial dependence of these quantities we solve the resulting Lane-Emden equation.

Equation of State
The background temperatureT as well as the necessary thermodynamic derivatives ∂ρ ∂s P , ∂ρ ∂P s are determined using the publicly available Helmholtz code (Timmes & Swesty 2000).

Energy Generation and Neutrino Losses
It is important to capture the full temperature and density dependent nature of the energy generation within the convective region, and in the case of the oxygen shell also of the neutrino losses. For this reason we have added not only radially varying reference state nuclear energy generation and neutrino loss rates, but have also kept the perturbations to allow for differential heating between cold and hot convective bubbles. The nuclear energy generation rate in the relevant zones is fit to the conventional ǫ nuc ∝ ρT n law, and the neutrino loss rate to ǫ ν ∝ ρ m T n .

Boundary Conditions
The top and bottom boundaries are modeled as impermeable and stress-free (v ⊥ = 0 and ∂ ∂r v r = 0). This means that presently we are not able to study important phenomena such as convective overshoot.
Generally we fix the entropy gradient at the bottom boundary by specifying a base luminosity, taken from KEPLER. For the boundary condition at the top we either force the entropy perturbation to vanish or we specify the top luminsity, i.e. the entropy gradient at the outer boundary.

Turbulent Diffusion
The diffusion of momentum and entropy by small-scale eddies is expected to dominate the diffusion due to radiative and molecular processes. The effects of all unresolved motions are parametrized as turbulent diffusion. In this model we have chosen constant turbulent viscous and thermal diffusivities, which we have tried to minimize to the extent permitted by the finite resolution.

Hydrogen Core Convection
As our initial testcase for this code, in order to get familiar with it, we chose to model core convection of a 15M ⊙ star half-way through the main sequence. This stage of the star's life is not very violent; there are no neutrino losses, the energy generation doesn't have a very steep temperature dependence, there are no shocks, and density and pressure contrasts are expected to be small -a perfect case for an anelastic model. We view this model as a proof of principle, and as preparation for modeling more advanced stages of stellar evolution 15 M ⊙ hydrogen core convection model reference state.

Parameters
The dependence of the reference state radius, density, temperature, and pressure on the Lagrangian mass coordinate are shown in Figure 1. Table 1 below summarizes the remaining parameters of interest. M base is the mass of the cen- Table 1.
Parameters of the 15M ⊙ model of core convection during main sequence hydrogen burning.
tral portion of the star that is cut out, and M top is the total mass within the convective core of this 15M ⊙ star. Even without the central region we are still modeling 99.8% of the convective core, both by volume and mass. Table 2.

Results
Numerical results of 15M ⊙ hydrogen core convection simulations. Resolution: n max = 241, l max = 63, m max = 31  and hence about 6 convective turnovers. The first row of Table 2 summarizes the numerical results of this simulation, and Figure 2 provides some visualization of the convective motions. Note that without rotation there is no preferred orientation of the model, which explains why figures 2a) and 2b) look so similar. Both slices seem to indicate a somewhat dipolar flow field, with high entropy gas rising on one side and falling at the other. Such dipolar flow pattern have been observed in fully compressible calculations as well (Jacobs, Porter, & Woodward 1998).
All thermodynamic contrasts ((δρ ′ /ρ) max , (δP ′ /P ) max , etc.) are of the order 10 −5 , comfortably below 1%. The mass-averaged reference state sound speed is about 700km/s, so the flow is completely subsonic. We conclude that this problem is very suited for an anelastic treatment, as expected.
The Rotating Case The rotating model was run for about 120 days in the star's life, with the same resolution as the non-rotating case. This amount of star time corresponds to about 15 convective turnovers. As is immediately obvious by looking at Figure 3, rotation has a significant effect on the convective flow pattern. The rotational axis breaks the symmetry of the problem, so the flow behaves differently in the equatorial plane than in any meridional plane and the dipolar flow pattern is no longer visible. In the meridional slice the strong, high kinetic energy flows are concentrated close to the axis, with convection being strongly suppressed in the equatorial regions further away from the axis. We also see evidence of Taylor columns, elongated, columnar flow structures parallel to the rotation axis.

Oxygen Shell Convection
Of the advanced stages of stellar evolution the oxygen shell burning phase is of particular importance. For one, convection during this brief time before core collapse might produce an asymmetry, leading to non-uniform mixing and perhaps even an asymmetric supernova explosion. Several two-dimensional studies of this problem have been conducted (Arnett 1994;Bazan & Arnett 1998;Asida & Arnett 2000), and several of these simulations found close to sonic flow, large density and pressure contrasts, and hence significant acoustic fluxes.
As an alternative to these simulations, we present results of our anelastic model of a similar problem -oxygen shell convection in a 25M ⊙ star, shortly after core oxygen depletion. As shown below, we have found a solution within the anelastic approximation, with small density and pressure contrasts, and subsonic flows. Table 3.

Parameters
Parameters of the 25M ⊙ model of convection during oxygen shell burning, about 4 days before core collapse. During the modeled phase the rate of nuclear energy generation exceeds the rate of neutrino losses by about a factor two. The net energy generation rate is positive, top base 4πr 2 ρ(r)(ǭ nuc (r) +ǭ ν (r)) = 1.5 × 10 45 erg/s. KEPLER indicates that this energy is going into P dV work and into change of the gravitational potential. In the models presented here we have chosen to allow this excess energy to escape off the top, as we are currently unable to model expansion. Table 4.
non-rotating 2 × 10 −3 1 × 10 −3 6 × 10 −4 7 × 10 −4 rotating 3 × 10 −3 3 × 10 −3 1 × 10 −3 1 × 10 −3 The Non-rotating Case The non-rotating model was run for an equivalent star time of 6500 seconds. Given an average speed of about 50 km/s we determine the convective turnover time to be about 70 seconds, which implies that we followed    Figure 5 shows that the convective speed remains around two orders of magnitude below the sound speed throughout the entire modeled region. All thermodynamic contrasts are below the 1% level, as shown in Table 5. It is important to note that nothing in the model forces the resulting convective velocities to be subsonic and nothing prevents the thermodynamic perturbations from rising to a level comparable to the background. The fact that this solution remained subsonic and with small thermodynamic contrasts implies to us that the anelastic approximation is a valid one for the problem of the convective oxygen shell.
The entropy perturbations in the equatorial plane (see Figure 6) indicate a flow pattern dominated by four large cells. High entropy fluid is rising towards 90 • and 270 • , and low entropy fluid is falling along the 0 • and 180 • directions. We have not yet determined how much of this simple structure is dependent on Rayleigh and Reynolds numbers. It is possible that higher resolution simulations, with higher Rayleigh and Reynolds numbers, will destroy this pattern.
The Rotating Case The rotating case was run for 3300 seconds, which corresponds to about 50 convective turnovers. Again the convective speed is very subsonic and the thermodynamic contrasts are still below 1%, albeit somewhat higher than in the non-rotating case.
Although the reference state is rotating with a fixed angular speed (Ω = 0.08 rad/s) the flow is free to arrange itself into any form of differential rotation. The result can be seen in Figure 5. Our rotating model shows a lot of differential rotation, with the outer regions of the equatorial plane rotating nearly twice as fast as the inner ones. Just like in the hydrogen core convection model we see Taylor-column-like structures parallel to the rotation axis in the meridional planes.

Conclusions
Anelastic hydrodynamics appears to be a useful tool for studying convection in the interiors of massive stars in three dimensions. As long as the problem at hand satisfies the anelastic conditions of subsonic flow speeds and small ther-modynamic perturbations, a lot of progress can be made with an anelastic code in a shorter amount of time and with less computational resources than with a fully compressible code.
In this paper we have presented results of our efforts to model convection in the core of 15M ⊙ main sequence star as well as in the oxygen burning shell of a 25M ⊙ star at a much more advanced stage of stellar evolution. The resulting convective speeds are in agreement with results from KEPLER using mixing length theory. The resulting three dimensional flow pattern are not always uniform, and can show strong asymmetries, such as the dipolar flow in the hydrogen core convection model. We have found a significant difference in the flow pattern of our rotating models compared to the non-rotating ones. The presence of rotation breaks the spherical symmetry and allows more vigorous convection along the axis of rotation. The rotating oxygen shell model showed a significant degree of differential rotation.
This work was supported by the Scientific Discovery though Advanced Computing (SciDAC) Program of the DOE (DE-FC02-01ER41176).