A comparison of 2D Magnetohydrodynamic supernova simulations with the CoCoNuT-FMT and Aenus-Alcar codes

ABSTRACT

1 INTRODUCTION

It was first recognized by Baade & Zwicky (1934) that in supernova explosions of massive stars (⁠|$\mathord {\ge } 8 \, \mathrm{M_\odot }$|⁠) a small fraction of the energy released during the terminal collapse of the stellar iron core is somehow transferred to the stellar envelope to produce a powerful explosion known as a core-collapse supernova (CCSN). The ensuing quest for the details of the explosion mechanism has been ongoing for decades. Among the many suggested explosion scenarios, the neutrino-driven mechanism, initially proposed by Colgate & White (1966) has emerged as the favoured explanation for the majority of CCSN explosions. In the modern paradigm of the delayed neutrino-driven mechanism (Bethe & Wilson 1985), shock revival is accomplished thanks to the increase of the post-shock pressure due to partial reabsorption of neutrinos that stream out from the young proto-neutron star (PNS) and the cooling layer of accreted material on its surface. Detailed simulations have determined that additional support by multidimensional instabilities like turbulent convection (Burrows, Hayes & Fryxell 1995; Herant 1995) and the standing accretion shock instability (Blondin, Mezzacappa & DeMarino 2003) is critical to make the neutrino-driven mechanism work.

Partly due to increasing computing power, supernova simulations have been able to continually include more detailed physics such as multigroup neutrino transport, magnetic fields, and general relativity, and to abandon symmetry assumptions in favour of fully three-dimensional models, but even the most sophisticated simulations still need to make approximations (for recent reviews, see Janka 2012; Burrows 2013; Müller 2016, 2020; Mezzacappa et al. 2020). The complexity of the codes and the disparities in numerical methodology present a challenge for drawing robust conclusions on the workings of the supernova engine from the models. Code comparison studies are a useful approach for increasing the robustness and reliability of simulation results and have already addressed many important ingredients of supernova codes, some focusing more strictly on the neutrino transport (e.g. Yamada, Janka & Suzuki 1999; Ott et al. 2008; Richers et al. 2017; Chan, Müller & Heger 2020), whereas others have more broadly addressed the interplay of neutrino transport and hydrodynamics in 1D, 2D (Liebendörfer et al. 2005; Müller, Janka & Dimmelmeier 2010; Just et al. 2018; O’Connor et al. 2018) and to some extent in 3D (Cabezón et al. 2018) with different solvers for the hydrodynamics and the neutrino transport alike. Numerical uncertainties related to the effects of grid resolution have also received considerable attention (e.g. Hanke et al. 2013; Couch & O’Connor 2014; Handy, Plewa & Odrzywołek 2014; Abdikamalov et al. 2015; Radice et al. 2016; Nagakura et al. 2019; Melson, Kresse & Janka 2020).

A similarly broad exploration of numerical uncertainties is yet to be performed for magnetorotational explosions, which constitute the best-explored alternative to the neutrino-driven mechanism and are widely thought to operate in hyperenergetic supernovae (‘hypernovae’), some of which are also associated with long gamma-ray bursts (GRBs; Woosley & Bloom 2006). In this scenario, whose history can be traced to early pioneering work of LeBlanc & Wilson (1970), Bisnovatyi-Kogan, Popov & Samokhin (1976), Meier et al. (1976), and Müller & Hillebrandt (1979), rapid rotation and strong magnetic fields are the critical ingredients for the explosion; the magnetic fields help tap the rotational energy of the forming compact object and channel it into kinetic energy in outflows – usually in the form of collimated, magnetically dominated jets – which are driven either by magnetic forces and pressure or through viscous heating behind the shock.

With the inclusion of rotation and magnetic fields in models, many additional unknowns enter into the simulations and have received substantial attention. The magnetorotational explosion mechanism is intimately linked to the problem of rotation and magnetism in massive stars. Initial conditions for magnetorotational explosion simulations currently come from ‘1.5D’ stellar evolution models that assume shellular rotation and include effective recipes for magnetic field generation and angular momentum transport by hydrodynamic and magnetohydrodynamic processes (Heger, Langer & Woosley 2000; Heger, Woosley & Spruit 2005; Woosley & Heger 2006). Efforts to better understand the angular momentum distribution and magnetic fields in the cores of massive stars by means of multidimensional simulation have only started recently (McNeill & Müller 2021; Varma & Müller 2021; Yoshida et al. 2021). Since 1.5D stellar evolution models predict magnetic fields that are rather weak and pre-dominantly toroidal, the general notion has long been that field amplification processes after collapse are critical in magnetorotational explosions, although this has recently been challenged (Obergaulinger & Aloy 2017, 2021; Aloy & Obergaulinger 2021). These amplification processes have duly received considerable attention. Akiyama et al. (2003) pointed out the potential importance of the magnetorotational instability (MRI; Balbus & Hawley 1991), which has been addressed in the context of CCSNe in a host of numerical and analytic studies (Obergaulinger et al. 2009; Sawai, Yamada & Suzuki 2013; Guilet & Müller 2015; Guilet, Müller & Janka 2015; Masada, Takiwaki & Kotake 2015; Mösta et al. 2015; Sawai & Yamada 2016; Rembiasz et al. 2016a,b; Reboul-Salze et al. 2021). Field amplification by an α − Ω dynamo in the PNS (Thompson & Duncan 1993) has also been investigated in multi-D simulations as an avenue for generating strong magnetic fields on large scales from weak seed fields (Raynaud et al. 2020). While important issues like the initial fields and angular momentum distribution, the saturation of the MRI, and non-ideal effects still warrant further investigation, global simulations of magnetorotational explosions based on parameterized initial fields have been available for quite some time in 2D both with (Burrows et al. 2007; Dessart et al. 2007; Obergaulinger & Aloy 2017, 2020; Bugli et al. 2020; Aloy & Obergaulinger 2021; Jardine, Powell & Müller 2021) and without (Kotake et al. 2004; Takiwaki et al. 2004; Sawai, Kotake & Yamada 2005; Moiseenko, Bisnovatyi-Kogan & Ardeljan 2006; Obergaulinger, Aloy & Müller 2006; Takiwaki, Kotake & Sato 2009) detailed neutrino transport, and are becoming routine in 3D (Winteler et al. 2012; Mösta et al. 2014, 2018; Halevi & Mösta 2018), even though only few models (Kuroda et al. 2020; Bugli, Guilet & Obergaulinger 2021; Obergaulinger & Aloy 2021) also include neutrino transport.

As magnetohydrodynamic supernova simulations mature and are being used to understand hypernova light curves (Shankar et al. 2021) and explosion properties and the potential role of hypernovae as a site for the rapid-neutron capture process (Winteler et al. 2012; Mösta et al. 2018; Reichert et al. 2021), it is critical to better address numerical uncertainties in these models following a similar approach as for neutrino-driven explosions, where comparison studies with different codes and resolutions studies have done much to distinguish between robust and uncertain features of the mechanism. However, to date, no study has as yet compared simulations of magnetorotational supernovae using independent codes.

In this paper, we present the first comparison of axisymmetric (2D) magnetorotational explosions simulations between two independent codes, namely the CoCoNuT-FMT (Müller & Janka 2015; Müller & Varma 2020) and Aenus-Alcar (Obergaulinger 2008; Obergaulinger, Janka & Aloy 2014) codes. We also conduct a minimal resolution study with CoCoNuT-FMT models and explore the sensitivity of the results to the approximate Riemann solver. Our paper is organised as follows: In Section 2, we describe the magnetohydrodynamic (MHD) solvers in CoCoNuT-FMT and Aenus-Alcar and also provide a brief description of the neutrino transport schemes. We then outline the simulation set-up in Section 3. We compare the CoCoNuT-FMT and Aenus-Alcar results in Section 4 and present results for our minimal resolution study with CoCoNuT-FMT in Section 5. We conclude with a summary and evaluation in Section 6.

2 NUMERICAL METHODS

We perform axisymmetric core-collapse supernova simulations with the pseudo-Newtonian MHD version of CoCoNuT-FMT (Müller et al. 2010; Müller & Janka 2015) and the special relativistic Aenus-Alcar code (Obergaulinger 2008; Just, Obergaulinger & Janka 2015). In both codes, we use an effective gravitational potential (Case A from Marek et al. 2006) to approximately capture general relativistic effects. The nuclear equation of state of Lattimer & Swesty (1991) with an incompressibility of |$K=220 \, \mathrm{MeV}$| is used for all models at high densities. Here, we briefly outline the magnetohydrodynamic and neutrino transport solvers in the two codes. For more details, we refer to the relevant code papers.

