Supernova Anyflip

[Amit Bolds]

Multidimensional fluid flow plays a paramount role in the explosions of massive stars as core-collapse supernovae. In recent years, three-dimensional (3D) simulations of these phenomena have matured significantly. Considerable progress has been made towards identifying the ingredients for shock revival by the neutrino-driven mechanism, and successful explosions have already been obtained in a number of self-consistent 3D models. These advances also bring new challenges, however. Prompted by a need for increased physical realism and meaningful model validation, supernova theory is now moving towards a more integrated view that connects multi-dimensional phenomena in the late convective burning stages prior to collapse, the explosion engine, and mixing instabilities in the supernova envelope. Here we review our current understanding of multi-D fluid flow in core-collapse supernovae and their progenitors. We start by outlining specific challenges faced by hydrodynamic simulations of core-collapse supernovae and of the late convective burning stages. We then discuss recent advances and open questions in theory and simulations.

How precisely the envelope is ejected has remained one of the foremost questions in computational astrophysics ever since the first modeling attempts in the 1960s (Colgate et al. 1961; Colgate and White 1966). In this review we focus on the critical role of multi-dimensional (multi-D) fluid flow during the supernova explosion itself and the final pre-collapse stages of their progenitors.

For pedagogical reasons, it is preferable to commence our brief exposition of multi-D hydrodynamic effects with the supernova explosion mechanism rather than to follow the sequence of events in nature, or historical chronology.

In principle, the collapse of an iron core to a neutron star opens a reservoir of several \(10^53\, \mathrm erg\) of potential energy, which appears more than sufficient to account for the typical inferred kinetic energies of observed core-collapse supernovae of order \(10^51\, \mathrm erg\) (see, e.g., Kasen and Woosley 2009; Pejcha and Prieto 2015).

In reality, even non-rotating progenitors are not spherically symmetric at the onset of collapse. Outside the iron core, there are typically several active convective burning shells (Fig. 1a) that will collapse in the wake of the iron core within hundreds of milliseconds. It was realized in recent years that the infall of asymmetric shells can be important for shock revival (Couch and Ott 2013; Mller 2015; Couch et al. 2015; Mller et al. 2017a).

The multi-D structure of supernova progenitors is thus directly relevant for the neutrino-driven mechanism, but the potential ramifications of multi-D effects during the pre-collapse phase are in fact much broader: How do they affect mixing at convective boundaries, and hence the evolution of the shell structure on secular time scales? How do they affect the angular momentum distribution and magnetic fields in supernova progenitors?

Observations contain abundant clues about the multi-D nature of core-collapse supernova explosion. Large birth kicks of neutron stars (Hobbs et al. 2005; Faucher-Gigure and Kaspi 2006; Ng and Romani 2007) and even black holes (Repetto et al. 2012) cannot be explained by stellar dynamics alone and require asymmetries in the supernova engine. There is also evidence for mixing processes during the explosion and large-scale asymmetries in the ejecta from the spectra and polarization signatures of many observed transients (e.g., Wang and Wheeler 2008; Patat 2017), and from young supernova remnants like Cas A (Grefenstette et al. 2014).

In this review, we are primarily concerned with the numerical techniques for modeling multi-D fluid flow in core-collapse supernovae and their progenitors, and with our current understanding of the theory and phenomenology of the multi-D fluid instabilities. Although multi-D effects are relevant to virtually all aspects of core-collapse supernova theory, we can only afford cursory attention to many of them in order to keep this overview focused.

We shall start by discussing the governing equations for reactive, self-gravitating flow and their numerical solution in the context of core-collapse supernovae and convective burning in Sect. 2. We do not treat numerical methods for MHD in supernova simulations, although we occasionally comment on the role of MHD effects in the later sections. In the subsequent sections, we then review recent simulation results and progress in the theoretical understanding of convection during the late burning stages (Sect. 3), of supernova shock revival (Sect. 4), and the hydrodynamics of the explosion phase (Sect. 5).

Modeling the late stages of nuclear burning and the subsequent supernova explosion involves solving the familiar equations for mass, momentum, and energy conservation with source terms that account for nuclear burning and the exchange of energy and momentum with neutrinos. Viscosity and thermal heat conduction mediated by photons, electron/positrons, and ions can be neglected, and so we have (in the Newtonian limit and neglecting magnetic fields),

Since Godunov-based finite-volume solvers are now most commonly employed for simulating core-collapse supernovae and the late burning stages, we shall focus on the problem-specific challenges for this approach in this section.

