Three-dimensional simulations of rapidly rotating core-collapse supernovae: finding a neutrino-powered explosion aided by non-axisymmetric flows
We report results from a series of three-dimensional (3D) rotational core-collapse simulations for and 27 stars employing neutrino transport scheme by the isotropic diffusion source approximation. By changing the initial strength of rotation systematically, we find a rotation-assisted explosion for the 27 progenitor , which fails in the absence of rotation. The unique feature was not captured in previous two-dimensional (2D) self-consistent rotating models because the growing non-axisymmetric instabilities play a key role. In the rapidly rotating case, strong spiral flows generated by the so-called low instability enhance the energy transport from the proto-neutron star (PNS) to the gain region, which makes the shock expansion more energetic. The explosion occurs more strongly in the direction perpendicular to the rotational axis, which is different from previous 2D predictions.
keywords:stars: interiors – stars: massive – supernovae: general.
After a half-century of extensive research, theory and simulations are now converging to a point that multi-dimensional (multi-D) hydrodynamic instabilities play a crucial role in the neutrino mechanism of core-collapse supernovae (CCSNe, see Foglizzo et al. (2015); Mezzacappa et al. (2015); Janka (2012); Burrows (2013); Kotake et al. (2012) for reviews). Multi-D fluid motions associated with convective overturn and the standing-accretion-shock-instability (SASI, Blondin et al. (2003)) are the keys because buoyant and turbulent flows increase the neutrino heating efficiency in the gain region, triggering the runaway expansion of the shock. In fact, a growing number of neutrino-driven models have been recently reported in self-consistent two-dimensional (2D) simulations, which has strengthened our confidence in the multi-D neutrino-driven paradigm (e.g., Marek & Janka (2009); Suwa et al. (2010); Müller et al. (2012); Dolence et al. (2015); Pan et al. (2015); Nakamura et al. (2015); O’Connor & Couch (2015)).
This success, however, is now raising new questions. With the exception in Bruenn et al. (2013), the explosion energies so far reported in these 2D models are generally not sufficient to explain the canonical supernova kinetic energy of erg. Moreover, the most challenging self-consistent three-dimensional (3D) simulations with spectral neutrino transport have failed to produce explosions for , , and progenitors (Hanke et al., 2013; Tamborra et al., 2014), or, in a few successful cases, showed much delayed neutrino-driven shock revival in 3D than in 2D (e.g., Lentz et al. (2015) and Melson et al. (2015)), leading to even smaller explosion energies in 3D compared to 2D (Takiwaki et al. (2014)).
One of the prime candidates to predominantly affect the CCSN explodability is general relativity (GR, e.g., Müller et al. (2012); Kuroda et al. (2012)). Rotation (e.g., Marek & Janka (2009); Suwa et al. (2010)), magnetic fields (Endeve et al., 2012; Obergaulinger et al., 2014), and inhomogeneities in the progenitor core (Couch & Ott, 2015) are also attracting much attention to turn an unsuccessful multi-D model into a successful explosion.
In this Letter, we focus on the roles of rotation and report results from a series of 3D rotational core-collapse simulations with spectral neutrino transport for and stars. We briefly describe our numerical approach in Section 2. We discuss our results in Section 3, followed by a summary in Section 4
2 Numerical Setup and Progenitor model
Initial conditions are taken from the and pre-supernova progenitors of Woosley et al. (2002). The models, which have been used in Takiwaki et al. (2012, 2014); Hanke et al. (2013); Müller (2015), are useful to clearly explore the impacts of rotation The initially constant angular frequency of or rad/s is imposed inside the iron core with a cut-off () outside. Although these angular frequencies are close to the high-end of those from most recent stellar evolution models (e.g., Heger et al. (2000, 2005), see also discussions in Ott et al. (2006)), we assume such rapid rotation to clearly see the impacts of rotation in this study. The model name is labeled as ”s11.2-R1.0-3D”, which represents the 11.2 model with rad/s that is computed in 3D simulation.
Our numerical method is based on that in Takiwaki et al. (2014) except several points. We use the equation of state (EOS) by Lattimer & Douglas Swesty (1991) (incompressibility MeV). Our code employs a high-resolution shock capturing scheme with an approximate Riemann solver of Einfeldt (1988) (see Nakamura et al. (2015) for more details). For the calculation presented here, self-gravity is computed by a Newtonian monopole approximation111Our 3D rotating models with an improved multipole approximation of gravity (e.g., Couch et al. (2013)) explode more energetically than those only with the monopole contribution (see, more details in Takiwaki et al. in preparation).. Our fiducial 3D models are computed on a spherical polar grid with a resolution of = , in which non-equally spatial radial zones covers from the center to an outer boundary of 5000 km.222This choice of the outer boundary position was shown to be insignificant especially in the simulation timescale ( 300 ms postbounce) in this work (see section 2.3 in Nakamura et al. (2015)). Our spatial grid has a finest mesh spacing km at the center and is better than 2% at km. For a numerical resolution test, we compute high-resolution runs with = .
In total, we have computed nine 3D models, which consists of six models with the fiducial resolution (i.e., the two progenitors with rad/s) and three high-resolution runs for the 11.2 model. By using the fastest computer in Japan, it typically took 2 months (equivalently 15 Pflops-day computational resources) for each of the high-resolution runs.
Figure 1 summarizes the blast morphology for the (left panels) and star (right panels), which are helpful to compare the hydrodynamics features between the non-rotating (top) and rapidly rotating (bottom) models, respectively.
In the non-rotating models, s11.2-R0.0-3D (top left) shows typical features of neutrino-driven convection in the postshock regions. The rising plumes grow stronger and larger in angular size from the initial small mushroom-like Rayleigh-Taylor fingers.
In models with rapid rotation, a clear oblate explosion is obtained for model s27.0-R2.0-3D (bottom right), in which the revived shock expands more strongly in the equatorial plane. This feature is only weakly visible for model s11.2-R2.0-3D (bottom left) due to the early shock revival (see also, top panel of Figure 2). Later we present detailed analysis of the origin of the oblate explosion and point out a new aspect of rapid rotation for assisting explosions.
Before going into detail, let us shortly summarize the evolution of the shock and (diagnostic) explosion energy of all the computed models in Figure 2. The top panels are for the 11.2 series with different and different numerical resolution (with the high resolution being ended with H). The average shock radii of the standard resolution models (solid line) and high resolution models (dashed line) do not deviate from each other. It is important to present that our results do not strongly depend on the grid size333Apparently our resolution is not sufficient for reproducing realistic viscosity (Couch & Ott, 2015)). The convergence may be partly due to the diffusive feature of the HLLE scheme employed in this work (e.g., Radice et al. (2015))..
The bottom panel of Figure 2 shows that all the variations of the non-rotating 27 progenitor star do not trend toward an explosion very clearly during the simulation, whereas the rapidly rotating model does so (red solid line) with the diagnostic energy much bigger than those of the non-rotating models. Using the same progenitor (27 ), the hydrodynamics features in the non-rotating models are qualitatively consistent with Hanke et al. (2013).
From here we proceed to focus on the impacts of rotation on the shock dynamics. As shown in the top panel of Figure 2, the shock radius of the rotating models (red and green) does not deviate from that of the non-rotating model (blue). Regarding the diagnostic energies, the rotational model shows less energetic explosion. That is because the neutrino luminosities of the rotating models are slightly smaller than that of the non-rotating models (e.g., Marek & Janka (2009)). In the case of such an early shock revival observed in the light progenitor, the rotation does not help the onset of the explosion and even decrease the explosion energy. It should be noted that time derivative of the diagnostic explosion energies, , are also small in these models and they seem to be almost saturated in spite of the relatively short postbounce time ( ms).
For the heavier progenitor of 27.0, the situation is inverted. As shown in the bottom panel of Figure 2, only the rapidly rotating model explodes. From here, we focus on the 27.0 models since the effects of the rotation is most distinct. We try to provide a new interpretation of the rotating explosion. For later convenience, let us define the deviation from the spherical (angle-)averaged variable as , where represents the angle average of . In the panels of Figure 3, the deviation of the density, pressure, and radial velocity in the equatorial plane of R2.0-3D is shown, respectively. The spiral flows (colored by red in the panels) are associated with high density and pressure, which pushes the matter outward and assists the explosion.
At 80 ms after bounce, the value of the rotation to the gravitational energy () in the vicinity of the PNS exceeds 6% and the iso-density surface of the PNS begins to be deformed with the dominance of mode. This behavior is quite similar to Ott et al. (2005) who were the first to observe the growth of the low- instability in the PNS context. Figure 3 shows the amplitude of the dipole mass deformation in the equatorial plane, whose definition is , where represents the spherical harmonics. After ms postbounce, the deformation amplitude in the center ( at 20 km, top panel of Figure 3) approaches a percent level (as is consistent with Ott et al. (2005)), simultaneously, the spiral flows begin to extend outwards later on (e.g., left and center panels of Figure 3). The deformation amplitude near at the stalled shock ( at 130 km) peaks at around ms postbounce with the maximum amplitude of in density. It should be noted that the rotation energy of the mildly rotating model does not reach the threshold of the onset of low- instability. In this case the hydrodynamic evolution is similar to that in the non-rotating model.
The spiral flow is not merely “pattern” of the density and the pressure. This wave really brings the mass and the thermal energy from the central to the outer region. The right panel of Figure 3 shows the snapshot of the radial velocity. In red region, matter has a positive radial velocity and in the blue region vice versa. From these three panels of Figure 3, one can see that the matter with the high density and high thermal energy (pressure) goes outward and that with low density and low thermal energy goes inward.
To quantify the energy transport by the spiral flows, we utilize the concept of the (turbulent) energy flux as
where represents the density, pressure, and the radial velocity respectively and the dashed values denote the deviation from the average. This flux means energy transport caused by instabilities in the context of Reynolds decomposition (Murphy & Meakin, 2011; Murphy et al., 2013). The kinetic energy can be ignored since that is typically ten times smaller than the internal energy. The total energy transport is written as the sum of one-dimensional energy flux, , and the modification by the spiral flows in this case, , since the background flows (without the spiral flows) yield to the following relations . Note that we actually estimate the flux as because this treatment empirically suppresses artificial overestimation due to the steep density gradient near at the shock (e.g., Murphy et al. (2013)).
Figure 3 shows the energy flux and radial profile of the entropy. The blue and green line corresponds to the non-rotating model and mildly rotating model, respectively. As was previously known in the non-rotating model, the turbulent flux (blue line) makes a significant contribution to the energy transport only in the postshock region with negative entropy gradient (compare two panels in Figure 3). On the other hand, it is shown for the rapidly rotating model (red line) that the energy flux continuously contributes to the energy transport from the central region to behind the shock.
Although similar results are obtained by Ott et al. (2012), the detailed mechanism was not discussed since the work was dedicated to the gravitational wave emission. In addition, there are severe limitations in Ott et al. (2012). Since an octant symmetry was assumed in their work, mode which is dominant in our model, was not resolved. Their neutrino transport is based on the leakage scheme and the rotation aided explosion was obtained only in adiabatic models whose neutrino cooling and heating is switched off (see their Figure 7).
How much does that spiral activity power the shock ? To evaluate this, we estimate the rate of the energy change in the gain region () as
where is the gravitational potential. Note that can be clearly distinguished from the rate of the energy increased by the neutrino heating at that moment (see eq. (107) of Janka (2001) for more general description). As shown in the top panel of Figure 4, for the rapidly rotating model sharply grows from zero at ms postbounce and peaks around ms postbounce. The timescale coincides with the epoch when the developing spiral flow from the center (Figure 3) thrusts into the gain region. Note also that both of the time variability of (Figure 4) and at the radius of 130 km (bottom panel of Figure 3) is correlated. In the non-rotating or mildly rotating cases (blue and green lines), on the other hand, takes negative value after 130 ms postbounce. In fact, the shock revival has not yet been obtained for these models in the simulation timescale. These analyses naturally support that the rapidly rotating model trends toward an explosion because of the energy gain sustained by the spiral activity.
The energy gain due to the spiral activity is comparable to that of the neutrino heating. The bottom panel of Figure 4 shows the time evolution of that is an integrated neutrino heating rate in the gain region. The value is erg/s in all the models and the contribution of in the rapidly rotating model is the same order to that of the net heating rate. The rapid rotating model is surely energized by the spiral flows.
4 Summary and Discussion
We reported results from a series of 3D core-collapse simulations for and 27 stars using the IDSA scheme for spectral neutrino transport. By changing the initial strength of rotation systematically, we observed a new indication of rotation-assisted explosion for the model. This model fails to explode in the corresponding non- or mildly rotating models. In the rapidly rotating case, strong spiral flows generated by the low instability enhance the energy transport from the PNS to the gain region, which makes the shock expansion more energetic. These impacts of the rotation were also shown to be sensitive to the progenitor models. For the lighter progenitor, the early shock revival is little affected by the rotation. In this case, the explosion energy of the rotating model becomes weaker, as previously pointed out by Marek & Janka (2009), predominantly due to the smaller neutrino luminosity.
The major limitation of this study would be omission of the magnetic fields. At present, saturation level of the field amplification due to magnetorotational instability (MRI) is still under a hot debate (Rembiasz et al., 2016). Depending on the MRI’s growth rate, the magnetic fields could be amplified during our simulation time, possibly fostering the onset of explosion (Mösta et al., 2015; Sawai et al., 2013). To test this, a high-resolution 3D MHD model is needed, which is another major undertaking.
The final (averaged) angular velocity of model s27-R2.0 is in the vicinity of the PNS, which is apparently too fast to be reconciled with observations of canonical radio pulsars (e.g., Faucher-Giguère & Kaspi (2006), see also Ott et al. (2006)). Such rapid rotation (e.g., rad/s in this work), rare as it should be (see discussions in Woosley & Bloom (2006); Cantiello et al. (2007); Fuller et al. (2015); Chatzopoulos et al. (2016)), is attracting much attention as to their possible relevance to hyper-energetic explosions (e.g., Iwamoto et al. (1998); Mazzali et al. (2008)). These events are also hypothetically related to the formation of magnetars and collapsars (e.g., Mészáros (2006) for a review). We hope that this work could give a momentum for theorists to pay more attention to 3D models with rapid rotation (plus magnetic fields, e.g., Mösta et al. (2015); Masada et al. (2015) for recent discoveries), which could possibly illuminate the yet unexplored variety of the explosion dynamics where non-axisymmetric instabilities play a substantial role.
T.T. is grateful to K. Nitadori and J. Makino for the tuning of our code. The computations in this research were performed on the K computer of the RIKEN AICS through the HPCI System Research project (Project ID:hp120304) together with XC30 of CfCA in NAOJ. YS was supported by JSPS postdoctoral fellowships for research abroad. This study was supported in part by the Grants-in-Aid for the Scientific Research from the Ministry of Education, Science and Culture of Japan (Nos. 24103006, 24244036, 26104001, 26707013, 26870823, 15H01039 and 15H00789) and by HPCI Strategic Program of Japanese MEXT.
- Blondin et al. (2003) Blondin J. M., Mezzacappa A., DeMarino C., 2003, ApJ, 584, 971
- Bruenn et al. (2013) Bruenn S. W., et al., 2013, ApJ, 767, L6
- Burrows (2013) Burrows A., 2013, Reviews of Modern Physics, 85, 245
- Cantiello et al. (2007) Cantiello M., Yoon S.-C., Langer N., Livio M., 2007, A&A, 465, L29
- Chatzopoulos et al. (2016) Chatzopoulos E., Couch S. M., Arnett W. D., Timmes F. X., 2016, preprint, (arXiv:1601.05816)
- Couch & Ott (2015) Couch S. M., Ott C. D., 2015, ApJ, 799, 5
- Couch et al. (2013) Couch S. M., Graziani C., Flocke N., 2013, ApJ, 778, 181
- Dolence et al. (2015) Dolence J. C., Burrows A., Zhang W., 2015, ApJ, 800, 10
- Einfeldt (1988) Einfeldt B., 1988, SIAM Journal on Numerical Analysis, 25, 294
- Endeve et al. (2012) Endeve E., Cardall C. Y., Budiardja R. D., Beck S. W., Bejnood A., Toedte R. J., Mezzacappa A., Blondin J. M., 2012, ApJ, 751, 26
- Faucher-Giguère & Kaspi (2006) Faucher-Giguère C.-A., Kaspi V. M., 2006, ApJ, 643, 332
- Foglizzo et al. (2015) Foglizzo T., et al., 2015, Publ. Astron. Soc. Australia, 32, e009
- Fuller et al. (2015) Fuller J., Cantiello M., Lecoanet D., Quataert E., 2015, ApJ, 810, 101
- Hanke et al. (2013) Hanke F., Müller B., Wongwathanarat A., Marek A., Janka H.-T., 2013, ApJ, 770, 66
- Heger et al. (2000) Heger A., Langer N., Woosley S. E., 2000, ApJ, 528, 368
- Heger et al. (2005) Heger A., Woosley S. E., Spruit H. C., 2005, ApJ, 626, 350
- Iwamoto et al. (1998) Iwamoto K., et al., 1998, Nature, 395, 672
- Janka (2001) Janka H.-T., 2001, A&A, 368, 527
- Janka (2012) Janka H.-T., 2012, Annual Review of Nuclear and Particle Science, 62, 407
- Kotake et al. (2012) Kotake K., Sumiyoshi K., Yamada S., Takiwaki T., Kuroda T., Suwa Y., Nagakura H., 2012, Progress of Theoretical and Experimental Physics, 2012, 010000
- Kuroda et al. (2012) Kuroda T., Kotake K., Takiwaki T., 2012, ApJ, 755, 11
- Lattimer & Douglas Swesty (1991) Lattimer J. M., Douglas Swesty F., 1991, Nuclear Physics A, 535, 331
- Lentz et al. (2015) Lentz E. J., et al., 2015, ApJ, 807, L31
- Marek & Janka (2009) Marek A., Janka H.-T., 2009, ApJ, 694, 664
- Masada et al. (2015) Masada Y., Takiwaki T., Kotake K., 2015, ApJ, 798, L22
- Mazzali et al. (2008) Mazzali P. A., et al., 2008, Science, 321, 1185
- Melson et al. (2015) Melson T., Janka H.-T., Bollig R., Hanke F., Marek A., Müller B., 2015, ApJ, 808, L42
- Mészáros (2006) Mészáros P., 2006, Reports on Progress in Physics, 69, 2259
- Mezzacappa et al. (2015) Mezzacappa A., et al., 2015, preprint, (arXiv:1501.01688)
- Mösta et al. (2015) Mösta P., Ott C. D., Radice D., Roberts L. F., Schnetter E., Haas R., 2015, Nature, 528, 376
- Müller (2015) Müller B., 2015, MNRAS, 453, 287
- Müller et al. (2012) Müller B., Janka H.-T., Heger A., 2012, ApJ, 761, 72
- Murphy & Meakin (2011) Murphy J. W., Meakin C., 2011, ApJ, 742, 74
- Murphy et al. (2013) Murphy J. W., Dolence J. C., Burrows A., 2013, ApJ, 771, 52
- Nakamura et al. (2015) Nakamura K., Takiwaki T., Kuroda T., Kotake K., 2015, PASJ, 67, 107
- O’Connor & Couch (2015) O’Connor E., Couch S., 2015, preprint, (arXiv:1511.07443)
- Obergaulinger et al. (2014) Obergaulinger M., Janka H.-T., Aloy M. A., 2014, MNRAS, 445, 3169
- Ott et al. (2005) Ott C. D., Ou S., Tohline J. E., Burrows A., 2005, ApJ, 625, L119
- Ott et al. (2006) Ott C. D., Burrows A., Thompson T. A., Livne E., Walder R., 2006, Astrophys. J., Suppl. Ser., 164, 130
- Ott et al. (2012) Ott C. D., et al., 2012, Phys. Rev. D, 86, 024026
- Pan et al. (2015) Pan K.-C., Liebendörfer M., Hempel M., Thielemann F.-K., 2015, preprint, (arXiv:1505.02513)
- Radice et al. (2015) Radice D., Couch S. M., Ott C. D., 2015, preprint, (arXiv:1501.03169)
- Rembiasz et al. (2016) Rembiasz T., Obergaulinger M., Cerdá-Durán P., Müller E., Aloy M. A., 2016, MNRAS, 456, 3782
- Sawai et al. (2013) Sawai H., Yamada S., Suzuki H., 2013, ApJ, 770, L19
- Suwa et al. (2010) Suwa Y., Kotake K., Takiwaki T., Whitehouse S. C., Liebendörfer M., Sato K., 2010, PASJ, 62, L49
- Takiwaki et al. (2012) Takiwaki T., Kotake K., Suwa Y., 2012, ApJ, 749, 98
- Takiwaki et al. (2014) Takiwaki T., Kotake K., Suwa Y., 2014, ApJ, 786, 83
- Tamborra et al. (2014) Tamborra I., Raffelt G., Hanke F., Janka H.-T., Müller B., 2014, Phys. Rev. D, 90, 045032
- Woosley & Bloom (2006) Woosley S. E., Bloom J. S., 2006, ARA&A, 44, 507
- Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015