2.1 CoCoNuT-FMT code

In the CoCoNuT models the innermost |$10\, \mathrm{km}$| are treated in 1D to avoid an overly stringent Courant time-step limit near the origin of the grid. The thermodynamic variables and the radial velocity are assumed to be spherically symmetric. The non-radial velocity components and the magnetic fields are not set to zero in the core, however. Instead, we assume uniform rotation at a rate corresponding to the total angular momentum in the core, and the magnetic field in the core is represented by vector spherical harmonics of degree ℓ = 1, whose coefficients are evolved according to the induction equation.

The fast multigroup transport (FMT) scheme of Müller & Janka (2015) is used for the neutrino transport. The FMT scheme solves the energy-dependent stationary zeroth-momentum equation using a hybrid closure based on a two-stream Boltzmann solution and an analytic Eddington closure at flux factors h > 0.5. The scheme accounts for gravitational redshift, but mostly disregards velocity-dependent terms, except for an effective frame correction of the absorption opacity at low optical depth (Müller et al. 2019). The set of neutrino rates includes emission by and absorption on nucleons and nuclei, neutrino scattering on nucleons and nuclei, and, for heavy-flavour neutrinos, Bremsstrahlung in an effective one-particle approximation. Recoil energy transfer in neutrino–nucleon scattering is included approximately. Nucleon potentials are taken into account in charged-current neutrino–nucleon interactions (Martínez-Pinedo et al. 2012). Nucleon correlations are included following Horowitz et al. (2017). Neutral-current scattering cross-sections are computed following Horowitz (2002) with a somewhat large value for the nucleon strangeness of g_A = −0.1. The deleptonization scheme of Liebendörfer (2005) is used during the collapse phase.

2.2 Aenus-Alcar code

The Aenus-Alcar code solves the special relativistic version of the MHD equations. Different from CoCoNuT, the divergence constraint is accounted for by using a constrained transport scheme (Evans & Hawley 1988). The equations of MHD are also solved using the HLLC solver, but combined with MP5 reconstruction (Suresh & Huynh 1997). The code can accommodate different coordinate systems, but is also run in spherical polar coordinates in this study. Aenus-Alcar treats the innermost |$10\, \mathrm{km}$| in 2D, but with a coarser grid resolution.

Neutrinos are treated by a spectral scheme for two-moment transport with an analytic closure (see Just et al. 2015 for an overview). Neutrinos of different energies are coupled via the effects of Doppler shifts, aberration, and gravitational redshift treated to |$\mathcal {O}(v/c)$|⁠. The set of reactions included in the Aenus-Alcar simulations is given in detail in, e.g. Obergaulinger & Aloy (2017), Obergaulinger, Just & Aloy (2018), and consists of emission/absorption and scattering reactions with nucleons and nuclei, inelastic scattering off electrons, electron–positron annihilation, and nucleon–nucleon Bremsstrahlung.

We do not aim to reproduce perfectly identical neutrino microphysics in both codes. We rather pit the two codes against each other as used in normal production runs so that the comparison is more reflective of ‘real-life’ uncertainties encountered in magnetorotational supernovae. To give a conservative assessment of possible errors, we view it as one useful strategy to overestimate rather than to underestimate code differences by not making the models more similar than they would be in independent studies. Of course, this also entails some risk of accidental error cancellation in some (but not all) metrics of the code comparison. This strategy is complementary to the opposite approach of reproducing all of the physics as closely as possible in different codes, which will be a useful undertaking for the future, perhaps within a comparison broader with a broader community base.

3 MODEL SET-UP

All our simulations use the stellar progenitor model 35OC of Woosley & Heger (2006). This model was computed using the Kepler code including a prescription for angular momentum transport by magnetic torques (Spruit 2002; Heger et al. 2005). The 35OC stellar model is a rapidly rotating model with an initial total angular momentum of |$J_\mathrm{init} = 10^{52}\, \mathrm{erg}\, \mathrm{s}$| and a surface rotational velocity |$v_\mathrm{rot} = 380\, \mathrm{km/s}$| on the zero-age main sequence (ZAMS). The progenitor has a ZAMS mass of M_ZAMS = 35M_⊙, but at the time of core collapse mass-loss has reduced its total mass to M_final = 28.07M_⊙.

All the CoCoNuT-FMT models (denoted by a prefix ‘coco’) have a computational domain that extends out to |$R_\mathrm{out} = 10^{5}\, \mathrm{km}$| using a non-uniform radial grid of 550 zones and spans colatitudes θ in the range [0^○, 180^○] on a spherical polar grid with 128 zones in θ. The Aenus-Alcar models (denoted by a prefix ‘Alcar’) implement the same angular resolution, but have a computational domain of 400 radial zones that extends to |$R_\mathrm{out} =1.4 \times 10^{5}\, \mathrm{km}$|⁠.

We simulate three pairs of models with CoCoNuT-FMT. For each strong- and weak-field model, we conduct runs using the MHD version of the HLLC Riemann solver (Gurski 2001) or the HLLE solver (Einfeldt 1988). These runs are labelled as coco-HLLC-Rs/w and coco-HLLE-Rs/w, respectively. In addition, we conduct two high-resolution runs coco-Res-Rs/w with the HLLC solver with 256 zones in angle for a minimal resolution test. The parameters of all runs are summarized in Table 1.

4 RESULTS – CODE COMPARISON

In the following, we present a comparison of the standard resolution models coco-HLLC-Rs, coco-HLLC-Rw, coco-HLLE-Rs, and coco-HLLE-Rw to models Alcar-Rs and Alcar-Rw as described in Table 1. We discuss the differences between the Aenus-Alcar and CoCoNuT-FMT results, as well as between the HLLC and HLLE Riemann solvers implemented in CoCoNuT-FMT.

4.1 Shock propagation

We begin by comparing the propagation of the shock (Fig. 1) between the different models, as this is one of the most basic metrics for the post-bounce dynamics. Due to the strong asymmetry of the explosion, the dynamics is best captured by the maximum shock radius rather than the angle-averaged shock radius, especially in the case of the strong magnetic field models, which develop a veritable jet-driven explosion with strongly collimated polar outflows.

Figure 1.

Evolution of the maximum shock radius at standard resolution for CoCoNuT-FMT using the HLLC and HLLE Riemann solver models along with the Aenus-Alcar models. Both the strong-field (-Rs) and weak-field (-Rw) models are shown.

Open in new tabDownload slide

In Alcar-Rs, shock expansion starts promptly after core bounce and does not exhibit a stalled shock phase. By contrast, both coco-HLLC-Rs and coco-HLLE-Rs stall for a short period of |$\mathord {\approx } 50\, \mathrm{ms}$| before the shock is revived. After shock revival, the shock in models coco-HLLC-Rs and coco-HLLE-Rs expands rapidly, achieving similar maximum shock radii as Alcar-Rs at the end of the simulations.

The different time of explosion in Alcar and CoCoNuT may be due to different dynamics of magnetic field amplification in the PNS surface region that result in the earlier emergence of a jet in model Alcar-Rs. There may, however, be a compounding factor in the form of the collapse dynamics. Both in the strong-field case and the weak-field case, the Alcar models take |$30{-}35\, \mathrm{ms}$| longer to collapse. The mass accretion rate at |$45{-}130\, \mathrm{ms}$| differs noticeably between CoCoNuT and Alcar with CoCoNuT exhibiting higher mass accretion rates early on after collapse (Fig. 2). The collapse time is known to be sensitive to the microphysics (Hüdepohl 2014; Powell, Müller & Heger 2021), and although both codes use the same high-density equation of state, they differ in their treatment of the low-density equation of state, nuclear burning, and deleptonization during collapse. Factors that affect the collapse dynamics may be exacerbated by transient adjustment of the stellar models after mapping the 1D stellar evolution model with shellular rotation into the supernova simulation codes. Although the resulting difference of the infall dynamics are mostly confined to the collapse phase, it is known that differences in the mass accretion rate can persist for some time after collapse even in simulations of non-rotating progenitors (Hüdepohl 2014). In the subsequent analysis, one must bear in mind that this will inevitably impact the shock propagation as well as the evolution of the PNS properties, and, more importantly, the neutrino emission at early post-bounce times to some degree.