Grid-induced perturbations Cartesian grids have the virtue of algorithmic simplicity and do not suffer from coordinate singularities, but also come with disadvantages as they are not adapted to the approximate symmetry of the physical problem. The unavoidable non-spherical perturbations from the grid make it impossible to reproduce the spherically symmetric limit in multi-D even for perfectly spherical initial conditions, or to study the growth of non-spherical perturbations in a fully controlled manner. The former deficiency is arguably an acceptable sacrifice, though it can limit the possibilities for code testing and verification, but the latter can introduce visible artifacts in simulations. For example, Cartesian codes sometimes produce non-vanishing gravitational wave signals from the bounce of non-rotating cores (Scheidegger et al. 2010), and often show dominant \(\ell =4\) modes during the growth phase of non-radial instabilities (Ott et al. 2012).

In spherical coordinates, the multi-scale nature of the problem can be accommodated to a large degree by employing a non-uniform radial grid that transitions to roughly equal spacing in \(\log r\) at large radii. Radial resolution can be added selectively in strongly stratified regions like the PNS surface (Buras et al. 2006b), or one can use an adaptive moving radial grid (Liebendrfer et al. 2004; Bruenn et al. 2018). However, some care must be exercised in using non-uniform radial grids. Rapid variations in the radial grid resolution \(\varDelta r/r\) can produce artifacts such as artificial waves and disturbances of hydrostatic equilibrium.

It is also straightforward to implement a moving radial grid to adapt to changing resolution requirements or the bulk contraction/expansion of the region of interest; see Winkler et al. (1984), Mller (1994) for an explanation of this technique. The MPA/Monash group routinely apply such a moving radial grid in quasi-Lagrangian mode during the collapse phase (Rampp and Janka 2000), and, with a prescribed grid function, in parameterized multi-D simulations of neutrino-driven explosions (Janka and Mller 1996; Scheck et al. 2006) and in simulations of convective burning (Mller et al. 2016b). The Oak Ridge group uses a truly adaptive radial grid in their supernova simulations with the Chimera code (Bruenn et al. 2018). A moving radial mesh might also appear useful for following the expansion of the ejecta and the formation of a strongly diluted central region in simulations of mixing instabilities in the envelope, but the definition of an appropriate grid function is non-trivial. Most simulations of mixing instabilities in spherical polar coordinates have therefore relied on simply removing zones continuously from the evacuated region of the blast wave to increase the time step (Hammer et al. 2010) rather than implementing a moving radial mesh (Mller et al. 2018).

Both fixed mesh refinement and spherical grids with non-uniform radial mesh spacing only provide limited adaptability to the structure of the flow. Truly adaptive mesh refinement can provide superior resolution in cases where very tenuous, non volume-filling flow structures emerge. Mixing instabilities in the envelope are a prime example for such a situation, and have often been studied using AMR in spherical polar coordinates (Kifonidis et al. 2000, 2003, 2006) in 2D and Cartesian coordinates (Chen et al. 2017).

Even with a 1D treatment for the innermost grid zones, one is still left with a severe time-step constraint at the grid axis in 3D. A number of alternatives to spherical polar grids with uniform spacing in latitude \(\theta \) and longitude \(\varphi \) can help to remedy this. The simplest workaround is to adopt uniform spacing in \(\mu =\cos \theta \) instead of \(\theta \). In this case, one has \(\sin \theta = (2N_\theta -1)^1/2/N_\theta \approx \sqrt2 N_\theta ^-1/2\) in the zones adjacent to the axis for \(N_\theta \) zones in latitude instead of \(\sin \theta \approx N_\theta ^-1/2\), so the time step limit scales as \(\varDelta t \propto \sqrt2 N_\theta ^-1/2 N_\phi ^-1\) instead of \(\varDelta t \propto N_\theta ^-1 N_\phi ^-1/2\), where \(N_\varphi \) is the number of zones in longitude. Alternatively, one can selectively increase the \(\theta \)-grid spacing in the zones close to the axis. However, the time step constraint at the axis is still more restrictive than at the equator in this approach, and the aspect ratio of the grid cells becomes extreme near the pole, which can create problems with numerical stability and accuracy.

Filtering in Fourier space, which has long been used in the meteorology community (Boyd 2001), provides another means of curing the restrictive time step constraint near the axis and has been implemented in the CoCoNuT-FMT code (Mller et al. 2019). This can also be implemented with minimal interventions in a solver for spherical polar grids, and is attractive because the amount of smoothing that is applied to the solution increases more gradually towards the axis than with mesh coarsening schemes. Mller et al. (2019) suggests that this eliminates the problem of axis-aligned flows. On the downside, simulations with Fourier filtering may occasionally encounter problems with the Gibbs phenomenon at the shock.

3a8082e126