3D simulations of Rayleigh-Taylor mixing in core-collapse SNe with CASTRO
aa r X i v : . [ a s t r o - ph . H E ] S e p Draft version June 4, 2018
Preprint typeset using L A TEX style emulateapj v. 03/07/07
3D SIMULATIONS OF RAYLEIGH-TAYLOR MIXING IN CORE-COLLAPSE SNE WITH CASTRO
C.C. Joggerst,
A. Almgren, , S. E. Woosley Draft version June 4, 2018
ABSTRACTWe present multidimensional simulations of the post-explosion hydrodynamics in three different15 M ⊙ supernova models with zero, 10 − Z ⊙ , and Z ⊙ metallicities. We follow the growth of theRayleigh-Taylor instability that mixes together the stellar layers in the wake of the explosion. Modelsare initialized with spherically symmetric explosions and perturbations are seeded by the grid. Calcu-lations are performed in two-dimensional axisymmetric and three-dimensional Cartesian coordinatesusing the new Eulerian hydrodynamics code, CASTRO . We find as in previous work, that Rayleigh-Taylor perturbations initially grow faster in 3D than in 2D. As the Rayleigh-Taylor fingers interactwith one another, mixing proceeds to a greater degree in 3D than in 2D, reducing the local Atwoodnumber and slowing the growth rate of the instability in 3D relative to 2D. By the time mixing hasstopped, the width of the mixed region is similar in 2D and 3D simulations provided the Rayleigh-Taylor fingers show significant interaction. Our results imply that 2D simulations of light curves andnucleosynthesis in supernovae (SNe) that die as red giants may capture the features of an initiallyspherically symmetric explosion in far less computational time than required by a full 3D simulation.However, capturing large departures from spherical symmetry requires a significantly perturbed explo-sion. Large scale asymmetries cannot develop through an inverse casacde of merging Rayleigh-Taylorstructures; they must arise from asymmetries in the initial explosion.
Subject headings: INTRODUCTION
Supernova 1987A has furnished astronomers with per-haps the best opportunity to study a core-collapse ex-plosion to date. One of the most exciting discoveries toemerge from this event is the evidence for large-scale,extensive mixing in the supernova ejecta. γ − ray linesemitted by the decay of Co were detected half a yearearlier than expected for a spherically symmetric ex-plosion model (Matz et al. 1988). Modeling the lightcurve in 1D requires a large amount of Ni to be mixedoutward and H and He to be mixed inward (Utrobin2004). The Bochum event, in which fine-structure H α lines are observed 2 weeks after the explosion, is ex-plained by the ejection of 10 − M ⊙ of Ni moving witha velocity in excess of 4000 km s − into the far hemi-sphere of 1987A (Hanuschik et al. 1988). The Rayleigh-Taylor (RT) instability was posited as a possible explana-tion for the large-scale mixing implied for this explosion.Many groups using both SPH and grid-based codes sim-ulated the evolution of 87A-like progenitor models, firstin 2D, then in 3D. It soon became apparent that morethan just the Rayleigh-Taylor instability is needed to ex-plain the high Ni velocities as well as the large degreeof outward mixing of heavy elements and inward mix-ing of lighter elements. Simulations that posit modestinitial perturbations are unable to replicate these fea-tures of the explosions. Kifonidis et al. (2006) find thatby following the explosion from early times, and usinglarge, low-order perturbations in the inner layers, theyare able to replicate the high Ni velocities and large Department of Astronomy and Astrophysics, University of Cal-ifornia at Santa Cruz, Santa Cruz, CA 95060; [email protected] Theoretical Astrophysics (T-2), Los Alamos National Labora-tory, Los Alamos, NM 87545 Computational Research Division, Lawrence Berkeley NationalLab, Berkeley, CA 94720 amounts of mixing seen in the explosion. The stand-ing accretion shock instability (SASI) (Blondin et al.2003; Blondin & Mezzacappa 2006; Scheck et al. 2006;Burrows et al. 2006, 2007; Marek & Janka 2009) pro-vides a mechanism by which the inner parts of the ejectacould be deformed.Most studies of the Rayleigh-Taylor instability in corecollapse supernovae use 1987A-like progenitor models.However, supernovae in the high-redshift universe, aswell as many modern day supernovae which die asred, not blue, giants, explode with presupernova struc-tures which differ significantly from 1987A. This maygreatly influence the evolution of the Rayleigh-Taylor in-stability in their outer layers. Simulations of rotatingzero-metallicity supernovae (Ekstr¨om et al. 2008; Hirschi2007) indicate that they may die as extended red, ratherthan compact blue, giants. Rotation mixes CNO ele-ments to the base of the hydrogen-burning shell, wherethey catalyze the CNO cycle, greatly boosting the burn-ing rate of H. This leads to the formation of a large con-vective envelope, effectively turning a compact blue starinto a giant red one. Models indicate that the shell boostinduced convection effectively mixes together the heliumshell and the hydrogen envelope. The dense helium shellhas been found to be important to the post-bounce evo-lution of 1987A-like progenitors (Hammer et al. 2009).Rotation has less of an effect on early stars with verylow, but not zero, metallicity, such as the Z = 10 − Z ⊙ progenitor models studied in Joggerst et al. (2010)as well as the current work. These stars remain blue.Solar metallicity stars end their lives as red giants, butretain their helium shells. Red stars should evolve differ-ently from 1987A-type models after the supernova shockis launched. They are larger, so the shock takes longerto traverse the star, and the Rayleigh-Taylor instabilitieshave more time to develop (Chevalier 1976).Previous simulations of core-collapse explosions intwo dimensions (Arnett et al. 1989; Hachisu et al. 1990;Fryxell et al. 1991; Herant & Benz 1991; Mueller et al.1991; Herant & Benz 1992; Hachisu et al. 1992) have pri-marily focused on 87A-like progenitor models. Studies ofRayleigh-Taylor-induced mixing in three dimensions hasexclusively focused on attempts to replicate observationsof 1987A, and have found significant differences betweensimulations performed in two and three dimensions. Themost successful of these 2D attempts, Kifonidis et al.(2003, 2006) found that to reproduce the high Ni ve-locities observed in 1987A large initial asymmetries oflow mode order were necessary. In studies of mixingin both red and blue supergiants(Joggerst et al. 2009;Herant & Woosley 1994) mixing is found to be more vig-orous in red stars than in blue stars. Joggerst et al.(2010) evolve 36 models spanning three rotation rates,three explosion energies, three masses, and two metal-licities in two dimensions, finding that mixing is morevigorous and goes on longer in stars that die as red asopposed to blue giants, in higher energy explosions, andin less massive stars. The literature on 3D simulationsis more sparse. Kane et al. (2000) find that for a single-mode perturbation of 10% amplitude, Rayleigh-Taylorfingers grow 30 −
40% faster in 3D than in 2D. Themost recent paper on the post-explosion hydrodynamicsin 1987A (Hammer et al. 2009) use a successful explosionmodel from Scheck (2007) to initialize a 3D simulation,along with several 2D simulations derived from merid-ional slices through the initial explosion model. Theyfind that for this model, which is highly asymmetric inthe inner regions, the Rayleigh-Taylor instability growsfaster and ballistically moving clumps of Ni reach thehydrogen envelope in 3D. These clumps are effectivelystopped by the dense He shell in the suite of 2D modelsinitialized from slices through the Scheck (2007) model.This is the first paper to simulate the post-explosionhydrodynamics of models with a diversity of presuper-nova structures in three dimensions. All models havethe same mass (15 M ⊙ ) and explosion energy (1 . x ergs at infinity), but two models explode as red giantsand one as a more compact blue star. The red giantwith zero metallicity effectively lacks a He shell becauseof convection arising from the hydrogen shell boost. TheZ = 10 − Z ⊙ and Z= Z ⊙ models both have helium shells,though the former is blue and the latter red at the timeof explosion. This diversity of models allows us to see ifpresupernova structure has an impact on the subsequentevolution of nearly isotropic 2D as opposed to 3D mod-els. These calculations set the stage for further studyof asymmetric explosions in these diverse presupernovamodels.The numerical methods employed in this study are dis-cussed in Section §
2, and initial models and problemconfiguration are discussed in Section §
3. Results arepresented in Section §
4, and are discussed and comparedwith previous 3D simulations in §
5. We offer some con-clusions in § NUMERICAL ALGORITHMS
The multidimensional simulations described in this pa-per are performed with
CASTRO (Almgren et al. 2010),an Eulerian adaptive mesh refinement (AMR) hydrody-namics code. Time integration of the hydrodynamics equations in
CASTRO is based on an unsplit higher-orderpiecewise parabolic method (PPM); the additional fea-tures are described below. All simulations described hereused either axisymmetric coordinates in 2D or Cartesiancoordinates in 3D.
Equation of State
We followed sixteen elements, from hydrogen through Ni, so that our elemental abundances could later beused to compute spectra. The atomic weights andamounts of the elements are used to calculate the meanmolecular weight of the gas required by the equation ofstate. The abundance of Ni was followed separatelyin order to calculate the energy deposited by radioactivedecay of Ni and Co.The equation of state in our simulations is as describedin Joggerst et al. (2010). It assumed complete ionizationand included contributions from both radiation and idealgas pressure: P = f ( ρ, T ) 13 aT + k B T ρm p µ (1) E = f ( ρ, T ) aT ρ + 1 . k B Tm p µ , (2)where P is the pressure, a is the radiation constant, k B is Boltzmann’s constant, T is the temperature, ρ is thedensity, m p is the proton mass, µ is the mean molecularweight, and E is the energy.The function f ( ρ, T ) is a measure of the contributionof radiation pressure to the equation of state that ac-counts for the fact that radiation ceases to be trapped atsome point. It is 1 in regions where radiation pressureis important, i. e. where gas is optically thick, and 0in regions where radiation pressure is unimportant, i. e.where gas is optically thin, with a smooth transition inbetween. The function f ( ρ, T ) takes the form f ( ρ, T ) = ρ ≥ − gm cm or T ≤ T neg f ( T ) = e Tneg − TTneg if ρ < − gm cm and T > T neg where T neg is the temperature at which contributionsto the pressure from radiation are 100 times less thanthat contributed by ideal gas pressure: T neg = 3 k b ρ m p µa / . (3)Radiation pressure will begin to dominate the equa-tion of state above T neg without some adjustment, eventhough this is unphysical, as radiation pressure is neg-ligible in optically thin regions. Damping the radiationcomponent of the EOS with the function f ( ρ, T ) in theseregions provides a more physical solution. Radioactive Decay of NiEnergy from the radioactive decay of Ni to Fe wasdeposited locally at each mesh point, in the same manneras described in Joggerst et al. (2010). The rate at whichenergy from the decay of Ni to Co was deposited inthe grid is given by: dE Ni = λ Ni X Ni e − λ Ni t q ( N i ) . (4)The decay rate of Ni, λ Ni , is 1 . × − s − , and theamount of energy released per gram of decaying Ni is q ( N i ), which we set to 2 . × erg g − . X Ni is themass fraction of Ni in the cell. The mass fraction of Co at a given time can be expressed in terms of themass fraction of initial Ni by X Co = λ Ni λ Co − λ Ni X Ni ( e − λ Ni t − e − λ Co t ) , (5)so that the energy deposition rate from Co is dE ø = λ Ni λ Co − λ Ni X Ni ( e − λ Ni t − e − λ Co t )) λ Co q ( Co ) . (6)We assumed a decay rate λ Co = 1 . × − s − andan energy per gram of decaying Co, q ( Co ), equal to6 . × erg g − . Gravity
Although
CASTRO supports several different approachesto solving for self-gravity, we used the monopole approx-imation for gravity in the calculations presented here, asdiscussed in Joggerst et al. (2010). Because the densityprofiles in the simulations depart very little from spher-ical symmetry, the gravitational field constructed usingthe monopole approximation is nearly identical to thatfound by solving the full Poisson equation for the grav-itational potential, and using the monopole approxima-tion significantly reduces the computational cost of thecalculations.Gravity from a point mass located at the origin wasalso included in the gravitational potential. The pointmass represents the compact remnant left behind bythe SN explosion. As infalling matter crosses the zero-gradient inner boundary near the origin, it is added tothis point mass.
Adaptive Mesh Refinement
CASTRO ’s AMR algorithm uses a nested hierarchy oflogically-rectangular grids with simultaneous refinementof the grids in both space and time. The integration al-gorithm proceeds recursively, advancing coarse grids intime, then advancing fine grids multiple steps to reachthe same time as the coarse grid, and finally synchroniz-ing the data at different levels.The regions of refinement evolve throughout the simu-lation based on user-specified refinement criteria appliedto the solution. In the simulations presented here, re-gridding of all grids at level ℓ + 1 and above occurredevery two level ℓ time steps.Our refinement criteria are based on the “error esti-mator” of L¨ohner (1987), which is essentially the ratio ofthe second derivative to the first derivative at the pointat which the error is evaluated. Details of the implemen-tation are discussed more fully in Almgren et al. (2010).The result is a dimensionless, bounded estimator, whichallows arbitrary variables to be subjected to the samecriterion for error estimation. In these calculations, den-sity, pressure, velocity, and the abundances of Ni, He,and O were used as refinement variables, with the el-emental abundances only used in a particular region iftheir abundance was greater than 10 − . SIMULATIONS SETUP AND EXECUTION
The simulations presented here were carried out in twodistinct stages. First each stellar model was evolved inone dimension using the code
KEPLER to the point wherethe core became unstable to collapse. It was then artifi-cially exploded by means of a piston located at the baseof the oxygen shell with enough energy that the explo-sions had 1 . × ergs of energy at infinity. The modelswere evolved forward in one dimension for 20 seconds formodels z15 and u15, and 100 seconds for s15, until allnuclear burning had ceased.These profiles were then mapped onto a two-dimensional r-z or three-dimensional Cartesian grid in CASTRO , and the calculation was evolved until the shockreached the edge of the simulation domain. The calcula-tion was then stopped, the domain enlarged, and the cal-culation restarted within the larger domain (see § ⊙ models, at zero, 10 − Z ⊙ , andsolar metallicities. Progenitor Models
Three different supernova models, representing threedifferent metallicities, are presented in this paper. Modelz15 has zero initial metallicity, and represents a Popula-tion III (Pop III) star. Model u15 has Z = 10 − Z ⊙ ,and model s15 has solar metallicity. All models are de-rived from rotating 15 M ⊙ progenitors. The zero and lowmetallicity models presented in this paper are the sameas presented in Joggerst et al. (2010) and Zhang et al.(2008). The solar metallicity model, model s15, is sim-ilar to the model presented in Woosley & Heger (2007).For moderate values of rotation, the amount of rotationwas found to have little effect on its post-explosion evo-lution (Joggerst et al. 2010). There are, however, pro-found differences between rotating and non-rotating zerometallicity supernovae. Several studies (Ekstr¨om et al.2008; Hirschi 2007; Joggerst et al. 2010)) have found thata small amount of rotation in zero metallicity massivestars dredges up enough carbon, nitrogen, and oxygenfrom the helium burning shell to catalyze a CNO cy-cle at the base of the hydrogen envelope. This leads torapid burning and a shell boost that triggers convection,mixing most of the helium and hydrogen shells together,and “puffing up” the outer envelope of the star. This isapparent in Figure 1. Models with a small amount ofmetals (i. e. u15 in this paper) do not show this shellboost behavior, and die as blue supergiants. Rotationalso has little effect on the presupernova structure of thesolar metallicity model s15, which dies as a red giant. Initialization of Multidimensional Data
In mapping the radial data from
KEPLER onto the mul-tidimensional grids in
CASTRO , special care was taken toproperly resolve the key elements of the simulations: theshock, the elemental shells, and the Ni at the centerof the explosion. In particular, both the Ni and theO shell were resolved with a minimum of 16 cells at
Fig. 1.—
Shell structure for H, He, O, Si, and Ni for our initial models. Other elements were included in our simulations but are left outof the figure for clarity. Note that while s15 and u15 (top and middle figures) show an He shell, this shell is not present in model z15 becauseof convection arising from a hydrogen shell boost. The models are 15 solar masses, but the outer layers are essentially a continuation ofthe levels at the right hand side of the plot, and have been omitted to better show the detail in the center. the highest level of refinement. This mapping resultsin explosions that are spherically symmetric in
CASTRO ;low-order departures from spherical symmetry are sup-pressed. Our models therefore only capture higher-orderasymmetries in the explosions.Perturbations that give rise to the Rayleigh-Taylor in-stabilities are seeded from the axisymmetric or Cartesiangrid. We performed calculations in 2D at twice the res-olution of the simulations considered in this paper andcompared the results of the high resolution simulations tothe low resolution simulations. The width of the mixedregion and velocities of the elements were essentially thesame, indicating that our results are numerically con-verged and do not depend on the initial scale of pertur-bations imposed by the grid.
Enlarging the Domain
In order to minimize the amount of computational re-sources used, we implemented a strategy of enlarging thesimulation domain only as necessary as the size of the re-gion of interest increased in time. The simulations wereinitialized on a 128 n grid with 2 levels of refinement,where n is the dimensionality of the simulation. Thisgives an effective grid resolution of 512 n at the finestlevel. A given simulation was advanced in time to thepoint where the shock was near the edge of the domain.The simulation was then stopped and the domain wasdoubled in each coordinate dimension. The enlarged do-main was covered by grids with cell spacing twice that of the previous coarsest resolution, so that after each en-largement the base grid continued to be 128 n . The sim-ulation then had 3 levels of refinement; in most cases(except for u15) the highest level of refinement was re-moved at this point, reducing the finest resolution by afactor of 2 as well. This allowed us to tailor the size ofthe domain to the size of the region of interest, withoutwasting computational resources in regions where noth-ing of interest was happening. Each time the domain wasenlarged data from the new regions of the domain wereinterpolated from the original model.To ensure that the strategy described above did notaffect the results of the simulations relative to analogouscalculations performed on domains which were largerthroughout, we ran two-dimensional simulations in whichthe final domain size was imposed from the start. Theresults were essentially the same, although more complexmixing structures are visible in the case with the initiallylarger domain. The width of the mixed region remainsthe same, i.e. heavy elements are mixed out the thesame point in radius and mass, while lighter elementsare mixed inward to the same point in the both cases.The velocities obtained by the elements are the same.This is sufficient evidence that the strategy of enlargingthe domain throughout the calculation did not introducespurious artifacts. RESULTS
In each multidimensional simulation presented here,when the shock encountered the large region of increas-ing ρr at the base of the hydrogen envelope of thepresupernova star (or, in the case of model z15, thehelium-hydrogen envelope) it slowed, giving rise to a re-verse shock. The reverse shock then propagated inwardsin mass towards the center of the star, and left a re-versed pressure gradient in its wake. This triggered theRayleigh-Taylor instability first between the helium shelland the hydrogen envelope, and then between the heliumand oxygen layers of the supernova. For model z15, thisinstability first arose between the oxygen layer and thehelium-hydrogen envelope. The instability grew, mixingtogether layers of increasing atomic mass as the reverseshock propagated inward through the star. Eventuallythe reverse shock passed by, and mixing stopped once theRayleigh-Taylor fingers had dissipated their momenta.In each simulation, the domain was enlarged five times,bringing the final effective resolution of each simulationto 16 , n , where n is the dimensionality of the sim-ulation. For model u15, four levels of refinement wereretained at the final stage of the simulation; for modelss15 and z15, two levels of refinement were retained.Mixing had effectively ceased by 4 . x seconds formodel s15, 3 . x seconds for model u15, and 1 . x seconds for model z15. Mixing had stopped by similartimes in 2D and 3D simulations. The models were run totwice this time to make sure that no additional mixingwould occur.Shown in Figure 2 are snapshots at identical timesfor the two- and three-dimensional simulations of models(from top) s15, u15, and z15. The abundance of He, O,Si, and Ni are shown in a logarithmic scale, startingclockwise from top left. These snapshots correspond totimes when mixing had effectively stopped in the differ-ent models. On the left are snapshots of the 2D models;at right are snapshots showing slices in the YZ plane atthe X axis. Upon first inspection, there appears to bevery little difference between the two and three dimen-sional simulations of all three models. The width of themixed region appears the same between the 2D and 3Dmodels, although mixing in the 3D models appears tobe more complete in this region than in the 2D mod-els. The shape of the instability differs slightly betweenthe 2D and 3D simulations. The mushroom shape ofthe instability is more clearly defined in 2D than in 3D,especially for model u15, where the instability had lesstime to grow and thus retained its original shape moreclearly. In the 3D models, the Rayleigh-Taylor instabilityappears more elongated, in line with what was found inthe single mode study of Kane et al. (2000). The great-est difference between the 2D and 3D simulations arisesin model s15, shown at the bottom of Figure 2, wherethere is actually less mixing of Ni in 2D than in 3D.In all models the Rayleigh-Taylor fingers have grownto the point where they have begun to interact withone another. Models z15 and s15, in which the RT in-stability had the longest time to grow, show the great-est degree of interaction. This interaction transfers en-ergy and momentum in the transverse directions to theblast wave, leading to a transition to turbulence in 3D,and to a chaotic regime in 2D (Remington et al. 2006).Miles et al. (2005) state that this transition to turbulencebegins to happen when the Rayleigh-Taylor fingers areabout 5 to 6 times longer than their initial wavelength.This has clearly happened in all simulations. Models z15 and s15 appear strikingly similar to the fully turbulentsimulations presented in Miles et al. (2005).The differences between 2D and 3D can be comparedmore quantitatively by examining Figure 3. This figuredisplays the radial average of abundances for significantelements in both 2D and 3D simulations of the threemodels after all mixing has stopped. The times shown arethe same times as in Figure 2. Figure 3 shows clearly thatthe width of the mixing region in the 2D and 3D modelsis nearly the same. Ni is the exception, and is slightlymore mixed in 2D than in 3D in the simulations of models15. Previous simulations of the Rayleigh-Taylor insta-bility in which these fingers do not interact show ≈ . x seconds, 3 . x sec-onds, and 1 . x seconds for models s15, u15, and z15,respectively, and so direct comparison of one model withanother would be misleading, as the material will slowdown as time goes on. It is clear from Figure 4 thatthe ultimate velocities after mixing has stopped and theRayleigh-Taylor instability has frozen out are essentiallythe same in the 2D and 3D simulations. This is differentfrom what Hammer et al. (2009) found in their simula-tions with larger degrees of initial asymmetry. They sawa high velocity tail for Ni, O, and Si. We do not see ahigh velocity tail here.The 3D renderings shown in Figure 5 give a clearerpicture of the shape of the Rayleigh-Taylor instabilitiesin our simulations than the slice through the 3D planeshown in Figure 2. We show renderings of the outermostsurfaces of carbon, oxygen, silicon, and nickel. Model s15is shown at the left, model u15 in the center, and modelz15 at the right. These renderings were taken at differenttimes, but after mixing had frozen out. Only in modelz15 have bits of the oxygen shell broken off and begun to
Fig. 2.—
Snapshots of elemental abundance in the 2D (left) and 3D (right) versions of (from top to bottom) models s15, u15, and z15.The 3D snapshots are slices at the X axis through the YZ plane. Because of artifacts in CASTRO, the extent of mixing is likely to beslightly exaggerated in slices at an axis. Nevertheless, the width of the mixed region in 2D and 3D appear very similar in all models. Thereis even slightly less mixing of Ni in model s15 in 3D than in 2D. The shape of the instability differs slightly between 2d and 3d in theexpected manner–the 2D shapes are slightly more “mushroomed” and the 3D shapes are more elongated, with the 3D models appearingto have transitioned to a fully turbulent regime, but otherwise the two simulations appear very similar.
Fig. 3.—
Abundance by mass of individual elements as a function of mass coordinate. The width of the mixed region is essentially thesame between the 2D and 3D cases, though 3D is smoother. The higher mass coordinate layers of the simulations have been omitted fromthe plot to better show the details of their inner portions.
Fig. 4.—
Mass of individual elements as a function of velocity per 100 km s − velocity bin. The amount of heavier elements such as Oand Ni at high velocities is essentially the same in 2D and 3D. Where the two differ, as they do most significantly for Ni and Si in modelss15 and z15, the 2D simulations actually show material mixed out to higher velocities than the 3D simulations. move through the H-He shell that surrounds the heavyelement core of this star. For no models have clumpsof heavier elements penetrated shells of lighter elements,as happened in the asymmetric explosion modeled byHammer et al. (2009). Even the broken-off clumps ofoxygen visible in the right hand panel of Figure 5 areencased in a layer of carbon. It is also apparent by ex-amining this figure that mixing has progressed further inthe simulations where the progenitor star was a red giant,model s15 and z15, than for model u15, whose progeni-tor died as a more compact blue star. The original shapeof the instability can be seen in the carbon surface (topmiddle panel). In addition, the Ni at the center of themodel is less deformed in this model than in previousmodels. The most vigourous mixing has taken place inmodel z15, which lacked a dense helium shell and diedas a red giant. DISCUSSION
Previous simulations of the Rayleigh-Taylor instabilityin core collapse supernovae have found significant dif-ferences between similar calculations performed in twoand three dimensions. Kane et al. (2000) examined thegrowth of Rayleigh-Taylor fingers arising from a singlemode perturbation, and found that the instability grew30 −
35% faster in 3D than in 2D in a 1987A-type pro-genitor model. This faster growth rate was not sufficientto replicate observations of 1987A, however. In the morerecent paper of Hammer et al. (2009), the authors per-formed simulations of a 1987A-type model in two and three dimensions, starting with a more realistic initialmodel from Scheck (2007). Their initial model was in thefirst stages of the explosion, and showed large-scale asym-metry in the heavy element core. They found, again, thatsingle perturbations grew around 30% faster in 3D thanin 2D. The authors proposed that artificial drag forcesarising from 2D geometry slow down the growth of theinstability in 2D simulations relative to 3D simulations,where these artificial drag forces are not present. Thiswas also described in Kane et al. (2000), who explainthat the rising bubbles have less effective matter to pushout of the way in 3D than in 2D.The Hammer et al. (2009) study also found that in 2D,rising bubbles of heavy elements became entrained in thedense He shell left behind by the shock. In 3D, however,these bubbles were able to burst through the He shell toreach the H envelope, replicating one of the key observa-tions of 1987A.Studies of the Rayleigh-Taylor instability outside thecontext of supernova explosions have also found that asingle Rayleigh-Taylor instability grows about 30% fasterin 3D than in 2D (Anuchina et al. 2004). Previous simu-lations of the Rayleigh-Taylor instability in a laboratorycontext (Miles et al. 2005) have found that the width ofthe mixed region is very similar between 3D and 2D,provided that the individual Rayleigh-Taylor fingers in-teract. Initially, the Rayleigh-Taylor fingers grow fasterin 3D than in 2D. Once the fingers begin to interact,3D simulations transition to turbulence more completelythan 2D simulations. 3D simulations can transfer mo-
Fig. 5.—
3D renderings of the outermost surface of carbon, oxygen, silicon, and nickel after mixing has frozen out. At left is model s15;center is u15, and right is z15. Note that while the scale and time are consistent within models, the 3 models are shown at different scalesand time: 5 x cm and 4 . x seconds for s15, 5 x cm and 3 . x seconds for model u15, and 3 x cm and 1 . x seconds formodel z15. For no model did heavier, inner layers penetrate the lighter, outer layers, though for model z15 clumps of material have brokenoff and moved a small distance into the H-He envelope. λ and amplitude η , such that η ≪ λ , thegrowth rate of the instability is given by Rayleigh (1883): dηdt = (cid:18) πAgλ (cid:19) / t (7)where g is a constant acceleration and t is time. A , the Atwood number, is defined by A ≡ ( ρ − ρ )( ρ + ρ ) (8)where ρ and ρ are the densities of the light and heavyfluids, respectively.As the densities of the two fluids become compera-ble, the Atwood number becomes smaller and the growthrate decreases. The growth rate of the bubbles in 3D isthus reduced, neatly canceling out the reduced drag thesesame bubbles experience in 3D as opposed to 2D. Thewidth of the mixed region remains essentially the samebetween 2D and 3D for cases where the Rayleigh-Taylorfingers interact with one another.It could be supposed that the reason the final state ofthe supernova remnant appears the same in our 2D and3D models is that the perturbations arising from the gridwere somehow larger in 2D than in 3D, and this effectexactly cancelled out the faster growth rates we expectto see in 3D. But this is not the case. Perturbationsarising from the grid are larger in 3D than in 2D. Theseperturbations arise because we try to fabricate a spheri-cal object with boxes, as required by the RZ and Carte-sian geometries used in our simulations. The maximumamount one would expect to depart from a sphere (or cir-cle) in space is equivalent to the longest distance across asingle cell, which in 2D is √ x and in 3D is √ x . Tomake certain, we performed the same set of calculationsin 2D at twice the resolution. These simulations give thesame results as our less refined simulations, indicatingthat our simulations are numerically converged for theproblem presented here.As in the simulations of Miles et al. (2005), Kane et al.(2000), and Hammer et al. (2009), we find that theRayleigh-Taylor instability initially grows faster in 3Dthan in 2D. As the Rayleigh-Taylor fingers grow to thepoint where they begin to interact with one another thesimulations begin to transition to a turbulent (3D) orchaotic (2D) regime. In 3D, momentum is transferredfrom the fingers in two dimensions, as opposed to onlyone dimension in 2D. The tranisiton to turbulence in 3Dresults in more complete mixing within the mixing regionthan the 2D simulation, reducing the Atwood number in3D relative to 2D. This has the effect of reducing thegrowth rate of the instability in 3D relative to 2D. Thereduced growth rate in 3D relative to 2D at late timescompensates for the fact that the bubbles do indeed ex-perience less drag in 3D than in 2D, and grow faster atearly times. By the time mixing has ceased, the width of the mixed region is nearly identical in 2D and 3D inmodels where the Rayleigh-Taylor fingers have time tointeract with one another.It can be seen from Eq. 7 that instabilites seeded bylong wavelength perturbations grow more slowly thanthose seeded by short wavelength perturbations. Ob-viously, the interior of a supernova deviates from theideal case: the acceleration is not constant, the fluidsare compressible, and η soon becomes of the same or-der as λ . In order for long wavelength perturbationsto dominate before mixing ceases there must be signif-icant initial power at long wavelengths. Because thereis less time for the instability to grow in blue stars, weexpect to see more traces of large perturbations at fairlylow wavenumbers in these stars than in red stars. Thesimulations of Hammer et al. (2009) have a significantamount of power at long wavelengths, so these pertur-bations grow fast enough to dominate before mixing istruncated.The initial perturbations in our simulations are a greatdeal smaller than those examined in previous studies.The spectrum of these perturbations is also flatter andmade up of higher modes. The perturbations are morerandom, though, it should be noted, not entirely so, andare largest around the 30 and 60 ◦ marks on the grid. Thiscan be seen during the initial stages of the developmentof the instability, though not at later times, for models15. The enhanced perturbations at these points is likelybehind the placement of the clumps visible in the 3Drendering of model z15 in Figure 5.Mixing had more time to develop in our two red gi-ant stars, models z15 and s15, than it generally doesin the blue giant 1987A-type progenitor models studiedin Kane et al. (2000) and Hammer et al. (2009). Mix-ing froze out after 4 x and 1 . x seconds for thered models s15 and z15, respectively, over an order ofmagnitude later than in the simulation of Hammer et al.(2009), in which mixing froze out after 9 x seconds.In our blue u15 model, mixing froze out after 4 x sec-onds, which is still nearly 5 times longer than in theHammer et al. (2009) model. The individual Rayleigh-Taylor instabilities in all simulations grew for sufficienttimes that they interacted, though more interaction tookplace in the red models. This was not the case in thesingle-mode perturbation explored in Kane et al. (2000),or the few large perturbation modes in (Hammer et al.2009). Smaller wavelengths grow faster than longerwavelengths, so longer wavelength perturbations are lesslikely to become nonlinear. In these models, as wellas the single instability investigated in Anuchina et al.(2004), the Rayleigh-Taylor fingers did not interact withone another, ensuring that the relative drag on the fingerswas the only significant difference between 2D and 3D.This accounts for the around 30% faster growth ratesdiscovered in these simulations. The simplicity of theinitial perturbation spectrum can delay the transition toturbulence, and hence make mixing appear more efficientin 3D than in 2D (Miles et al. 2005).For long enough wavelength perturbations, sphericalgeometry will ensure that the perturbations diverge, sothey never interact and thus never enter a turbulentregeime. This will quickly lead to significant depar-tures from axisymmetry. There may exist a transitionmode perturbation which would be the longest wave-1length mode that could lead to turbulent interaction.The instabilities will diverge without interaction if thewavelength of an instability is greater than the height atwhich it would become nonlinear. This tranistion modeought to be dependent on presupernova structure: redsupergiants have longer time for the Rayleigh-Taylor in-stability to grow, and thus potentially a lower mode thanblue giants.The Rayleigh-Taylor instability may “forget” its initialconditions at late times. This can happen though an in-verse cascade, where shorter wavelength bubbles mergeto form longer wavelength structures. Under idealizedconditions, the flow may become self similar. This tran-sition to self-similarity likely happens during the earlystages of Rayleigh-Taylor growth in our simulations, ac-counting for the similar end appearance of our simula-tions despite increasing the resolution (and decreasingthe initial perturbations). This inverse cascade can onlylead to the instabilities becoming so big. A necessarycondion for self-similarity is that λ ≪ h ≪ L (9)where h is the height of the instability and L is thesize of the region in which the instability is growing.The merger of short wavelength instabilities to largerwavelength instabilites should be truncated long before λ becomes comperable to L , which would produce large-scale asymmetries. Large-scale departures from symme-try must arise from the supernova explosion mechanismitself, as they cannot arise through an inverse cascade ofthe Rayleigh-Taylor instability.Recent observations suggest that asymmetry incore-collapse type supernova explosions is common(Wang & Wheeler 2008). For SNe type II, this asym-metry is characterized by a dominant polarization an-gle, which is most likely due to a directed, jet-likeflow. Polarization data show that most SNe also exhibitcomposition-dependent, large scale departures from ax-isymmetry. The degree of departure from axisymmetrymay be less in SNeII-P; Chugai et al. (2005) and Chugai(2006) are able to match observations of one SNeII-P, SN2004dj, by using a bipolar distribution of Ni.This fits with the emerging picture of collapse supernovamechanisms, which are likely to be inherently mulAtidi-mensional and aspherical (Janka et al. 2007) The SASIexplosion mechanism (Burrows et al. 2006; Scheck et al.2006; Marek & Janka 2009) seems to be one of the bettercandidates for reviving the stalled shock and successfullyexploding a massive star, and it is inherently asymmet-ric. How relevant, then, are spherical explosion models?This is the first time the post-explosion hydrodynamicsof a wide range of supernova progenitors has been stud-ied in 3D. It is reasonable to model Pop III and Pop IIcore collapse SNe in 3D first as spherical explosions asno observations of these primordial supernovae have yetbeen made. This study will contextualize future studiesof supernovae progenitors with asymmetric explosions.For type II SNe in which departures from asymmetryare small and dominated by an axisymmetric component,2D simulations may do as well as 3D simulations in repro-ducing the axisymmetric features of the supernova rem-nant. In these stars, Rayleigh-Taylor fingers of all butthe lowest mode order are expected to grow far enough that they begin to interact with one another. Simulationsin 2D of such stars can be used to predict nucleosynthe-sis and fallback for such explosions with a fair degree ofaccuracy, as was done in Joggerst et al. (2010). 2D calcu-lations may also be used to more accurately predict lightcurves of SNe for which departures from axisymmetry donot dominate, which would use far less computational re-sources than similar calculations performed in 3D. CONCLUSIONS
We have performed two and three dimensional simu-lations of supernova models with three different metal-licities: model z15, with Z = 0 Z ⊙ , model u15, withZ= 10 − Z ⊙ , and model s15, with Z= Z ⊙ . The initialperturbations are seeded by the grid. Calculations per-formed at twice the resolution in 2D indicate that oursolutions are numerically converged. We find no signifi-cant difference in the width of the mixed region betweenour 3D and 2D simulations of the same supernova. Pre-vious studies of Rayleigh-Taylor mixing in 1987A-typesupernova models found that mixing was more efficientand individual Rayleigh-Taylor instabilities grew about30% faster in 3D than in 2D. Kane et al. (2000) stud-ied single-mode density perturbations with amplitudes of10% of mode l = 10, and found that perturbations grew30-35% faster in 3D than in 2D. The more recent studyof Hammer et al. (2009) examined the post-explosion hy-drodynamics of a 1987A-type model with a fully 3D ex-plosion model from Scheck (2007) with large scale asym-metries at early times. They found that the instabilitygrew about 30% faster in 3D than in 2D.While we also find a faster initial growth rate in 3Dthan in 2D in our simulations, the growth rate in 3D de-clines once the simulations transition to turbulence, andthe final width of the mixed region is the same in 2Dand 3D. Because the Rayleigh-Taylor instability in ourmodel is seeded by perturbations arising from the grid,our perturbations have a lower amplitude, and a flatterspectrum shifted towards higher modes than either theKane et al. (2000) or Hammer et al. (2009) studies. Inour progenitor models the Rayleigh-Taylor instability hasmore time to develop and grow than in previous simula-tions. With a growth time on the order of 10 seconds,instabilities in models s15 and z15 grow for more than10 times longer than the simulations of Hammer et al.(2009), in which mixing froze out at around 9 x sec-onds. In model u15, which died as a blue supergiant,mixing ceases after 4 x seconds, which is a shorter timethan the red giant models in this survey but is still nearly5 times longer than the the simulations of Hammer et al.(2009).In our progenitor models and with perturbationsseeded by the grid the RT instability has enough timeto grow that the fingers of the instability begin to inter-act with one another. This interaction transfers momen-tum and energy transverse to the shock, allowing the 3Dsimulations to mix more fully than 2D simulations. Themore efficient transfer of momentum and energy in 3Drenders the 3D simulations fully turbulent in the mixedregion. This higher degree of mixing in 3D than in 2Dreduces the Atwood number in the 3D simulations rel-ative to the Atwood number in the 2D simulations. Alower Atwood number leads to a reduction in the bub-ble growth rate in 3D, which compensates for the lower2effective drag in 3D. The result is that the width of themixed region is virtually identical between simulationsperformed in 2D and those performed in 3D. 3D simula-tions have transitioned to full turbulence and are morecompletely mixed than 2D simulations. This was foundin a laboratory context in the simulations of Miles et al.(2005). For simulations in which interactions betweenthe Rayleigh-Taylor fingers do not occur, either becauseonly one finger is modeled (Anuchina et al. 2004) or theRayleigh-Taylor fingers grow for a short enough time thatthey never interact (Kane et al. 2000), the greater effec-tive drag experienced in 2D versus 3D will continue todominate the simulation, ensuring that the width of themixed region will be larger in 3D than in 2D. For simula-tions in which the Rayleigh-Taylor instability has enoughtime to grow that the fingers begin to interact with oneanother, the width of the mixed region will be the samein 2D and 3D simulations.Instabilities arising from long wavelength perturba-tions take longer to grow than those arising from shortwavelength perturbations. Large amplitude perturba-tions of long wavelength are necessary for such instabili-ties to grow substantially. These large-scale fingers maydiverge before they can transition to turbulence, and thusmay give rise to departures from axisymmetry. The largeamplitude, long wavelength perturbations present in themore realistic explosion model of Hammer et al. (2009)ensure that such fingers don’t interact, which results inthe significant differences between their 2D and 3D sim-ulations. Rayleigh-Taylor theory predicts that an inversecascade, in which shorter wavelength instabilities mergeto form longer wavelength instabilites, should be trun-cated long before the wavelegth λ becomes comparableto the size of the region in which the instability is grow-ing. Large scale departures from symmetry do not arisethrough an inverse cascade in our 3D simulations; they must have their origins in the explosion itself.It is likely that many, if not most, core-collapse super-novae explosions are highly asymmetric. Observationsnot just of 1987A but of other supernova, as well the“kicks” observed in pulsars, indicate that asymmetry isa near-universal phenomenon amongst the explosions ofmassive stars. From the theoretical side, it seems in-creasingly likely that viable supernova explosions, likethe SASI mechanism, are inherently asymmetric. In su-pernovae that die as red giants the Rayleigh-Taylor in-stability will have enough time to develop that individualfingers will interact with one another. If the asymmetriesin such explosions are dominated by axisymmetric com-ponents, 2D models might well be good predictors of thewidth of the mixed region and the axisymmetric featuresof the remnant. Future work should investigate asym-metric explosions in a wide variety of progenitor models.Work at UCSC and LBL was supported in part by theSciDAC Program under contract DE-FC02-06ER41438.Work at LANL was carried out under the auspices ofthe National Nuclear Security Administration of the U.S.Department of Energy at Los Alamos National Labora-tory under Contract No. DE-AC52-06NA25396. StanWoosley was supported by the NSF and NASA as wellas SciDAC National Science Foundation (AST 0909129)and the NASA Theory Program (NNX09AK36G). Thesimulations were performed on the open cluster Coyoteat Los Alamos National Laboratory. Additional comput-ing resources were provided on the Pleiades computerat UCSC under NSF Major Research Instrumentationaward number AST-0521566. C.C.J. would like to thankAndrew Aspden for assistance with parallelzing a datareduction routine. Figure 5 was created using VisIt.becomes comparableto the size of the region in which the instability is grow-ing. Large scale departures from symmetry do not arisethrough an inverse cascade in our 3D simulations; they must have their origins in the explosion itself.It is likely that many, if not most, core-collapse super-novae explosions are highly asymmetric. Observationsnot just of 1987A but of other supernova, as well the“kicks” observed in pulsars, indicate that asymmetry isa near-universal phenomenon amongst the explosions ofmassive stars. From the theoretical side, it seems in-creasingly likely that viable supernova explosions, likethe SASI mechanism, are inherently asymmetric. In su-pernovae that die as red giants the Rayleigh-Taylor in-stability will have enough time to develop that individualfingers will interact with one another. If the asymmetriesin such explosions are dominated by axisymmetric com-ponents, 2D models might well be good predictors of thewidth of the mixed region and the axisymmetric featuresof the remnant. Future work should investigate asym-metric explosions in a wide variety of progenitor models.Work at UCSC and LBL was supported in part by theSciDAC Program under contract DE-FC02-06ER41438.Work at LANL was carried out under the auspices ofthe National Nuclear Security Administration of the U.S.Department of Energy at Los Alamos National Labora-tory under Contract No. DE-AC52-06NA25396. StanWoosley was supported by the NSF and NASA as wellas SciDAC National Science Foundation (AST 0909129)and the NASA Theory Program (NNX09AK36G). Thesimulations were performed on the open cluster Coyoteat Los Alamos National Laboratory. Additional comput-ing resources were provided on the Pleiades computerat UCSC under NSF Major Research Instrumentationaward number AST-0521566. C.C.J. would like to thankAndrew Aspden for assistance with parallelzing a datareduction routine. Figure 5 was created using VisIt.