Figure 2.

Mass accretion rate |$\dot{M}$| as a function of post-bounce time, measured as the time derivative of the mass of the PNS (defined as the region where |$\rho \gt 10^{11} \, \mathrm{g}\, \mathrm{cm}^{-3}$|⁠).

Open in new tabDownload slide

Unlike the strong-field models, the models with weaker magnetic fields (coco-HLLC-Rw, coco-HLLE-Rw, Alcar-Rw) all go through a prolonged accretion phase, similar to the dynamics in purely neutrino-driven explosions of massive progenitors. The shock trajectories during the accretion phase are similar, with Alcar-Rw exhibiting slightly larger stalled shock radii early on, but slightly slower shock expansion later. Most of the time, the relative difference between the runs remains below |$10{{\ \rm per\ cent}}$|⁠. All three models undergo shock revival around |$0.4\, \mathrm{s}$| post-bounce with minor differences. The Alcar-Rw model is the last to explode and continues to have a smaller shock radius than both CoCoNuT-FMT models until about |$0.8\, \mathrm{s}$| where its shock radii slightly surpasses that of coco-HLLC-Rw. At their largest divergence, the relative difference between the three models is |$\mathord {\approx } 12\, {{\ \rm per\ cent}}$|⁠. We also note that the CoCoNuT-FMT weak-field models show a much more pronounced and prolonged activity of the standing accretion shock instability (SASI; Blondin et al. 2003; Foglizzo et al. 2007; Guilet & Foglizzo 2012, visible as oscillations in the maximum shock radius in coco-HLLC-Rw and coco-HLLE-Rw in the phase leading up to shock revival. This is in line with the established finding that stronger shock retraction is conducive to more vigorous SASI activity (Scheck et al. 2008; Müller, Janka & Heger 2012).

Shock revival in model coco-HLLE-Rw occurs about 0.05s later than in model coco-HLLC-Rw, but the shock subsequently expands faster than in the HLLC model, exhibiting a larger maximum shock radius at the end of the simulation runs. By contrast, no significant differences are seen between coco-HLLE-Rs and coco-HLLC-Rs. Since shock trajectories differ between the HLLE and HLLC run only in the weak-field case, which develops a ‘hybrid’ explosion supported alike by magnetic fields, neutrino heating, and convection (Obergaulinger & Aloy 2017), we suspect the differences between coco-HLLE-Rw and coco-HLLC-Rw do not actually reflect a systematic dependence on the Riemann solver, but are due to stochastic variations, which have been shown to be considerable in explosion models without MHD in 2D (Cardall & Budiardja 2015; Müller 2015).

Overall, the differences in shock propagation between the corresponding CoCoNuT-FMT and Aenus-Alcar models are within the range of deviations between independent codes in comparisons of neutrino hydrodynamics simulations in 2D (Müller & Janka 2015; Cabezón et al. 2018; Just et al. 2018). Even in controlled 1D comparisons with identical neutrino rates (Liebendörfer et al. 2005; Müller et al. 2010; O’Connor et al. 2018), deviations of the order of |$10{-}20{{\ \rm per\ cent}}$| in shock radius are not uncommon. We suspect that the prompt explosion in Alcar-Rs that is not mirrored in the strong magnetic field models of CoCoNuT-FMT is an artefact of the different treatments of the inner core. As described in Section 3, Aenus-Alcar treats the inner core in coarse 2D geometry, which allows for rotational shear and hence more rapid magnetic field amplification immediately after bounce, and consequently, larger magnetic pressure to immediately drive the shock out. These difference in magnetic field amplification will be discussed further in Section 4.3.

4.2 Explosion energetics

Figure 3.

Evolution of the diagnostic explosion energy E_expl for CoCoNuT-FMT and Aenus-Alcar models with strong and weak initial magnetic fields.

Open in new tabDownload slide

Although the CoCoNuT-FMT and Aenus-Alcar models show qualitative agreement in the evolution of E_expl, Fig. 3 more clearly exhibits quantitative differences between the runs than could be seen in the shock propagation. Among the strong-field models, models coco-HLLE-Rs and coco-HLLC-Rs agree quite well with each other, but develop a positive diagnostic energy later than the Alcar-Rs run, which is consistent with the longer delay to shock revival discussed in the previous section. After this slight delay, model coco-HLLC-Rs in particular exhibits comparable rates of explosion energy growth, with Alcar-Rs increasing at |$\mathord {\approx } 3.12 \times 10^{51}\, \mathrm{erg}\, \mathrm{s}^{-1}$| while model coco-HLLC-Rs shows, on average, a rise at a rate of |$\mathord {\approx } 3.17 \times 10^{51}\, \mathrm{erg}\, \mathrm{s}^{-1}$|⁠. The explosion in coco-HLLE-Rs initially grows at a similar rate as in coco-HLLC-Rs, but the model becomes slightly less energetic by the end of the simulation. The cause of this divergence is unclear but since our simulations end by |$500 \, \mathrm{ms}$|⁠, it is possible that this is just a temporary deviation due to some element of stochasticity in the dynamics of the evolving explosion. The diagnostic energy is evidently a much more sensitive measure for the dynamics than the shock radius. The energetics of the narrowly collimated jets that drive shock expansion is sensitive to a number of details, e.g. how efficiently energy is extracted from differential rotation, and how much material is entrained by the collimated jets that drive the explosion. It should also be noted that relatively minor differences in the kinetic, magnetic, and thermal energy of the ejecta can translate into larger differences once the potential energy is subtracted. For these reasons, the rather large discrepancy in explosion energy between the Aenus-Alcar and CoCoNuT-FMT models is not surprising. In Section 4.3, we shall demonstrate that the two codes produce somewhat different PNS structures that can motivate (though not rigorously explain) significant differences in the energetics of the jet-driven explosions. We will also revisit the sensitivity of the explosion energy in our resolution analysis in Section 5.

The weak-field models appear to be qualitatively more disparate than their strong-field counterparts. All three models have significant delays of |$\mathord {\sim } 400 \, \mathrm{ms}$| to the explosion with a slightly longer delay phase for the Alcar-Rw model, mirroring the different time of shock revival (Fig. 1). Although the shock radii of all three models are comparable by the end of the simulations, we find pronounced quantitative differences in their explosion energies; the explosion energy of model Alcar-Rw is |$80{{\ \rm per\ cent}}$| lower than that of coco-HLLC-Rw, while model coco-HLLE-Rw ends up with a significantly higher explosion energy than the other two models by the end of the simulations. This does not indicate fundamental differences between the codes, but is likely just a reflection of the strongly stochastic explosion dynamics of 2D models in the ‘hybrid’ regime where neutrino heating, convection or SASI, and magnetic fields jointly drive the explosion. The unsteady, and even non-monotonic growth of the explosion energy in coco-HLLC-Rw and Alcar-Rw clearly reveals this stochasticity, and the variations between the different runs are similar to the variations seen in non-magnetic models for different realizations of the same model set-up with the same code (Cardall & Budiardja 2015; Müller 2015). This does not, of course, preclude that there are other systematic differences in the dynamics between the weak-field runs that are hidden behind the stochastic variations.

4.3 Proto-neutron star parameters and magnetic field amplification

We next consider the bulk parameters and structure of the PNS, and the amplification of magnetic fields inside it as key drivers of the overall dynamics. Fig. 4 shows the PNS mass M_PNS, angular momentum J_PNS, and magnetic energy E_mag, which we compute as volume integrals over the region where |$\rho \gt 10^{11} \, \mathrm{g}\, \mathrm{cm}^{-3}$| following the customary, though somewhat arbitrary definition of the PNS surface. The evolution of the angle-averaged PNS radius R_PNS for the different runs is plotted in Fig. 5. Angle-averaged radii are determined by identifying the PNS surface using the same fiducial density |$\rho =10^{11} \, \mathrm{g} \, \mathrm{cm}^{-3}$| based on angle-averaged density profiles.

Figure 4.

PNS mass M_PNS, angular momentum J_PNS, and magnetic energy E_mag, as a function of post-bounce time for the standard-resolution models. The PNS is defined as the region where |$\rho \gt 10^{11} \, \mathrm{g}\, \mathrm{cm}^{-3}$|⁠.

Open in new tabDownload slide

Figure 5.

Comparison of the time evolution of the angle-averaged PNS radius R_PNS between the CoCoNuT-FMT and Aenus-Alcar models.

Open in new tabDownload slide

The growth of the PNS mass (top panel in Fig. 4) is very similar in the weak-field models Alcar-Rw, coco-HLLC-Rw, and coco-HLLE-Rw runs. We do, however, see an extended period before 0.35s, where the CoCoNuT-FMT models have a more massive PNS. The Alcar-Rw model then catches up after |$0.35\, \mathrm{s}$|⁠, and eventually ends up with a slightly more massive PNS. Among the strong-field models, coco-HLLC-Rs and coco-HLLE-Rs exhibit more substantial accretion on to the PNS immediately after collapse compared to model Alcar-Rs, leading to a  |$9{{\ \rm per\ cent}}$| higher PNS mass in the CoCoNuT-FMT runs throughout. This is explained by the more rapid onset of the explosion seen in Alcar-Rs. The PNS mass difference largely accounts for a similar difference in PNS angular momentum between the CoCoNuT-FMT and Alcar strong-field models (Fig. 4, middle panel).

More dramatic differences between the two codes can be seen in the total PNS angular momenta after about |$0.3\, \mathrm{s}$| post-bounce for the weak-field models. Models coco-HLLC-Rw and coco-HLLE-Rw start to gradually lose angular momentum, whereas model Alcar-Rw continues to slowly gain angular momentum as it accretes mass. This difference is especially noteworthy since the PNS masses of the weak-field models are quite similar throughout. This behaviour in the CoCoNuT-FMT models is likely due to spurious non-conservation of angular momentum due to our numerical implementation in axisymmetry, where the velocity component v_φ is zeroed in the zones adjacent to the grid axis. Aside from these boundaries, we shall see that the different treatments of the inner PNS lead to significant differences in the angular momentum profiles between the two codes.

More pronounced differences between the CoCoNuT-FMT and Aenus-Alcar models than in the PNS mass and angular momentum are found in the PNS magnetic energy E_mag (the bottom panel of Fig. 4). In coco-HLLC-Rs and coco-HLLE-Rs, the magnetic energy in the PNS first peaks around |$1.25\times 10^{51}\, \mathrm{erg}$| at |$0.1\, \mathrm{s}$| and remains high for some time, decaying after a second peak. Although the overall time evolution of E_mag in these two models is quite similar, coco-HLLE-Rs is set apart by an earlier gradual drop in magnetic energy around 0.1s before coco-HLLC-Rs. This is noteworthy, considering that the choice of Riemann solver shows very little effect on the other PNS quantities. In model Alcar-Rs, the magnetic field energy peaks immediately after core bounce at |$\mathord {\approx } 1.05\times 10^{51}\, \mathrm{erg}$|⁠, but then decreases sharply to a much lower level.

The weak magnetic field models exhibit similarly large differences as their strong-field counterparts. In particular, we see a much earlier rise in PNS magnetic energy in Alcar-Rw than in coco-HLLC-Rw and coco-HLLE-Rw, followed by a slow, gradual drop after |$\mathord {\sim } 0.5\, \mathrm{s}$|⁠. By contrast, the weak-field runs with CoCoNuT-FMT show a gradual rise in magnetic energy over time that continues until the end of the simulation.

Figure 6.

Top to bottom: Profiles of the energy density contained in the radial, meridional, and zonal components B_r, B_θ, and B_φ of the magnetic field for the strong-field models coco-HLLC-Rs (solid lines) and Alcar-Rs (dashed lines) at various post-bounce times.

Open in new tabDownload slide

Figure 7.

Top to bottom: Profiles of the energy density contained in the radial, meridional, and zonal components B_r, B_θ, and B_φ of the magnetic field for the weak-field models coco-HLLC-Rw (solid lines) and Alcar-Rw (dashed lines) at various post-bounce times.

Open in new tabDownload slide

Figure 8.

Spherically averaged specific angular momentum as a function of radius at different post-bounce times for the strong-field models Alcar-Rs (dotted lines) and coco-HLLC-Rs (solid lines).

Open in new tabDownload slide

Figure 9.

Spherically averaged specific angular momentum as a function of radius at different post-bounce times for the weak-field models Alcar-Rw (dotted lines) and coco-HLLC-Rw (solid lines).

Open in new tabDownload slide

Both for the weak- and strong-field case, the profiles of e_B,r, e_B,θ, and e_B,φ are very similar at bounce, i.e. field amplification during the collapse phase in the inner core is consistent between the two codes despite the different treatment of the induction equation. Later, however, the situation changes for two reasons.

In the strong-field case, we first see a steady drop of magnetic field energy in the innermost few kilometres in the ACLAR-Rs run because the inner core is treated in low-resolution 2D, leading to more numerical diffusion. In model coco-HLLC-Rs, strong field dissipation in the core is avoided by enforcing uniform rotation and spherical symmetry in the meridional flow components in the innermost |$10\, \mathrm{km}$|⁠. Secondly, we see stronger field amplification in model coco-HLLC-Rs around a radius of |$20\, \mathrm{km}$|⁠, primarily in e_B,φ and to a lesser extent in e_B,r, which is indicative of strong shear amplification at the PNS surface. This is consistent with the finding that the angular momentum profiles show a hump around |$20\, \mathrm{km}$| in model coco-HLLC-Rs, whereas the specific angular momentum plateaus outside the PNS than in Alcar-Rs (Fig. 8). This implies a steeper negative angular velocity gradient, i.e. stronger shear, in around and outside |$20\, \mathrm{km}$| in coco-HLLC-Rs. The stronger fields generated around these radii in coco-HLLC-Rs also translate into stronger fields in the jet in the long term; the components e_B,r and e_B,φ are about a factor of 10 higher at the maximum shock radius at late times. The profiles of e_B,r reveal a numerical disadvantage of the divergence cleaning scheme as originally implemented in CoCoNuT-FMT; choosing the maximum cleaning speed allowed by the CFL limit results in spurious transport of magnetic energy into the pre-shock region. The different treatment of the inner core in both codes has further implications for the long-term evolution of magnetic fields in the PNS. While Alcar-Rs initially loses magnetic field energy in the innermost zones, it also shows stronger field amplification in e_B,θ just inside |$10\, \mathrm{km}$| in contrast to CoCoNuT-HLLC-Rs. This is due to the inclusion of differential rotation of the inner core in the Alcar code, which allows for field amplification in this region that is responsible for the rise in the PNS magnetic energy E_mag after the very rapid decreases in magnetic energy in e_B,r, e_B,θ, and e_B,φ from numerical dissipation in the innermost zones (cp. Fig. 4).

The different evolution of the PNS magnetic energy in the weak-field case (with a faster rise in Aenus-Alcar than in CoCoNuT-FMT) can also be understood as a consequence of the different core treatment. Model Alcar-Rw again shows field amplification by differential rotation inside the core, whereas CoCoNuT-HLLC-Rw only shows a slow decline in e_B,r and e_B,θ and a slow rise in e_B,φ by compression inside the innermost |$10 \, \mathrm{km}$|⁠. Interestingly, the (putatively numerical) decline of the magnetic field energy in the innermost zones is less pronounced in model Alcar-Rw than in model Alcar-Rs, suggesting that this phenomenon is very sensitive to the dynamics near the origin of the grid. Field amplification in the differentially rotating PNS core generates substantial toroidal fields, which dominate the PNS magnetic energy; ultimately, the toroidal component of the magnetic energy density is increased relatively uniformly out to a radius of about |$20{-}30\, \mathrm{km}$|⁠. By contrast, field amplification in model coco-HLLC-Rw is confined to the PNS surface region outside the uniformly rotating core at radii of |$10{-}20\, \mathrm{km}$|⁠. Similar to the strong-field case, the amplification of the toroidal field by shear in this region is stronger than in the Alcar model. Due to the stronger field amplification at the PNS surface, the PNS magnetic energy in model CoCoNuT-HLLC-Rw eventually overtakes model Alcar-Rw at late times, overcompensating for the spurious lack of field amplification inside the PNS.

The angular momentum profiles of models coco-HLLC-Rs and Alcar-Rs (Fig. 8) and of models coco-HLLC-Rs and Alcar-Rs (Fig. 9) yield more insights into the causes of the different patterns of field amplification in the CoCoNuT and Alcar models. The strong-field model coco-HLLC-Rs shows the build-up of a hump in specific angular momentum outside the PNS centred at radii of |$40{-}50\, \mathrm{km}$| in the late-time profiles. The presence of such a hump implies considerable shear in and above the PNS surface region, which explains the strong magnetic fields generated in this layer in model coco-HLLC-Rs (Fig. 6). In model Alcar-Rs, the hump is less pronounced, more smeared out, and located further outside. In the early profiles, where the hump is yet absent, the transition to a plateau in specific angular still occurs further outside in model Alcar-Rs compared to coco-HLLC-Rs, implying less shear and less field amplification in the PNS surface region. Further analysis reveals that the different angular momentum distribution in the two runs is related to the different PNS radii seen in Fig. 5; effectively the PNS has more fluffy, disc-like surface in model Alcar-Rs. It is difficult to decide a priori which of the codes captures the structure of the disc-like surface region more accurately. We shall see in Section 5 that CoCoNuT-FMT also tends to produce more a more fluffy PNS surface at higher resolution. This tentatively suggests that the more fluffy PNS surface structure with a flatter angular momentum distribution is more realistic.

The angular momentum profiles deeper in the PNS also reveal an interesting phenomenon. In model Alcar-Rs, there is a clear secular drop in angular momentum at radii of less than |$10\, \mathrm{km}$|⁠, whereas model coco-HLLC-Rs shows little changes in the angular momentum profiles in this region. One might be tempted to connect the slow-down of PNS rotation with numerical diffusion of angular momentum into the innermost zone, where the angular momentum drops dramatically. However, the fact that the secular slow-down of PNS rotation is much less pronounced in model Alcar-Rw (Fig. 9) suggests that numerical diffusion of angular momentum only plays a secondary role. Instead, the slow-down of core rotation in model Alcar-Rs appears to be connected to the action of magnetic stresses that initially build up due to field winding and then become strong enough to make the inner core spin more slowly (akin to the winding and unwinding of a spring). This is reflected by the decreasing toroidal field strength in the innermost zones in model Alcar-Rs at late times (Fig. 6, bottom panel), and the process of winding and unwinding also contributes to the long-period oscillation in the PNS magnetic energy in this model (Fig. 4, bottom panel). In the angular momentum profile at |$0.05\, \mathrm{s}$|⁠, we find an interesting trace of such long-term torsional motions in the form of a dip in the specific angular momentum outside a radius of |$10\, \mathrm{km}$| (Fig. 8). The long-period oscillations due to the winding and unwinding of the field are absent in the CoCoNuT-FMT models because the inner core is forced to rotate uniformly.

The angular momentum profiles of the weak-field models Alcar-Rw and coco-HLLC-Rw are much more similar (Fig. 9). Both models develop humps in their angular momentum profiles with sharp angular momentum gradients at the outer flank of the hump. The humps are, however, located at different radii, which reflects the aforementioned tendency towards a more inflated PNS surface structure in Alcar (cp. Fig. 5). The evolution of the angular momentum profile in the PNS core is quite similar despite different numerical treatments for this region. Different from the strong-field model, magnetic stresses remain too weak to substantially affect the rotational evolution of the PNS core over the time-scales of the simulation. Model coco-HLLC-Rw has higher specific angular momentum at radii less than about |$3\, \mathrm{km}$|⁠, possibly avoiding a slight loss of angular momentum into the central zone in Alcar-Rw. On the other hand, it does not capture the moderate amount of differential rotation in the PNS core due to its core treatment, while Alcar-Rw maintains the profile of differential rotation quite stably after an initial drop of specific angular momentum in the innermost zones in the first |$50\, \mathrm{ms}$| after bounce.

4.4 Neutrino emission

Even in a magnetorotational supernova, neutrinos still play a major role by influencing the nucleosynthesis conditions in the outflows (e.g. Janka 2012), and in some cases also by providing subsidiary energy input into the outflows. Since the fast-multi group transport scheme used in CoCoNuT-FMT (Müller & Janka 2015) and the two-moment method implemented in Aenus-Alcar (Just et al. 2015) represent quite different approaches to multigroup neutrino transport, and since the codes do not use exactly the same set of neutrino interaction rates, differences in the neutrino emission between the two codes are expected. However, the impact of the transport treatment on the neutrino emission can hardly be disentangled from the effect of the PNS structure, which is sensitive to the multidimensional dynamics of the models as outlined in the previous section. In contrast to 1D code comparisons, it would be a futile undertaking to fully track down the cause for all the differences in the neutrino emission. We therefore compare the neutrino emission in the two codes at a rather descriptive level, only occasionally highlighting phenomena that can be traced either to differences in the neutrino transport or to the model dynamics with greater confidence.

Neutrino luminosities for the various Alcar models and the standard-resolution CoCoNuT model are shown in Fig. 10 and reveal substantial differences. Neutrino mean energies are shown in Fig. 11.

Figure 10.

Neutrino luminosities of electron neutrinos (ν_e, top row), electron antineutrinos (⁠|$\bar{\nu }_\mathrm{e}$|⁠, middle), and heavy-flavour neutrinos (ν_x, bottom), measured at a radius of |$2000\, \mathrm{km}$|⁠.

Open in new tabDownload slide

Figure 11.

Neutrino mean energies of electron neutrinos (ν_e, top), electron antineutrinos (⁠|$\bar{\nu }_\mathrm{e}$|⁠, middle), and heavy-flavour neutrinos (ν_x, bottom), measured at a radius of |$2000\, \mathrm{km}$|⁠.

Open in new tabDownload slide

Especially at early times, the neutrino luminosities of all flavours tend to be higher in the CoCoNuT models, while the corresponding CoCoNuT models with the HLLC and HLLE solver show very similar luminosities. The higher early luminosities in CoCoNuT reflect differences in the collapse dynamics and early post-bounce accretion rate (Fig. 2) that were discussed in Section 4.1.

The temporal evolution of the neutrino mean energies is qualitatively similar, but CoCoNuT models are characterized by somewhat higher mean energies throughout, except for the heavy-flavour neutrinos.

In the case of the strong-field models, the electron neutrino and antineutrino luminosity is higher by up to |$2\times 10^{52}\, \mathrm{erg} \, \mathrm{s}^{-1}$| in models coco-HLLC-Rs and coco-HLLE-Rs compared to Alcar-Rs during the first |$\mathord {\sim }100\, \mathrm{ms}$| after bounce and then gradually become more similar to the Alcar model. The higher electron-flavour mean energies are consistent with smaller PNS radii in CoCoNuT; the combination of smaller radii and higher surface temperatures fortuitously leads to similar neutrino luminosity as in Alcar after |$\mathord {\sim }200\, \mathrm{ms}$|⁠. Note that unlike the electron-flavour luminosity, which is affected by differences in the early post-bounce accretion rate, the heavy-flavour luminosity remains consistently higher in the strong-field CoCoNuT models throughout the simulation. The faster contraction of the PNS in CoCoNuT alone cannot account for the stronger emission since the heavy-flavour neutrino mean energy is not substantially higher in CoCoNuT. Different interaction rates likely play a substantial role here. The omission of electron–positron annihilation in CoCoNuT can theoretically push the luminosity and mean energy up or down depending on the gradient of the temperature and chemical potentials near the neutrinosphere, the effective treatment of neutrino–nucleon scattering will tend to push the neutrino luminosity and mean energy down, and the strangeness corrections and nucleon correlations will likely push them upwards.

It is noteworthy that there is a marked divergence between models coco-HLLC-Rs and coco-HLLE-Rs after |$\mathord {\sim }0.4 \, \mathrm{s}$|⁠, which is when coco-HLLC-Rs develops a more fluffy PNS surface structure than coco-HLLE-Rs (Fig. 4). At this point, the electron-flavour mean energies in model coco-HLLC-Rs drop below those in coco-HLLE-Rs by about |$1\, \mathrm{MeV}$|⁠. Model coco-HLLC-Rs also exhibits a noticeably smaller heavy-flavour luminosity than coco-HLLE-Rs from |$0.3 \, \mathrm{s}$| onward. While the variations between these two models are considerably smaller than the differences between the CoCoNuT and Alcar models, this illustrates that neutrino luminosities and mean energies show considerable sensitivity to subtle dynamical variations in multidimensional MHD supernova models of rapidly rotating progenitors.

Among the weak-field models, we see similar differences in the early neutrino emission that can be ascribed to the same effects as for the strong-field models. In contrast to the strong-field models, there is, however, a cross-over in the electron flavour luminosities, with a phase of higher luminosities in Alcar-Rw compared to coco-HLLC-Rw and coco-HLLE-Rw. This is again consistent with the mass accretion rates in the two codes. At late times, the CoCoNuT models again exhibit higher electron-flavour luminosities; at this junction the more compact PNS surface structure is likely part of the explanation. Interestingly, the heavy-flavour luminosities and mean energies become similar between Alcar and CoCoNuT at late times. As heavy-flavour neutrinos originate from higher densities, these are less affected by the different ‘fluffiness’ of the PNS surface structure in the two codes.

The above results suggest that agreement between different codes in the neutrino luminosities and mean energies can only be achieved at a level of |$10{-}20{{\ \rm per\ cent}}$| in MHD models of rapidly rotating supernova progenitors. This is decidedly worse than what was achieved in the comparison between the fmt transport scheme and the Vertex neutrino transport code in Müller & Janka (2015) in the non-rotating case. The interplay between the multidimensional dynamics and the neutrino emission proves a considerable obstacle in analysing and eliminating code-dependent differences in the neutrino emission.

4.5 Nucleosynthesis conditions

Because of the ejection of neutron-rich material of relatively high entropy in the MHD-powered jets, magnetorotational supernovae have been proposed as an alternative site for the rapid neutron capture process (r-process, e.g. Winteler et al. 2012; Halevi & Mösta 2018; Grimmett et al. 2021; Reichert et al. 2021). This is supported by observational indications for r-process enrichment in low-metallicity environments that may be difficult to explain by neutron stars mergers (Nishimura, Takiwaki & Thielemann 2015; Kobayashi, Karakas & Lugaro 2020). In this context, it is important to analyse whether the different neutrino emission in Alcar and CoCoNuT has a substantial impact on the nucleosynthesis conditions in the outflows.

To this end, we visualise the distribution of the electron fraction Y_e in the ejecta in the different models using histograms for all the weak- and strong-field models (Fig. 12). The unbound ejecta are taken to comprise material with positive total specific energy, positive radial velocity, and a density lower than |$10^{11}\, \mathrm{g}\, \mathrm{cm}^{-3}$| in order not to accidentally include material in the PNS. Histograms are plotted at a post-bounce time of |$\mathrm{\approx } 0.4 \, \mathrm{s}$| for coco-HLLC-Rs, coco-HLLE-Rs, and Alcar-Rs, and at |$\mathrm{\approx } 0.85\, \mathrm{s}$| for coco-HLLC-Rw, coco-HLLE-Rw, and Alcar-Rw.

Figure 12.

Histograms of the electron fraction Y_e of the ejecta for the strong-field models at |$\mathord {\approx } 0.40\, \mathrm{s}$| and at |$\mathord {\approx } 0.85\, \mathrm{s}$| for the weak-field models. The four panels show, from top to bottom, the Y_e-distribution in the coco-HLLC, coco-HLLE, coco-Res and Alcar models, respectively.

Open in new tabDownload slide

All the strong-field models show a neutron-rich tail in the ejecta distribution that extends down to values of Y_e as low as |$0.25{-}0.3$|⁠, as well as a smaller amount of proton-rich ejecta. Neutron-rich conditions arise in matter that is quickly ejected in the magnetically collimated jets from close to the PNS so that Y_e freezes out at low values before it can be raised by neutrino irradiation (Wanajo, Janka & Müller 2011; Müller 2016). In slowly moving ejecta, neutrino absorption tends to establish proton-rich conditions because of similar ν_e and |$\bar{\nu }_e$| luminosities and small differences between the ν_e and |$\bar{\nu }_e$| mean energies (Fröhlich et al. 2006; Pruet et al. 2006). In all the weak-field models, which do not develop fast jets, the ejecta are pre-dominantly proton-rich, similar to recent multi-D models of neutrino-driven explosion using multigroup neutrino transport (Eichler et al. 2018; Wanajo et al. 2018; Burrows et al. 2020), although all models contain small amounts of neutron-rich ejecta except for Alcar-Rw.

Despite the differences in the neutrino-emission between CoCoNuT and Alcar, the qualitative shape of the Y_e-distribution is quite similar in all models. This is not as surprising as it might first seem, since Y_e is less sensitive to the absolute values of the neutrino luminosities and mean energies than to the differences between electron neutrinos and antineutrinos, which are similar for both codes.

Nonetheless, we still find noticeable quantitative differences in the Y_e-distribution. Once more, disentangling how the numerical methods and approximations employed in the two codes shape the distribution proves difficult, but a few plausible effects can be discerned.

Among the strong-field models, coco-HLLE-Rs has the least amount of material with Y_e < 0.4. This is consistent with the more diffusive nature of the HLLE solver, which will tend to raise Y_e in the jets by numerical mixing with less neutron-rich material. However, the lowest Y_e in model coco-HLLE-Rs is even lower than in coco-HLLC-Rs, even though the total amount of material below Y_e < 0.4 is substantially smaller. This is not a contradiction because the HLLE solver sometimes produces more violent, unsteady flow with bigger flow structures in the subsonic regime (Müller 2009). The amount of strongly neutron-rich material in Alcar-Rs is higher than in the CoCoNuT models and Y_e reaches lower minimum values (0.23 in Alcar-Rs as opposed to 0.28 in coco-HLLC-Rs). However, it is not possible to disentangle whether this is due to less severe numerical dissipation, different jet dynamics, or the neutrino transport treatment. The neutrino transport likely plays some role because model Alcar-Rs shows a rather remarkable proton-rich tail up to Y_e = 0.8, whereas the strong-field CoCoNuT models only eject moderately proton-rich material. It is unlikely that numerical diffusion alone could compress the distribution of proton-rich material so significantly as to explain the narrow distribution in the CoCoNuT model on its own.

The Y_e-distributions in the weak-field models agree inasmuch as the bulk of the ejecta fall into the moderately proton-rich range of |$Y_\mathrm{e}=0.5{-}0.6$|⁠. However, one notices a more extended proton-rich tail in Alcar-Rw compared to the CoCoNuT-FMT models. Unlike model Alcar-Rw, the CoCoNuT models also have a small admixture of slightly neutron-rich ejecta. Interestingly, the Y_e-distribution in model coco-HLLE-Rw is broader than in coco-HLLC-Rw and shifted to slightly smaller values. From the viewpoint of nucleosynthesis, the impact of the Riemann solver appears quite substantial. A shift of the minimum Y_e from 0.48 in coco-HLLC-Rw to about 0.39 in coco-HLLE-Rw will imply a considerable change in iron-group nucleosynthesis from NSE freeze-out from a predominant contribution of ⁵⁶Ni to a substantial contribution of neutron-rich nuclides like |$^{60}\, \mathrm{Fe}$| and |$^{48}\, \mathrm{Ca}$| (Hartmann, Woosley & El Eid 1985; Wanajo et al. 2018) and others. In summary, our code comparison so far already suggests that yields from the tail of the Y_e-distribution in MHD supernova models are expected to be quite sensitive to numerical details, whereas the composition of the bulk of the ejecta close to Y_e = 0.5 appears more robust. This suggests that the models can be more confidently used to predict ⁵⁶Ni yields that are relevant for supernova light curves, as ⁵⁶Ni is produced in the regime close to and above Y_e = 0.5, whereas r-process yields from simulations may be prone to considerable uncertainties.

5 RESULTS – RESOLUTION STUDY

To supplement the code comparison in Section 3, we compare models coco-HLLC-Rs and coco-HLLC-Rw to the corresponding high-resolution models coco-Res-Rs and coco-Res-Rw, which have been calculated with an angular resolution of 0.7^○ (256 zones) instead of 1.4^○ (128 zones). By revisiting parts of the analysis in Section 4 we shall not only explore a further dimension of uncertainties in MHD supernova models, but also obtain clues about the reasons underlying the differences between Alcar and CoCoNuT models.

5.1 Shock propagation and explosion energy

Figs 13 and 14 show the shock trajectories and explosion energies of the standard- and high-resolution models as the simplest global metrics of the dynamics. Among the strong-field models, we note that the high-resolution model coco-Res-Rs exhibits somewhat faster shock expansion than the baseline model coco-HLLC-Rs. The explosion energy is lower in coco-Res-Rs, however, and grows only at a rate of |$\mathord {\approx } 1.7 \times 10^{51}\, \mathrm{erg}\, \mathrm{s}^{-1}$| as opposed to |$\mathord {\approx } 2 \times 10^{51}\, \mathrm{erg},\mathrm{s}^{-1}$| in coco-HLLC-Rs. This suggests that uncertainties of order |$\mathord {\sim }10{{\ \rm per\ cent}}$| in these global quantities are expected due to limitations in resolution, or possibly from stochastic model variations (cp. Cardall & Budiardja 2015) alone. We note that the small explosion energy of the high-resolution model is not at odds with the resolution study of Sawai & Yamada (2016), who found a faster growth of the explosion energy for higher resolution in MHD-driven supernova models, but considered a somewhat different scenario. Sawai & Yamada (2016) started their simulations with smaller seed fields and attempted to follow field amplification by the MRI using very high resolution inside the PNS, whereas our simulations start with sufficiently strong fields to promptly trigger a magnetorotational explosion.

Figure 13.

Evolution of the maximum shock radius for the coco-HLLC models at standard resolution compared to the corresponding high-resolution coco-Res models. Both the weak- and strong-field case are shown. field cases. The radius given in a log10 scale.

Open in new tabDownload slide

Figure 14.

Comparison of the evolution of the diagnostic explosion energy E_expl in the strong- and weak-field models coco-HLLC at standard resolution to the corresponding high-resolution coco-Res models.

Open in new tabDownload slide

For the weak-field models, we observe small differences in the shock propagation and explosion energy between models coco-Res-Rw and coco-HLLC-Rw as well. Towards the end of the simulations, coco-Res-Rw has a larger maximum shock radius, but the difference in shock propagation is quite variable over time. Similarly, coco-Res-Rw has a lower explosion energy until |$0.9\, \mathrm{s}$| but then overtakes model coco-HLLC-Rw. The variations between these two models may thus be purely stochastic in nature.

5.2 Proto-neutron star structure and neutrino emission

Fig. 15 shows the PNS mass, angular momentum, and magnetic energy for models coco-Res-Rs, coco-Res-Rw, coco-HLLC-Rs, and coco-HLLC-Rw. In the weak-field case, the high-resolution and standard-resolution runs show almost perfect agreement in these quantities. There is only a hint of a slightly faster rise of the PNS magnetic energy in coco-HLLC-Rw close to the end of the simulation. However, more substantial differences emerge between the strong-field models coco-HLLC-Rs and coco-Res-Rs. After about |$0.3 \, \mathrm{s}$|⁠, the PNS angular momentum in model coco-Res-Rs starts to decrease faster than in coco-HLLC-Rs and ends up with a |$5{{\ \rm per\ cent}}$| smaller value. More intriguingly, the PNS magnetic energy in model coco-Res-Rs peaks at almost double the value in coco-HLLC-Rs. The evolution of the PNS radius also differs significantly (Fig. 16). Model coco-Res-Rs shows significant radius inflation at about |$0.3\, \mathrm{s}$|⁠, and while a similar phenomenon occurs in the standard-resolution run, this happens about |$0.1 \, \mathrm{s}$| later. This is noteworthy because a similar divergence of the PNS radii was found in our comparison between the strong-field Alcar and CoCoNuT models in Section 4. This is a further indication that the PNS surface region is subject to some form of tipping process that results in the adjustment to a more stable structure under conditions that are highly sensitive to the numerical implementation and grid resolution. Models coco-HLLC-Rs and coco-Res-Rs offer an opportunity to visualize the compact and inflated state of the PNS more clearly.

Figure 15.

Comparison of the PNS mass M_PNS, angular momentum J_PNS, and magnetic energy E_mag, as a function of post-bounce time for the standard-resolution coco-HLLC models and the high-resolution coco-Res models. The PNS is defined as the region where |$\rho \gt 10^{11} \, \mathrm{g}\, \mathrm{cm}^{-3}$|⁠.

Open in new tabDownload slide

Figure 16.

Comparison of the time evolution of the angle-averaged PNS radius R_PNS between CoCoNuT-FMT models with different resolution.

Open in new tabDownload slide

In Fig. 18 (left-hand panel), we show 2D plots of the entropy in the innermost |$100\, \mathrm{km}$| of the grid in models coco-Res-Rs and coco-HLLC-Rs. In both cases, the PNS has a strongly oblate structure, and its surface even forms something of a torus or thick disc. In the high-resolution model coco-Res-Rs, the disc-like structure is considerably more extended than in model coco-HLLC-Rs. We note that the edge of the PNS accretion disc at about |$50\, \mathrm{km}$| in coco-HLLC-Rs corresponds to the steep gradient in angular momentum profile at late times in this model (Fig. 17, cp. also Fig. 8). This suggests that the inflation of the disc is caused by outward angular momentum transport from the inner regions of the PNS by strong magnetic fields to larger radii. This is consistent with our interpretation that the different angular momentum profiles in the Alcar and CoCoNuT discussed in Section 4 are due a bifurcation instability that leads to a different PNS surface structure when triggered.

Figure 17.

Spherically averaged specific angular momentum as a function of radius at different post-bounce times for the strong-field models coco-Res-Rs (dotted lines) and coco-HLLC-Rs (solid lines).

Open in new tabDownload slide

We can see a similar effect, albeit at smaller level, in the comparison between models coco-HLLE-Rs and coco-HLLC-Rs with the latter exhibiting a more fluffy disc structure than the HLLE model (Fig. 18, right-hand panel). Intriguingly, the more dissipative Riemann solver produces a disc structure that is less inflated than with the more accurate HLLC solver. Less numerical dissipation on the grid scale – either due to the use of a more accurate solver or higher resolution – seems to enhance rather than to inhibit angular momentum transport. This appears counterintuitive at first glance, but there are other known instances of instabilities being suppressed by insufficient numerical dissipation, e.g. spurious behaviour of the Kelvin–Helmholtz instability in smoothed-particle hydrodynamics can be cured by introducing appropriate artificial conductivity at contact discontinuities (Price 2008). On the other hand, the instability could still arise more easily at higher resolution due to increased small-scale numerical noise. Alternatively, the instability might be triggered by stronger magnetic field amplification in model coco-Res-Rs (Fig. 15, bottom panel).

Figure 18.

Left: 2D plots of the entropy s in the vicinity of the PNS for the standard-resolution model coco-HLLC-Rs (left half) and the high-resolution model coco-Res-Rs (right half) at |$\mathord {\approx } 0.46\, \mathrm{s}$|⁠. Right: Corresponding plot of the entropy s in model coco-HLLE-Rs (left half) and coco-HLLC-Rs (right half) at |$\mathord {\approx } 0.45\, \mathrm{s}$| The entropy in both plots is capped at 12k_b/baryon to make the low-entropy, disc-like PNS surface region more clearly visible.

Open in new tabDownload slide

Regardless of the reason, the divergent disc structures in models coco-HLLC-Rs, coco-Res-Rs, and coco-HLLE-Rs further illustrate that the dynamics of MHD supernova models is highly sensitive to details of the numerical implementation such as the Riemann solver, and to the grid resolution, even though some bulk metrics like the explosion energy do not immediately show the same strong dependency.

Unsurprisingly, the variations in PNS surface structure among the strong-field CoCoNuT models are also reflected in the neutrino emission. Neutrino luminosities and mean energies for all neutrino species for the weak- and strong-field HLLC models at different resolution are shown in Figs 19 and 20, respectively. The neutrino luminosities hardly show any resolution dependence, and neither do the mean energies for the weak-field models. However, the onset of disc inflation in model coco-Res-Rs translates into a marked divergence of the neutrino energies of all flavours from model coco-HLLC-Rs. At the end of the simulation, the mean energies of the electron flavour neutrinos are about |$2\, \mathrm{MeV}$| lower, which is comparable to the differences between Alcar and CoCoNuT, and consistent with the significant expansion of the PNS radius in coco-Res-Rs.

Figure 19.

Comparison of the luminosities of electron neutrinos (ν_e, top row), electron antineutrinos (⁠|$\bar{\nu }_\mathrm{e}$|⁠, middle) and heavy-flavour neutrinos (ν_x, bottom) for the standard-resolution coco-HLLC models and the high-resolution coco-Res models. Luminosities are measured as spherical averages at a radius of |$2000\, \mathrm{km}$|⁠.

Open in new tabDownload slide

Figure 20.

Comparison of the mean energies of electron neutrinos (ν_e, top row), electron antineutrinos (⁠|$\bar{\nu }_\mathrm{e}$|⁠, middle) and heavy-flavour neutrinos (ν_x, bottom) for the standard-resolution coco-HLLC models and the high-resolution coco-Res models, measured as spherical averages at a radius of |$2000\, \mathrm{km}$|⁠.

Open in new tabDownload slide

5.3 Nucleosynthesis conditions

Finally, we show histograms of the Y_e-distribution of the ejecta in model coco-Res-Rs and coco-Res-Rw in Fig. 12. As in Section 4.5, the distribution is computed for the unbound outflowing material at |$\mathrm{\approx } 0.40\, \mathrm{s}$| for the strong-field case and at |$\mathord {\approx } 0.85\, \mathrm{s}$| for the weak-field case.

We notice much more extended tails in the high-resolution models. In model coco-Res-Rs, the neutron-rich tail reaches down about as far as in the Alcar model to 0.21, even though the amount of very neutron-rich material is still smaller than in Alcar. In model coco-Res-Rw, the proton-rich tail is now also quite flat as in model Alcar-Rw, only that a numerical cut-off is applied at Y_e = 0.7, corresponding to the upper boundary of the low-density NSE table that is used above temperatures of |$4.5\, \mathrm{GK}$|⁠. However, the ejecta in model coco-Res-Rw still cluster more strongly at a slightly proton-rich value of 0.7 than in Alcar-Rw.

It is noteworthy that the impact of the grid resolution and the Riemann solver on the Y_e-distribution is of a similar magnitude as the differences between the CoCoNuT and Alcar models, i.e. between codes with independent and rather dissimilar neutrino transport schemes. While the impact on actual yields would have to be evaluated by detailed post-processing of the nucleosynthesis, the sensitivity of the Y_e-distribution suggests that r-process yields from current magnetorotational supernova models may still be subject to substantial quantitative uncertainties as these hinge on the neutron-rich tail of the ejecta distribution.

6 CONCLUSIONS

We conducted the first code comparison of 2D MHD supernova models with the two independent codes CoCoNuT-FMT and Aenus-Alcar. Starting from the same initial conditions, we simulated the collapse and explosion of a rapidly rotating helium star progenitor model with an initial helium star mass of 35M_⊙ (Woosley & Heger 2006) for two initial magnetic field configurations with the same dipole geometry but different initial field strengths. The initial poloidal and toroidal field strengths are chosen as |$\mathrm{B_\mathrm{pol,tor} = 10^{12}G}$| in the strong-field case and as |$\mathrm{B_\mathrm{pol,tor} = 10^{10}G}$| in the weak-field case. For the CoCoNuT-FMT code, we also investigated the impact of different Riemann solvers (HLLC and HLLE) and the grid resolution CoCoNuT-FMT in a minimal resolution study.

We find overall qualitative agreement between the codes. Both in the strong-field case and in the weak-field case CoCoNuT and Alcar yield explosions at similar times with broadly similar shock trajectories. None the less, quantitative differences emerge upon closer inspection. Even within the same code, there is an appreciable sensitivity to numerical details like the Riemann solver and the grid resolution, although the various CoCoNuT models are still more similar to each other than to the corresponding Alcar run.

Notable differences include a more rapid explosion in the strong-field Alcar model compared to CoCoNuT-FMT, which translates into a significantly higher explosion energy. At the end of the simulations, the explosion energies still differ by tens of per cent, even though the relative difference may shrink on longer time-scales. Similarly, large differences in explosion energy are found for the weak-field models, although these may be impacted by stochastic model variations. Closer inspection reveals differences that can be traced to, or plausibly be associated with the disparate numerical schemes used in the two codes. The treatment of the innermost |$10\, \mathrm{km}$| in CoCoNuT-FMT as a quasi-spherical, uniformly rotating region leads to systematic differences in field amplification in the PNS. The two codes also show systematic differences in the PNS surface structure, with Alcar models tending to have a more ‘fluffy’ PNS surface or accretion disc around the PNS core. The specific angular momentum profiles and the field amplification in the PNS surface region are also significantly different. The different PNS structure is reflected in the neutrino emission and explains the tendency towards smaller electron-flavour neutrino mean energies in Alcar. However, difference in the neutrino emission also arise due to subtle differences in the infall dynamics and hence the mass accretion rate, and due to the disparate neutrino transport methodology.

We also compared the composition of the ejecta in the various models, specifically the distribution of the electron fraction Y_e as a crucial factor that decides, e.g. about the production of r-process material in MHD-driven supernovae. The CoCoNuT and Alcar models are qualitatively in agreement in that they produce considerable neutron-rich material in the strong-field case and pre-dominantly proton-rich material in the weak-field case. However, the width of the Y_e-distribution and the amount of very neutron-rich or very proton-rich material varies noticeably between CoCoNuT and Alcar, and is also quite sensitive to the numerical resolution and the Riemann solver.

In contrast to 1D comparisons of supernova codes (Liebendörfer et al. 2005; Müller et al. 2010; O’Connor et al. 2018) and multi-D code comparisons for non-rotating, non-magnetized progenitors (Müller & Janka 2015; Cabezón et al. 2018; Just et al. 2018), disentangling and eliminating the causal factors behind discrepant results proves much more difficult. Therefore, no attempt has been made to achieve perfect agreement between the two codes. It is perhaps rather astonishing that despite all the differences, the most coarse-grain metrics of the dynamics like the shock propagation and the explosion energy, and the Y_e-distribution still end up in qualitative agreement.

Our comparison study suggests that MHD supernova models in axisymmetry are quite sensitive to details of the numerical implementation, resolution, and the employed microphysics, and therefore still subject to significant uncertainties. While they appear robust on a qualitative level, there is clearly not the same level of convergence between codes that has been reached in 1D supernova simulations (Liebendörfer et al. 2005; Müller et al. 2010; O’Connor et al. 2018).

Of necessity, our work is only a first step in bracketing and reducing the uncertainties in MHD supernova models. Follow-up studies are highly desirable. Comparison studies with reduced model complexity, e.g. without neutrino transport in the post-bounce phase and an identical treatment of the collapse phase, may be useful to specifically identify and iron out differences between MHD solvers, rather than dealing with a coupled radiation MHD problem. Focusing on MHD only would also make it easier to conduct a comparison in 3D rather than in 2D, which will be critical to study aspects like dynamo field amplification and the kink instability in jet-driven explosions (Mösta et al. 2014; Kuroda et al. 2020). On the other hand, the impact of uncertainties in MHD supernova models with neutrino transport should also be explored in more depth. In particular, it would be desirable to investigate the implications for detailed nucleosynthesis yields rather than just considering the Y_e-distributions in our current study. While MHD supernova simulations still need to further mature in many respects, e.g. by considering more realistic initial magnetic field configurations (Varma & Müller 2021), further code comparisons will be an important tool on the long road towards a robust understanding of the role of magnetic fields in the explosions of massive stars.

ACKNOWLEDGEMENTS

We thank A. Heger, D. Price, and S. Campbell for useful discussions. BM acknowledges support by ARC Future Fellowship FT160100035. MO acknowledges support from the Spanish Ministry of Science, Education and Universities (PGC2018-095984-B-I00) and the Valencian Community (PROMETEU/2019/071), the European Research Council under grant EUROPIUM-667912, and from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 279384907 – SFB 1245 as well as from the Spanish Ministry of Science via the Ramón y Cajal programme (RYC2018-024938-I). Part of this work has been performed using computer time allocations from Astronomy Australia Limited’s ASTAC scheme, the National Computational Merit Allocation Scheme (NCMAS), and an Australasian Leadership computing grant on the NCI NF supercomputer Gadi. This research was supported by resources provided by the Pawsey Supercomputing Centre, with funding from the Australian Government and the Government of Western Australia.

DATA AVAILABILITY

The data from our simulations will be made available upon reasonable requests made to the authors.

Footnotes

REFERENCES