# 3D Simulations of Strongly Magnetised Non-Rotating Supernovae: Explosion
Dynamics and Remnant Properties

Vishnu Varma1,2, Bernhard Müller1, Fabian R. N. Schneider3,4
1 School of Physics and Astronomy, 10 College Walk, Monash University,
Clayton, VIC 3800, Australia
2 Astrophysics Group, Lennard-Jones Laboratories, Keele University, Keele ST5
5BG, UK
3 Heidelberger Institut für Theoretische Studien, Schloß-Wolfsbrunnenweg 35,
69118 Heidelberg, Germany
4 Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität
Heidelberg, Mönchhofstraße 12–14, 69120 Heidelberg, Germany E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We investigate the impact of strong initial magnetic fields in core-collapse
supernovae of non-rotating progenitors by simulating the collapse and
explosion of a $16.9\,M_{\odot}$ star for a strong- and weak-field case
assuming a twisted-torus field with initial central field strengths of
$\mathord{\approx}10^{12}\,\mathrm{G}$ and
$\mathord{\approx}10^{6}\,\mathrm{G}$. The strong-field model has been set up
with a view to the fossil-field scenario for magnetar formation and emulates a
pre-collapse field configuration that may occur in massive stars formed by a
merger. This model undergoes shock revival already black $100\,\mathrm{ms}$
after bounce and reaches an explosion energy of $9.3\times
10^{50}\,\mathrm{erg}$ at $310\,\mathrm{ms}$, in contrast to a more delayed
and less energetic explosion in the weak-field model. The strong magnetic
fields help trigger a neutrino-driven explosion early on, which results in a
rapid rise and saturation of the explosion energy. Dynamically, the strong
initial field leads to a fast build-up of magnetic fields in the gain region
to 40% of kinetic equipartition and also creates sizable pre-shock ram
pressure perturbations that are known to be conducive to asymmetric shock
expansion. For the strong-field model, we find an extrapolated neutron star
kick of $\mathord{\approx}350\,\mathrm{km}\,\mathrm{s}^{-1}$, a spin period of
$\mathord{\approx}70\,\mathrm{ms}$, and no spin-kick alignment. The dipole
field strength of the proto-neutron star is $2\times 10^{14}\,\mathrm{G}$ by
the end of the simulation with a declining trend. Surprisingly, the surface
dipole field in the weak-field model is stronger, which argues against a
straightforward connection between pre-collapse fields and the birth magnetic
fields of neutron stars.

###### keywords:

stars: massive – stars: magnetic fields – supernovae: general

††pubyear: 2022††pagerange: 3D Simulations of Strongly Magnetised Non-Rotating
Supernovae: Explosion Dynamics and Remnant Properties–3D Simulations of
Strongly Magnetised Non-Rotating Supernovae: Explosion Dynamics and Remnant
Properties

##  1 Introduction

In the last few decades, there has been significant progress in our
understanding of the role of magnetic fields in core-collapse supernovae. This
is owed, in part, to rapid improvements in numerical modelling and computing
capabilities.

Since traditional stellar evolution models (Heger et al., 2005) have predicted
weak pre-collapse magnetic field strengths, most of the progress in
understanding magnetic effects in core-collapse supernovae have been focused
on the so-called magnetorotational mechanism (LeBlanc & Wilson, 1970;
Bisnovatyi-Kogan et al., 1976; Müller & Hillebrandt, 1979; Burrows et al.,
2007; Obergaulinger et al., 2009; Sawai et al., 2013; Masada et al., 2015;
Mösta et al., 2015; Bugli et al., 2020; Obergaulinger & Aloy, 2021). This
mechanism relies on a rapidly differentially rotating star to amplify the weak
magnetic fields, e.g., by the magnetorotational instability (MRI; Balbus &
Hawley 1991; Akiyama et al. 2003) or an $\alpha$-$\Omega$-dynamo in the proto-
neutron star (Duncan & Thompson, 1992; Thompson & Duncan, 1993; Raynaud et
al., 2020). Strong magnetic stresses can then drive the explosion, usually in
the form of collimated jets.

However, stellar evolution models also show magnetic braking causes the
majority of supernova progenitor cores to rotate slowly (Heger et al., 2005;
Meynet et al., 2011) as the magnetic torques transport angular momentum to the
envelope very efficiently. Even without magnetic braking, the winds of
Galactic stars will prevent rapid rotation at core collapse. These predictions
are bolstered by the white dwarf spin period distribution, which can only be
explained if similarly efficient angular momentum transport mechanisms like
magnetic torques operate in low-mass stars (Suijs et al., 2008). It is thus
expected that the magnetorotational mechanism likely only operates in rare
hyperenergetic supernovae explosions (known as “hypernovae”), which are
sometimes linked to long-duration gamma-ray bursts (GRBs) as well (Woosley &
Bloom, 2006). These scenarios require special evolutionary pathways to limit
the spin-down of the progenitor cores (e.g., Petrovic et al., 2005; Woosley &
Heger, 2006; Georgy et al., 2012; Georgy et al., 2009; Cantiello et al., 2007;
Aguilera-Dena et al., 2020).

The role of magnetic fields in supernovae outside the magnetorotational
mechanism has so far received less attention. Without strong differential
rotation, there is no obvious avenue for generating dynamically relevant
fields that could become the main driver of an explosion. Nonetheless, some
studies have explored the generation and impact of magnetic fields in non-
rotating core-collapse supernovae in recent years. Early studies by Endeve et
al. (2010); Endeve et al. (2012) looked to understand the possible role of
standing accretion shock instabilities (SASI) to amplify magnetic fields after
core-collapse. They found that while the SASI oscillations can amplify the
magnetic fields on the neutron star to magnetar strengths, this does not
affect global shock dynamics. Global studies in axisymmetric (2D) simulations
by Obergaulinger et al. (2014), Matsumoto et al. (2020) and Jardine et al.
(2022) with MHD and neutrino transport, find that even very strong magnetic
fields ($\mathord{\approx}10^{12}\,\mathrm{G}$), without rotation, only play
an auxiliary in neutrino-driven explosions. However, recently, Müller & Varma
(2020) have shown that in 3D, even very weak initial magnetic fields as low as
$10^{6}\,\mathrm{G}$ can play a subsidiary role in aiding the neutrino-driven
mechanism thanks to efficient amplification by a turbulent dynamo. A 3D study
by Matsumoto et al. (2022) also showed that in more strongly magnetised
models, the amplified magnetic field on large hot bubbles just behind a
stalled shock can lead to a faster and more energetic explosion.

Another possibility for dynamically relevant magnetic fields, is that they
might already be present in the progenitor. This is known as the fossil-field
hypothesis for the generation of magnetic fields in neutron stars and, in
particular, in magnetars. If the initial fields are strong enough, the
compression via core-collapse alone would be sufficient to amplify the
magnetic fields to magnetar strengths of $10^{14\texttt{-}15}\,\mathrm{G}$,
through magnetic flux conservation (Braithwaite & Spruit, 2004; Spruit, 2009;
Ferrario et al., 2016).

To explain magnetars by the fossil-field hypothesis, a remarkably strong field
already would need to be present in the progenitor core in a subset of
supernova progenitors. Based on a rough coincidence of the magnetar fraction
and the fraction of roughly $7\%$ of massive stars with unusually strong
surface fields (Donati & Landstreet, 2009; Grunhut et al., 2017; Schöller et
al., 2017), a connection is often made between strongly magnetised massive
stars on the main sequence and magnetars with the underlying assumption that
the interior field is likewise strong and not destroyed through the subsequent
evolutionary stages up to iron core collapse (see, e.g., Quentin & Tout, 2018,
for the evolution of magnetised, low-mass stars up to white-dwarf formation).
Connecting the origin of magnetar fields back to early evolutionary stages is
appealing because stellar mergers might provide one natural avenue for
producing a population of strongly magnetised stars roughly at the rate
required to explain magnetars (Schneider et al., 2019); other channels for
forming strongly magnetised OB stars are also conceivable (see, e.g., Mestel,
1999, 2001; Moss, 2001; Donati & Landstreet, 2009), and this work should be
viewed as being agnostic about the exact origin of the strong magnetic fields.

Rate arguments in favour of the fossil-field hypothesis are not uncontested,
however. Some population synthesis models indicate that the simple fossil-
field hypothesis based on the distribution of surface magnetic field strengths
among massive stars cannot simultaneously reproduce both pulsar and magnetar
populations (Ferrario & Wickramasinghe, 2006; Makarenko et al., 2021). It has
also been found that the magnetar birth rate is comparable to that of neutron
stars (and could be as high as the neutron star birth rate), which may leave
too few known stars with surface $B\,>\,10\,\mathrm{kG}$ to account for the
number of magnetars (Woods, 2008; Keane & Kramer, 2008) and raises doubts
whether fossil fields are the only means for generating magnetar-strength
fields.

With many open questions about the origin of neutron star magnetic fields and
recent theoretical results that suggest a broader role of magnetohydrodynamic
effects in supernova explosions, further research into the impact of magnetic
fields on supernova explosion dynamics and compact remnant birth properties is
called for. In this study, we present new simulations of a strongly
magnetised, non-rotating neutrino-driven supernova in 3D and compare this to a
weakly magnetised model to better understand how the strength of the magnetic
field aids the explosion mechanism even without any initial rotation. While 1D
stellar evolution models indicate only weak magnetic fields developing in
supernova progenitors, our decision to explore the strongly magnetised regime
is motivated by the fossil-field hypothesis and also by the 3D shell burning
simulation of Varma & Müller (2021), which suggests the shell convection can
generate rather strong magnetic fields of $10^{10\texttt{-}11}\,\mathrm{G}$ in
the oxygen shell.

Our paper is structured as follows: In Section 2, we describe the progenitor
model as well as the initial conditions implemented in our simulations, this
is followed by a description of the numerical methods of the code we used,
CoCoNuT-FMT in Section 3. The results of the simulations are presented in
Section 4 where we first look at the large-scale explosion dynamics, and
analyse the role of the magnetic fields in the explosion, followed by an
investigation of the resulting proto-neutron star (PNS) properties. We
summarise our results and discuss their implications in Section 5.

##  2 Progenitor Model and Initial Conditions

We simulate the core collapse of a strongly magnetised, slowly-rotating
massive star of $\mathord{16.9\,\mathrm{M_{\odot}}}$ that is the result of a
binary merger of $8\,\mathrm{M_{\odot}}$ and $9\,\mathrm{M_{\odot}}$ stars as
modelled in Schneider et al. (2019); Schneider et al. (2020) in an attempt to
black reproduce the magnetic star $\tau$ Sco. The merger was simulated using
the AREPO code, followed by 1D stellar evolution of the merger product in the
MESA code. We map this MESA progenitor model to the CoCoNuT-FMT supernova code
at the onset of core collapse when the maximum infall velocity has reached
$1000\,\mathrm{km}\,\mathrm{s}^{-1}$. black With a maximum angular velocity of
$2.2\times 10^{-8}\,\mathrm{rad}\,\mathrm{s}^{-1}$ the rotation is, in fact,
so slow that it can be neglected altogether, hence we do not map the angular
velocity from the MESA progenitor model to CoCoNuT-FMT, and instead simulate
the non-rotating case.

The MESA model is evolved until core collapse with the same setup as in
Schneider et al. (2021). The entropy and density profiles mapped are shown in
Figure 1.

Three-dimensional progenitor models of stars with strong fossil fields in the
core are presently unavailable and unrealistic because they would require
multi-dimensional stellar evolution calculations over secular time scales. We
therefore manually specify a field configuration that is compatible with the
expectations for a strong fossil field. Prior to mapping the merger progenitor
into MESA, the magnetic field geometry of the star from Schneider et al.
(2019) is $80\texttt{-}85\%$ toroidal, but with a sizable poloidal component.
This is consistent with the stable twisted-torus structure found in radiative
zones found by Braithwaite & Nordlund (2006).

We hence initialise the magnetic field geometry using an analytical
prescription from Kamchatnov (1982) for such twisted-torus fields. Although
Kamchatnov (1982) did not prove the stability of this configuration, this is
not essential for our purpose. The fields have no time to adjust before core
collapse after the beginning of the simulation, and during the collapse even
an initially stable magnetic field would become unstable anyway. The key point
is that the configuration exhibits substantial magnetic helicity as one would
expect from any realistic strong fossil field.

In spherical polar coordinates, the components $A_{r}$, $A_{\theta}$, and
$A_{\varphi}$ of the vector potential for the initial field configuration are
given by

$\displaystyle
A_{r}=B_{0}\frac{(1+\mathcal{R}^{2})\cos\theta}{4(1+2\mathcal{R}^{2}+\mathcal{R}^{4})},$
(1) $\displaystyle
A_{\theta}=B_{0}\frac{(-1+\mathcal{R}^{2})\sin\theta}{4(1+2\mathcal{R}^{2}+\mathcal{R}^{4})},$
(2) $\displaystyle
A_{\phi}=B_{0}\frac{\mathcal{R}\sin\theta}{2(1+2\mathcal{R}^{2}+\mathcal{R}^{4})},$
(3)

where $B_{0}$ is the central (maximum) magnetic field strength, and
$\mathcal{R}=r/R_{0}$ is a rescaled radial coordinate.

As we are investigating the role of magnetic fields in aiding the explosion
mechanism for non-rotating supernova progenitors, we run two simulations using
the same progenitor with different initial magnetic field strengths, namely a
weak-field model, with $B_{0}=10^{6}\,\mathrm{G}$ and a strong-field model
with $B_{0}=10^{12}\,\mathrm{G}$ and $R_{0}=5\times 10^{8}\,\mathrm{cm}$.

The strong-field model serves as the representative case for the fossil-field
magnetar progenitor scenario, whereas the weak-field model provides a control
case. Since the magnetic field prescription from Equations (1–3) leads to a
very sharp drop-off in magnetic field strength with radius, we scale the inner
core to have magnetic field strengths of $B_{0}$, such that for the strong
magnetic field case, the pre-supernova oxygen shell has magnetic field
strengths of $B=10^{10}-10^{11}G$, (similar to Varma & Müller 2021). This is
similar to the MESA progenitor model which should have a core field strength
of $\mathord{\gtrsim}10^{11}\,\mathrm{G}$, given magnetic flux is conserved
from the original stellar merger calculation.

Figure 1: Inverse plasma-$\beta$ at the onset of collapse for the strong-field
model (top), and initial density and entropy profiles (bottom) for the
progenitor.

Without 3D progenitor models close to magnetohydrostatic equilibrium, a
delicate balance must be struck in the choice of the initial field strength.
The initial fields must not be too strong because this would render the
initial pre-collapse structure highly inconsistent, but to realistically mimic
a strong fossil field that would not have been overwhelmed by turbulent
convective motions during the late burning stages, the magnetic stresses ought
to be larger than the expected turbulent Reynolds stresses from convection.
The configuration chosen in this study places the ratio of magnetic to thermal
pressure (inverse plasma-$\beta$) in the oxygen shell at 0.05 as a reasonable
compromise (see Figure 1).

##  3 Numerical Methods

For our 3D simulations we employ the Newtonian magnetohydrodynamic (MHD)
version of the CoCoNuT-FMT code as described in detail in Müller & Varma
(2020); Varma & Müller (2021); Varma et al. (2021). The MHD equations are
solved in spherical polar coordinates using the HLLC (Harten-Lax-van Leer-
Contact) Riemann solver (Gurski, 2001; Miyoshi & Kusano, 2005) and piecewise
parabolic reconstruction (Colella & Woodward, 1984). The divergence-free
condition $\nabla\cdot\mathbf{B}=0$ is maintained using a modified version of
the original hyperbolic divergence cleaning scheme of Dedner et al. (2002) and
follows ideas of Tricco et al. (2016) to maintain energy conservation.
Compared to the original cleaning method, we rescale the Lagrange multiplier
$\psi$ to $\hat{\psi}=\psi/c_{\mathrm{h}}$, where $c_{\mathrm{h}}$ is the
hyperbolic cleaning speed. Details of this approach and differences to Tricco
et al. (2016) are discussed in the Appendix of Varma & Müller (2021). The
extended system of MHD equations for the density $\rho$, velocity
$\mathbf{v}$, magnetic field $\mathbf{B}$, the total energy density $\hat{e}$,
mass fractions $X_{i}$, and the rescaled Lagrange multiplier $\hat{\psi}$
reads,

$\displaystyle\partial_{t}\rho+\nabla\cdot\rho\mathbf{v}$ $\displaystyle=$
$\displaystyle 0,$ (4)
$\displaystyle\partial_{t}(\rho\mathbf{v})+\nabla\cdot\left(\rho\mathbf{v}\mathbf{v}-\frac{\mathbf{B}\mathbf{B}}{4\pi}+P_{\mathrm{t}}\mathcal{I}\right)$
$\displaystyle=$
$\displaystyle\rho\mathbf{g}-\frac{(\nabla\cdot\mathbf{B})\mathbf{B}}{4\pi},$
(5)
$\displaystyle\partial_{t}{\hat{e}}+\nabla\cdot\left[(e+P_{\mathrm{t}})\mathbf{v}-\frac{\mathbf{B}(\mathbf{v}\cdot\mathbf{B}){-c_{\mathrm{h}}\hat{\psi}\mathbf{B}}}{4\pi}\right]$
$\displaystyle=$
$\displaystyle\rho\mathbf{g}\cdot\mathbf{v}+\rho\dot{\epsilon}_{\mathrm{nuc}},$
(6)
$\displaystyle\partial_{t}\mathbf{B}+\nabla\cdot(\mathbf{v}\mathbf{B}-\mathbf{B}\mathbf{v})+\nabla(c_{\mathrm{h}}\hat{\psi})$
$\displaystyle=$ $\displaystyle 0,$ (7)
$\displaystyle\partial_{t}\hat{\psi}+c_{\mathrm{h}}\nabla\cdot\mathbf{B}$
$\displaystyle=$ $\displaystyle-\hat{\psi}/\tau,$ (8)
$\displaystyle\partial_{t}(\rho X_{i})+\nabla\cdot(\rho X_{i}\mathbf{v})$
$\displaystyle=$ $\displaystyle\rho\dot{X}_{i},$ (9)

where $\mathbf{g}$ is the gravitational acceleration, $P_{\mathrm{t}}$ is the
total pressure, $\mathcal{I}$ is the Kronecker tensor, $\tau$ is the damping
timescale for divergence cleaning, and $\dot{\epsilon}_{\mathrm{nuc}}$ and
$\dot{X}_{i}$ are energy and mass fraction source terms from nuclear
reactions. This system conserves the volume integral of a modified total
energy density $\hat{e}$, which also contains the cleaning field $\hat{\psi}$,

${\hat{e}}=\rho\left(\epsilon+\frac{v^{2}}{2}\right)+\frac{B^{2}+\hat{\psi}^{2}}{8\pi},$
(10)

where $\epsilon$ is the mass-specific internal energy.

To reduce numerical dissipation near the grid axis in our simulations, we have
modified the mesh coarsening algorithm of Müller et al. (2019a) by
implementing a third-order accurate slope-limited prolongation scheme.
Neutrinos are treated using the FMT (fast multi-group transport) scheme of
Müller & Janka (2015a), which solves the energy-dependent zeroth moment
equation for three neutrino species in the stationary approximation using a
closure obtained from a two-stream solution of the Boltzmann equation

The two models are run with a grid resolution of $550\times 128\times 256$
zones in $r$, $\theta$ and $\varphi$ (corresponding to $1.4^{\circ}$ in
angle), an exponential grid in energy spaces with $21$ zones from
$4\,\mathrm{MeV}$ to $240\,\mathrm{MeV}$. The simulations use the equation of
state of Lattimer & Swesty (1991) with an incompressibility modulus of
$K=220\,\mathrm{MeV}$ at high densities, and an ideal gas consisting of
photons, electrons, positions, and non-relativistic nucleons and nuclei, in
conjunction with an NSE solver above $5\,\mathrm{GK}$ and a flashing treatment
at lower temperatures (Rampp & Janka, 2000).

##  4 Results

###  4.1 Explosion Dynamics

Both the strong- and weak-field models collapse to a PNS within the same time,
$\mathord{\approx}0.1\,\mathrm{s}$ after the models are mapped to 3D. In
Figure 2, we see that while their initial shock trajectories are very similar,
they begin to deviate significantly after about $0.1\,\mathrm{s}$ post-bounce,
at which point the strong magnetic field model undergoes much more rapid shock
expansion. In the weak-field case, the shock is also revived, but shock
expansion sets in about $0.05\texttt{-}0.1\,\mathrm{s}$ later. Up to the time
of explosion in the strong-field model, the mass accretion rate agrees
perfectly in both runs (Figure 3), indicating that the strong initial field
does not spuriously affect the _spherically-averaged_ collapse dynamics. After
the onset of the explosion, the accretion rate drops visibly in both cases. In
the strong-field model, net accretion on the PNS black (defined as the radius
where the angle-averaged density reaches
$10^{11}\,\mathrm{g}\,\mathrm{cm}^{-3}$) ceases shortly after
$0.2\,\mathrm{s}$ post-bounce. The subsiding accretion is responsible for a
change in PNS surface structure that is reflected in the evolution of the PNS
radius and the gain radius in Figure 2. While the two models match very
closely in terms of PNS radius and gain radius for the first
$0.225\,\mathrm{s}$, they then start to diverge, with the PNS radius and gain
radius of the strong magnetic field models dropping more rapidly.

Shock revival also has consequences for the neutrino luminosities and mean
energies, which are shown in Figure 4. Up until $0.1\,\mathrm{s}$ after
bounce, the neutrino-emission is very similar in both models. black Small
differences arise because of stochastic variation in the early post-bounce
entropy profiles that are imprinted on the PNS by prompt convection and
slightly change the thermodynamics conditions at the neutrinosphere (cp.
Müller & Varma, 2020). Once the shock starts to expand in the strong-field
model, the electron neutrino and antineutrino luminosities drop significantly
faster than in the weak-field model, and the rise of their mean energies is
slower. Interestingly, the heavy-flavour neutrino mean energy remains higher
in the strong-field model, which is consistent with faster contraction of the
PNS. It is also noteworthy that there are small differences in the neutrino
emission between the strong- and weak-field models prior to shock revival,
which would actually suggest slightly better heating conditions in the weak-
field model. These can be ascribed to minute differences in the collapse phase
and the PNS structure at early post-bounce times that ultimately stem from the
slight perturbation of the spherical collapse dynamics by the introduction of
the initial fields.

Figure 2: Evolution of the maximum shock radius (solid), PNS radius (dotted),
and gain radius (dash-dotted) for the strong-field (black) and weak-field
(red) models.  Figure 3: Mass accretion rate $\dot{M}$ for the strong-field
model (black) and for the weak-field model (red), calculated at a radius of
$250\,\mathrm{km}$.  Figure 4: Neutrino luminosities (top) and mean energies
(bottom) of electron neutrinos ($\nu_{e}$), electron antineutrinos
($\bar{\nu}_{e}$) and heavy-flavour neutrinos ($\nu_{\mathrm{x}}$), measured
at a radius of $2000\,\mathrm{km}$.

Figure 5 shows the diagnostic explosion energy $E_{\mathrm{expl}}$ (Buras et
al., 2006) for both models. The diagnostic explosion energy is defined as an
integral over the region that is nominally unbound,

$E_{\mathrm{expl}}=\int\limits_{e_{\mathrm{tot}}>0}\rho
e_{\mathrm{tot}}\,\mathrm{d}V,$ (11)

where $e_{\mathrm{tot}}$ is the total energy density, i.e., the sum of the
internal, kinetic, gravitational, and magnetic energy density.

The diagnostic explosion energy of the strong-field model begins to rise
around the time when the shock trajectories of the two models deviate, and
reaches more tan $9\times 10^{50}\,\mathrm{erg}$ by the end of the simulation.
The explosion energy already appears to level off at this stage, which is
rather unusual when compared to current 3D hydrodynamic simulations of core-
collapse supernova explosions. Long-term 3D explosion models of massive
progenitors typically show a sustained rise of the explosion energy over up to
several seconds (Müller, 2015; Bruenn et al., 2016; Bollig et al., 2021), but
at a significantly smaller growth rate, more reminiscent of the weak-field
model, which has only reached about $10^{50}\,\mathrm{erg}$. The unusually
rapid growth and saturation of the explosion energy is further emphasised by
Figure black 6, which shows that the rate of growth
$\mathrm{d}E_{\mathrm{expl}}/\mathrm{d}t$ peaks at
$10^{52}\,\mathrm{erg}\,\mathrm{s}^{-1}$ about $0.22\,\mathrm{s}$ post-bounce,
and drops by more than a factor of five afterwards.

Figure 5: Evolution of the diagnostic explosion energy $E_{\mathrm{expl}}$ for
the strong-field model (black) and the weak-field model (red).  Figure 6:
Comparison of the growth rate ${\dot{E}_{\mathrm{expl}}}$ of the diagnostic
explosion energy (solid black and red curves) to the approximation for the
rate of energy ejection into the ejecta from Equation (13) (solid blue curve),
which accounts for the release of
$\epsilon_{\mathrm{rec}}=6\,\mathrm{MeV}/\mathrm{baryon}$ by nucleon
recombination for the measured mass outflow rate ${\dot{M}_{\mathrm{out}}}$
and the measured magnetohydrodynamic energy flux $F_{\mathrm{t}}$ at the base
of the gain region of the strongly magnetised model (blue). The nuclear
recombination power $\epsilon_{\mathrm{rec}}\dot{M}_{\mathrm{out}}$ is also
plotted on its own for both models (dash-dotted blue and red curves). The
volume-integrated neutrino heating rate $\dot{Q}_{\nu}$ is also shown (dotted
curves).  Figure 7: Radial profiles of the total magnetohydrodynamic energy
flux $F_{\mathrm{t}}$ with rest masses subtracted (black, solid), the purely
hydrodynamic total enthalpy flux $F_{\mathrm{h}}$ (black, dotted), and the
mass outflow rate $\dot{M}_{\mathrm{out}}$ (blue) for the strong-field model
at $0.20\,\mathrm{s}$ post-bounce. The energy and mass fluxes are computed as
integrals over the area covered by outflows with $v_{r}>0$. The vertical
dashed line at $\approx{\color{black}100}\,\mathrm{km}$ indicates where does
the angle-averaged total energy in the outflows become positive.

Strong initial fields thus evidently have a major impact on the explosion
dynamics, allowing the model to reach a typical supernova explosion energy
(Kasen & Woosley, 2009; Murphy et al., 2019) within only $0.3\,\mathrm{s}$
after bounce. As we shall show, this effect is largely of an indirect nature.
It comes about because the strong initial fields help trigger the explosion
early, rather than from direct delivery of the explosion energy by magnetic
stresses as in the magnetorotational paradigm. The explosion in energised in a
similar manner as in existing hydrodynamics models of neutrino-driven
explosions with a few subtle differences related to early shock revival.

Previous studies by Marek et al. (2006); Müller (2015); Bruenn et al. (2016)
have shown that for the neutrino-driven explosion mechanism, the growth of the
diagnostic explosion energy is largely determined by the mass outflow rate,
${\dot{M}_{\mathrm{out}}}$ of the ejecta, which are first lifted to
$e_{\mathrm{tot}}\approx 0$ by neutrino heating, before obtaining most of
their net positive energy from nucleon recombination. Ignoring secondary
effects like the scooping-up of bound material by the shock and losses of
explosion energy by turbulent braking of the expanding neutrino-heated
bubbles, this results in a growth rate of the explosion energy of
$\dot{E}_{\mathrm{expl}}\approx\epsilon_{\mathrm{rec}}\dot{M}_{\mathrm{out}}$,
where $\epsilon_{\mathrm{rec}}$ is the recombination energy per unit mass.

Comparing the rate of energy release
$\epsilon_{\mathrm{rec}}\dot{M}_{\mathrm{out}}$ by nucleon recombination to
the rate of growth of the explosion energy is rather tricky for the strong-
field model since $\dot{M}_{\mathrm{out}}$ is non-stationary and varies
strongly with radius (Figure 7), unlike in models with later shock revival
where $\dot{M}_{\mathrm{out}}\approx\mathrm{const.}$ over a large range in
radius (see, e.g., Figure 18 of Müller 2015). To capture the energy release by
nuclear recombination with reasonable accuracy during a rather non-stationary
outflow, it is critical to evaluate $\dot{M}_{\mathrm{out}}$ not too far away
from the actual recombination radius. In our subsequent analysis, we evaluate
$\dot{M}_{\mathrm{out}}$ at a radius of $400\,\mathrm{km}$ by integrating over
all directions where the radial velocity component $v_{r}$ is positive,

$\dot{M}_{\mathrm{out}}=\int\limits_{v_{r}>0}\rho
v_{r}r^{2}\,\mathrm{d}\Omega.$ (12)

Even in this case, however, nucleon recombination can only account for a
little more than half of the growth rate of the explosion energy.

The extra energy input into the explosion can be explained by $P\,\mathrm{d}V$
work (and to some extent, work by magnetic stresses) by the underlying layers
on the ejecta region. It is noteworthy that a more significant role of the
magnetohydrodynamic energy flux into the ejecta region compared to
recombination was also found for early explosions in axisymmetry (Bruenn et
al., 2016) and for some mass shells in rapidly developing electron-capture
supernova explosions (Janka et al., 2008).

Taking into account the hydrodynamic energy flux into the ejecta region we
expect the explosion energy to grow as

$\frac{\mathrm{d}E_{\mathrm{expl}}}{\mathrm{d}t}\approx
F_{\mathrm{t}}+\epsilon_{\mathrm{rec}}\dot{M}_{\mathrm{out}},$ (13)

where $F_{\mathrm{t}}$ is the total magnetohydrodynamic energy flux,

$F_{\mathrm{t}}=\int\limits_{v_{r}>0}r^{2}\left[(\rho
e+P_{\mathrm{t}}+\rho\Phi)v_{r}+\frac{\mathbf{B}(\mathbf{v}\cdot\mathbf{B})}{4\pi}\right]\,\mathrm{d}\Omega.$
(14)

In Figure 6, we compare
$F_{\mathrm{t}}+\epsilon_{\mathrm{rec}}\dot{M}_{\mathrm{out}}$ to
$\dot{E}_{\mathrm{expl}}$ using a value of
$\epsilon_{\mathrm{rec}}=6\,\mathrm{MeV}/m_{\mathrm{b}}$ (where
$m_{\mathrm{b}}$ is the atomic mass unit), which accounts for incomplete
recombination as well as some turbulent mixing. The growth of
${E}_{\mathrm{expl}}$ is described reasonably well by the assumption of
delivery by nuclear recombination and an additional, subdominant contribution
by $P\,\mathrm{d}V$-work and magnetic stresses, although
$\epsilon_{\mathrm{rec}}\dot{M}_{\mathrm{out}}$ alone somewhat underestimates
the actual growth rate.

To quantify the direct role of magnetic fields in powering the explosion, we
also compare the total magnetohydrodynamic energy flux $F_{\mathrm{t}}$ to the
purely hydrodynamic total enthalpy flux $F_{\mathrm{h}}$,

$F_{\mathrm{h}}=\int\limits_{v_{r}>0}r^{2}(\rho
e+P+\rho\Phi)v_{r}\,\mathrm{d}\Omega,$ (15)

where the magnetic pressure and magnetic stresses are excluded and only the
gas pressure $P$ enters aside from purely advective terms. At the base of the
ejecta region at black $100\,\mathrm{km}$ (where the angle-averaged total
energy in the outflows becomes positive), the black hydrodynamic enthalpy flux
$F_{\mathrm{h}}$ accounts for most of the total energy flux $F_{\mathrm{t}}$
into the ejecta. The rest of the energy flux into the ejecta is due to work by
magnetic pressure and magnetic stresses. Their contribution is non-negligible,
but smaller than the hydrodynamics flux $F_{\mathrm{h}}$ at the base of the
ejecta region.

The rate of growth of the explosion energy is therefore well explained by
nuclear recombination and some additional magnetohydrodynamic energy flux into
the ejecta region. The unusually rapid growth of the explosion energy is
merely the result of a rather high mass outflow rate, and this in turn is
simply due to the early onset of the explosion at a time when the mass
accretion rate is still high (close to $1M_{\odot}\,\mathrm{s}^{-1}$) and the
mass of the gain region is large. At first glance, the fact that
$\dot{E}_{\mathrm{expl}}$ even exceeds the volume-integrated neutrino heating
rate $\dot{Q}_{\nu}$ in the gain region (Figure 6) might suggest that
additional powering other than neutrino heating is needed to explain the high
mass outflow rate. For a steady-state outflow one would expect a ratio
$\dot{E}_{\mathrm{expl}}/\dot{Q}_{\nu}\approx\epsilon_{\mathrm{nuc}}/|e_{\mathrm{bind}}|$,
where $|e_{\mathrm{bind}}|$ is the typical binding energy of neutrino-heated
matter before it turns away from the PNS to be ejected, and if
$|e_{\mathrm{bind}}|$ is taken as the net binding energy at the gain radius,
one usually obtains $\dot{E}_{\mathrm{expl}}/\dot{Q}_{\nu}<1$. The rather
small magnetohydrodynamic energy flux of only $\mathord{\sim}3\times
10^{51}\,\mathrm{erg}\,\mathrm{s}^{-1}$ deeper in the gain region at radii of
less than $100\,\mathrm{km}$ in Figure 7 argues against an additional power
source for lifting the ejecta to roughly zero net energy before recombination.
The large ratio $\dot{E}_{\mathrm{expl}}/\dot{Q}_{\nu}$ is instead explained
by other effects. First, Müller et al. (2017) pointed out that the effective
turnaround radius of the ejecta can become quite large in a developing
explosion, so that material need not be lifted out of the potential well all
the way from the gain radius. Second, the average binding energy in the gain
black region decreases quite quickly in the strong-field model as shown in
Figure 10. Third, for a non-stationary outflow, $\dot{E}_{\mathrm{expl}}$ and
$\dot{Q}_{\nu}$ cannot be compared at the same time because there is a delay
between peak neutrino heating of an ejected mass element in the gain region
and the time when it acquires contributes to the explosion energy after
undergoing recombination.

The dynamics of the outflows in the strong-field model thus clearly point
towards an _indirect_ effect of the initial field on the explosion energy. The
key role of the magnetic fields is to _trigger_ shock revival early, which
provides the conditions for a rapid growth of the explosion energy. We
therefore need to elucidate how the strong initial fields precipitate shock
revival.

###  4.2 Dynamics of the Gain Region and Evolution Towards Shock Revival

Contrary to the explosion energetics, the strong initial magnetic fields play
a clearly recognisable role in the dynamics of the gain region prior to and
around shock revival. The earlier onset of explosion in the strong-field model
appears to be partly due to similar effects as in Müller & Varma (2020), where
the magnetic fields in the gain region were generated by a turbulent dynamo,
but our analysis also reveals a new phenomenon that is peculiar to stellar
cores with strong pre-collapse fieldsblack .

To analyse the dynamics of the gain region, we first consider the critical
ratio of the advection and heating time scales
$\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}$, (Figure 8), which quantifies the
conditions for runaway shock expansion by neutrinos (Buras et al., 2006). The
advection and heating timescale $\tau_{\mathrm{adv}}$ and
$\tau_{\mathrm{heat}}$ are defined as

$\tau_{\mathrm{adv}}=\frac{M_{\mathrm{g}}}{\dot{M}},$ (16)
$\tau_{\mathrm{heat}}=\frac{|E_{\mathrm{g}}|}{\dot{Q}_{\nu}},$ (17)

where $E_{\mathrm{g}}$ and $M_{\mathrm{g}}$ are the binding energy and mass of
the gain region, $\dot{M}$ is the mass accretion rate, and $\dot{Q_{\nu}}$ is
the volume-integrated neutrino heating rate.

Figure 8: Ratio of the advection timescale $\tau_{\mathrm{adv}}$ to the
heating timescale $\tau_{\mathrm{heat}}$ in the gain region of the strong-
field model (black) and the weak-field model (red). A horizontal dotted line
indicates the threshold criterion $\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}=1$
for runaway shock expansion in the neutrino-driven mechanism.

It is noteworthy that by the time the shock starts to expand in the strong-
field model and deviates from the shock trajectory in the weak-field model at
$0.1\,\mathrm{s}$ after bounce, $\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}$
still has not reached the critical threshold
$\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}\approx 1$ and is actually slightly
_lower_ than in the weak-field model. The lower ratio
$\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}$ can be explained by the slightly
lower neutrino heating rates in the strong-field model (Figure 6; cp. the
discussion of the neutrino luminosities in Section 4.1). While
$\tau_{\mathrm{adv}}/\tau_{\mathrm{heat}}\approx 1$ is reached in the strong-
field model about $0.02\,\mathrm{s}$ later, this provides strong evidence that
dynamical effects of magnetic fields initially aid shock expansion (even
though their importance as a subsidiary driver of the explosion diminishes
once the runaway criterion is met).

We evidence for important dynamical effects of the strong seed fields both
inside the gain region (as in Müller & Varma 2020) and outside the shock.
Their impact inside the gain region can be gauged by the convective efficiency
parameter for the conversion of neutrino heating to turbulent kinetic and
magnetic energy in the gain region, (Müller & Janka, 2015b; Müller & Varma,
2020),

$\eta_{\mathrm{conv}}=\frac{(E_{\mathrm{kin}}+E_{\mathrm{\mathbf{B}}})/M_{\mathrm{gain}}}{[(r_{\mathrm{sh}}-r_{\mathrm{gain}})\dot{Q}_{\nu}/M_{\mathrm{gain}}]^{2/3}},.$
(18)

As shown by Figure 9, $\eta_{\mathrm{conv}}$ is $\mathord{\lesssim}50\%$
higher than in the weak-field model until about $0.1255\,\mathrm{s}$ post-
bounce. The convective efficiency computed from the turbulent kinetic energy
alone roughly equals the total convective efficiency in the weak-field model,
in which the turbulent magnetic energy is negligible during this phase. The
situation is thus analogous to that in Müller & Varma (2020); the turbulent
magnetic energy is filled as an extra energy reservoir on top of the kinetic
turbulent energy, which is fed by the neutrino heating. Different from Müller
& Varma (2020), the increase in convective efficiency is observed as soon as
convection develops in the strong-field model, however.

Figure 9: Evolution of the convective efficiency $\eta_{\mathrm{conv}}$ in the
gain region as defined by Equation (18). The solid lines show the convective
efficiency based on the total kinetic and magnetic turbulent energy
$E_{\mathrm{kin}}+E_{\mathrm{mag}}$, while only the turbulent kinetic energy
$E_{\mathrm{kin}}$ is included for dash-dotted lines. Black and red curves are
used for the strong- and weak-field model, respectively. Only the solid line
is visible for the weak-field model because the magnetic energy in the gain
region is negligible.

The early boost to the convective efficiency is due to the strong seed fields,
which allow the magnetic fields in the gain region to reach levels close to
equipartition with little amplification, compared to the situation with weak
seed fields, where this takes an appreciable amount of time despite
exponential amplification by a turbulent dynamo (Müller & Varma, 2020). This
is illustrated by Figure 11, which shows the turbulent kinetic energy
$E_{\mathrm{kin}}$ and the magnetic field energy $E_{B}$ in the gain region.
$E_{\mathrm{kin}}$ and $E_{\mathrm{kin}}$ are computed as

$\displaystyle E_{\mathrm{kin}}$
$\displaystyle=\int\frac{1}{2}\rho|\mathbf{v}^{\prime}|^{2}\,\mathrm{d}V,$
(19) $\displaystyle E_{B}$
$\displaystyle=\int\frac{|\mathbf{B}|^{2}}{8\pi}\,\mathrm{d}V,$ (20)

where $\mathbf{v^{\prime}}$ denotes the fluctuations of the velocity around
its spherical Favre average black (i.e.
$\mathbf{v^{\prime}}=\mathbf{v}-\langle\mathbf{v}\rangle$). We also show the
ratio $E_{B}/E_{\mathrm{kin}}$ of the magnetic and kinetic turbulent energy,
in Figure 12. The strong initial magnetic field strength model undergoes a
very brief exponential magnetic field amplification process, placing the ratio
of $E_{B}/E_{\mathrm{kin}}\approx 40\%$ when runaway expansion occurs, similar
to the saturation value found in Müller & Varma (2020). By contrast, although
the weak magnetic field model undergoes a rapid magnetic field amplification
process in the gain region after $\mathord{\approx}0.20\,\mathrm{s}$, the
neutrino-driven mechanism launches the shock before the magnetic energies are
strong enough to affect the dynamics of the gain region. black We note that
there is rough energy equipartition between the radial magnetic fields
(generated predominantly by radial stretching of field lines) and the non-
radial components (generated predominantly by non-radial shear flows). The
energies $E_{B,r}$ and $E_{B,\theta\varphi}$ in the radial and non-radial
components,

$\displaystyle E_{B,r}$
$\displaystyle=\int\frac{B_{r}^{2}}{8\pi}\,\mathrm{d}V,$ (21) $\displaystyle
E_{B,\theta\varphi}$
$\displaystyle=\int\frac{B_{\theta}^{2}+B_{\varphi}^{2}}{8\pi}\,\mathrm{d}V,$
(22)

are shown in Figure 12, as a ratio with the kinetic energy in the gain region.
For approximate equipartition, the dynamical effect of magnetic fields on
shock propagation through tension and pressure _forces_ is arguably quite
neutral, but storing extra _energy_ in magnetic fields will still make it
easier to unbind material and release additional energy back into the
neutrino-heated material as it expands. The detailed dynamical impact of the
magnetic fields on the quasi-hydrostatic structure of the gain region still
remains to be investigated more rigorously; we recall that the chain of
effects through which turbulence facilitates shock expansion is highly non-
trivial even in the hydrodynamics case (e.g., Mabanta & Murphy, 2018).

In summary, a similar mechanism operates in the strong-field model as in the
MHD-aided explosion of Müller & Varma (2020); magnetic fields in the gain
region become sufficiently strong to boost the overall turbulent energy in the
gain region and thereby aid shock revival. The key difference is that this
process occurs significantly earlier due to the strong seed fields. As in
Müller & Varma (2020), it is not possible to disentangle precisely which
turbulent terms (e.g., the magnetic stresses in the momentum equation, the
flux terms in the energy equation, etc.) in the MHD equations are responsible
for precipitating shock expansion. Based on the ratio
$E_{B}/E_{\mathrm{kin}}$, the magnetic stress terms are clearly non-
negligible. There is somewhat better evidence than in Müller & Varma (2020)
that magnetic fields also palpably reduce the binding energy of the gain
region (Figure 10). $E_{\mathrm{bind}}$ in the strong-field models starts to
visibly deviate from the weak-field model from $\approx 0.09\,\mathrm{s}$
after bounce, i.e. already slightly before the shock trajectories diverge.

Figure 10: Average specific binding energy $e_{\mathrm{bind}}$ in the entire
gain region of the strong-field model (black) and the weak-field (model)
black.  Figure 11: Evolution of the total magnetic (solid) and turbulent
kinetic (dotted) energy within the gain region (Equations 19 and 20) for the
strong-field model (black) and the weak-field model (red).  Figure 12:
Evolution of the ratio between the total magnetic and turbulent kinetic energy
within the gain region for strong-field model (black) and the weak-field model
(red). black For both models, the ratio is plotted for the total magnetic
energy (solid), the radial only component (dashed), and the angular components
(dotted) of the magnetic energy.

(a) Strong:0.01s

(b) Weak:0.04s

(c) Strong:0.05s

(d) Weak:0.10s

(e) Strong:0.09s

(f) Weak:0.13s

Figure 13: Relative ram pressure perturbations $P^{\prime}$ ahead of the shock
for the strong-field model (left column) and the weak-field model (right
column) at a radius of $200\,\mathrm{km}$ at different times (denoted below
each panel). The final row shows the ram pressure perturbations amplitudes
immediately before the maximum shock radius reaches $200\,\mathrm{km}$.

(a) (a) (b) (b)

Figure 14: Spectrum $\hat{P}^{\prime}_{\mathrm{ram},\ell}$ of the ram pressure
perturbations (Equation 24) at $200\,\mathrm{km}$, ahead of the shock, (a) for
the strong-field model and (b) for the weak-field model.

(a) Strong

(b) Weak

Figure 15: black Entropy on meridional slices through the strong (left) and
weak-field (right) model at $0.30\,\mathrm{s}$ after bounce. The strong-field
model has a weakly bipolar explosion geometry with a prominent plume along the
South polar axis and a weaker outflow along the North polar axis. The biggest
neutrino-heated plume in the weak-field model is also aligned with the axis,
but there are smaller plumes in other directions, most notably near the
equatorial plane in the 4 o’clock direction.

In addition to increasing the turbulent in the gain region by jump-starting
turbulent field amplification, the strong initial fields affect the dynamics
by modifying the pre-shock conditions, similar to what occurs for
“perturbation-aided explosions” in progenitors with vigorous convection in the
oxygen or silicon shell (Couch et al., 2015; Müller & Janka, 2015b; Müller et
al., 2017). In the perturbation-aided mechanism, vortical convective motions
are translated to sizable perturbations of the pre-shock density and ram
pressure by advective-acoustic coupling (Müller & Janka, 2015b; Abdikamalov &
Foglizzo, 2020; Abdikamalov et al., 2021), which then facilitate asymmetric
shock expansion. We demonstrate that the initial fields assumed in our models
also lead to differential infall and large-scale asymmetries in the pre-shock
ram pressure. To this end, we consider the relative ram pressure perturbations
$P^{\prime}_{\mathrm{ram}}$,

$P^{\prime}_{\mathrm{ram}}=\frac{P_{\mathrm{ram}}-\langle
P_{\mathrm{ram}}\rangle}{\langle P_{\mathrm{ram}}\rangle},$ (23)

where $P_{\mathrm{ram}}=\rho v_{r}^{2}$ is the angle-dependent ram pressure,
and $\langle P_{\mathrm{ram}}\rangle$ is its average over solid angle. We
evaluate $P^{\prime}_{\mathrm{ram}}$ at a radius of $200\,\mathrm{km}$, which
is small enough to capture the pre-shock conditions sufficiently well but
larger than the maximum shock radius before shock revival. 2D Hammer
projection of the ram pressure perturbations at this radius are shown for the
strong-and weak-field models at different times in Figure 13. In the strong-
field model, large-scale asymmetries in the pre-shock ram pressure clearly
emerge already at $0.05\,\mathrm{s}$ and grow considerably stronger shortly
before shock revival at $0.09\,\mathrm{s}$. The contrast between the
directions of high and low ram pressure reaches about $20\%$, and there is a
clear black large-scale asymmetry with a ram pressure deficit in the Southern
hemisphere, an excess at mid-latitudes in the Northern hemisphere, and a
deficit in a smaller region above latitudes of $\gtrsim 60^{\circ}$ in the
Northern hemisphere. Such large-scale perturbations are most conducive to
perturbation-aided shock revival (Müller & Janka, 2015b). The large-scale
component of the perturbations is in fact more clearly visible than in
simulations based on 3D progenitor models with oxygen shell convection (cp.
Figure 8 in Müller et al. 2017), while the magnitude of relative pressure
perturbation are of a similar magnitude. In the weak-field model, the magnetic
fields are not strong enough to dynamically affect the infall. Weaker pre-
shock ram pressure perturbations emerge from numerical noise in the infall
region in this case. These become weaker in time, and are dominated by small-
scale structures. Figure 14 shows a spherical harmonics decomposition of
$P^{\prime}_{\mathrm{ram}}$ to illustrate the different magnitude and scale of
the ram pressure perturbations in the two models more quantitatively. We
compute the spectrum $\hat{P}_{\mathrm{ram},\ell}$ as

$\hat{P}_{\mathrm{ram},\ell}=\sqrt{\sum_{m=-\ell}^{\ell}\left|\int Y_{\ell
m}^{\star}P^{\prime}_{\mathrm{ram}}\,\mathrm{d}\Omega\right|^{2}}.$ (24)

The spherical harmonic functions are orthonormalised and calculated using the
SHTOOLS package (Wieczorek et al., 2016). The presence of a dipole black and
quadrupole of growing amplitude in the strong-field model is clearly visible
in Figure 14. black The precise relation of the pre-shock ram pressure
perturbations to the initial magnetic field structure is not yet fully
transparent, and needs to be investigated further using linear perturbation
theory in future (see discussion below). In particular, the emergence of an
appreciable North-South asymmetry, which could be related to an instability of
the initial twisted-torus field configuration, warrants further attention.

black As shown by Figure 15, the forcing by asymmetries in ram pressure is
reflected in the explosion geometry of the strong-field model. The most
prominent high-entropy plume aligns with the region of lower ram pressure in
the Southern hemisphere. There is a somewhat weaker neutrino-heated plume
around the North polar axis as well, corresponding to the ram pressure deficit
above latitudes of $60^{\circ}$ in the bottom left panel of Figure 13. The
prominent downflows at moderate latitudes in the Northern hemisphere are
aligned with the regions of excess ram pressure. In the weakly magnetised
model, the explosion geometry is characterised by a big plume along the North
polar axis and several less prominent plumes. In this case, the emergence of a
prominent polar plume is likely facilitated by the grid geometry; a tendency
towards grid-aligned explosion in weakly perturbed models has been observed
before with the CoCoNuT code, though it is not universal (see, e.g., model
he3.5 in Müller et al. 2019b). There is no clear relation to the pre-shock ram
pressure perturbations.

It should be noted, however, that the mechanism whereby magnetic fields are
converted to acoustic perturbations during collapse still needs to be
investigated in future work along the lines of existing linear theory for the
amplification of hydrodynamic perturbations (Lai & Goldreich, 2000; Takahashi
& Yamada, 2014; Abdikamalov & Foglizzo, 2020; Abdikamalov et al., 2021). black
The nature of the magnetohydrodynamic perturbations in the progenitor warrants
further scrutiny. black In the magnetohydrodynamic case, the linear evolution
during collapse can be expected to be richer with interactions between between
different waves families beyond the mechanism of vortical-acoustic and
entropic-acoustic coupling in the hydrodynamic case.

black There is also the question whether the twisted-torus initial conditions
adequately approximate the magnetic field structure in the progenitor. For the
hydrodynamic case, it is important to take into account that the initial
perturbation spectrum from weakly compressible convective flow is dominated by
vorticity and entropy perturbations and not by acoustic perturbations (Müller
& Janka, 2015b; Abdikamalov & Foglizzo, 2020). In the absence of consistent 3D
progenitor models, this can adequately be captured by imposing a solenoidal
flow geometry in convective regions in the progenitor. In the fossil-field
scenario, it is less clear whether imposing a reasonable twisted-torus field
geometry already captures the relevant physical constraints on the wave
structure of the pre-collapse perturbations. black This is particularly
important for the strong-field case, where the magnetic fields are strong
enough to noticeably perturb hydrostatic equilibrium with initial values of a
few $10^{-2}$ for the ratio of magnetic to thermal pressure in the oxygen
shell. For this reason, MHD models of quasi-equilibria with strong fossil
fields in supernova progenitors are an important desideratum for the future.

###  4.3 Neutron Star Properties

We next consider the birth properties of the neutron star produced in the two
simulation. These are particularly relevant for evaluating the prospects for
the fossil-field scenario for magnetar formation, although it must be borne in
mind that the magnetic field and spin rate of the neutron star can still be
affected by processes that occur on longer time scales than we can simulate
here, such as spin-up and field burial due to fallback.

We show the bulk properties of the PNS in Figure 16, namely is radius
${R_{\mathrm{PNS}}}$, mass ${M_{\mathrm{PNS}}}$, and magnetic energy
${E_{\mathrm{PNS}}}$, which we compute as volume integrals over the region
where $\rho>10^{11}\mathrm{g\,cm^{-3}}$. Initially, the PNS mass and radius
for both models are largely identical until
$\mathord{\approx}0.15\,\mathrm{s}$ aside from minimal differences resulting
from small differences during the collapse phase. After this point,
$M_{\mathrm{PNS}}$ continues to grow in the weak-field model and the model
maintains a larger $R_{\mathrm{PNS}}$ since accretion continues longer than in
the strong-field model. The final baryonic masses are $1.63\,M_{\odot}$ and
$1.69\,M_{\odot}$ for the strong- and weak-field model, respectively. Using
the approximation of Lattimer & Yahil (1989); Lattimer & Prakash (2001) for
the neutron star binding energy, the gravitational mass $M_{\mathrm{grav}}$
can be obtained from

$M_{\mathrm{by}}=M_{\mathrm{grav}}+0.084M_{\odot}(M_{\mathrm{grav}}/M_{\odot})^{2},$
(25)

which results in $M_{\mathrm{grav}}=\,1.46M_{\odot}$ for the strong-field
model and $M_{\mathrm{grav}}=\,1.50M_{\odot}$ for the weak-field model.

The evolution of the interior and surface field of the PNS is of particular
interest in the context of magnetar and pulsar birth scenarios. Due to the
vast differences in initial magnetic field strength, the initial
$E_{\mathrm{PNS}}$ (bottom panel of Figure 16) differs by orders of magnitude
between the strong-and weak-field models. For the strong-field model
$E_{\mathrm{PNS}}$ remains rather steady at
$10^{47\texttt{-}48}\,\mathrm{erg}$ until about $0.25\,\mathrm{s}$ after
bounce. In the weak-field model, $E_{\mathrm{PNS}}$ starts at much smaller
values, but gradually grows due to compression of the field lines as the
$\mathrm{R_{PNS}}$ contracts, shear flows in the PNS surface region, and
amplification by PNS convection. Eventually $E_{\mathrm{PNS}}$ becomes
comparable to the strong-field model about $0.25\,\mathrm{s}$ after bounce. It
is noteworthy that by the end of the simulation, the magnetic energy in the
PNS is dropping in the strong-field model and ends up slightly lower than in
the weak-field model. The drop in magnetic energy for the strong-field model
is likely due to numerical reconnection and variations in field amplification
by PNS convection. While turbulent reconnection is bound to occur in the PNS
convection zone, the destruction of magnetic fields may be overestimated
numerically, especially at the modest resolution affordable inside the PNS.
Nonetheless, the simulations at least strongly indicate that the memory of the
initial fields is largely lost inside the PNS.111PNS fields originating from
field with significant helicity might still be protected better from field
decay in the long run, but once PNS convection has died down and differential
rotation has been eliminated, field decay will be extremely slow because of
the high conductivity of neutron star matter, regardless of the helicity of
the field configuration.

To elucidate the spatial structure of the PNS magnetic field, we plot profiles
of the dipole strength of the radial magnetic field and the root-mean-squared
average of the total field strength against density for the strong- and weak-
field models in Figure 17. The tentative “surface” of the PNS where
$\rho=10^{11}\,\mathrm{g}\,\mathrm{cm}^{-3}$ is indicated by dotted vertical
line, but it must be borne in mind that further accretion, ablation of surface
material by the neutrino-driven wind, and possibly field breakout from the PNS
may yet change the PNS surface structure and surface magnetic fields in
particular. Just as for the bulk magnetic energy of the PNS, the PNS surface
magnetic field has effectively lost memory of the initial condition after less
than $0.3\,\mathrm{s}$ after bounce, even though the initial magnetic field
strengths of both these models were six orders of magnitude apart. At the
final simulation time, the surface magnetic field strengths in both models are
very similar, with dipole fields of $5\times 10^{14}\,\mathrm{G}$ for the
weak-field case and $2\times 10^{14}\,\mathrm{G}$ for the strong field case.
Thus, both models reach magnetar-strength fields at the surface as in previous
3D MHD supernova simulations of non-rotating and slowly-rotating progenitors
(Müller & Varma, 2020; Matsumoto et al., 2022), at least transiently. It is
interesting, however, that the surface field in the strong-field model is
decreasing at the end of the simulation, which may reduce possible tensions
with the observed field strength distribution of young neutron stars. The
declining field strength is consistent with the development of net outflow
from the PNS surface, i.e., an incipient neutrino-driven wind, which implies a
change from flux compression to flux expansion. The final surface field
strength is obviously still beyond the scope of the current simulations.

Figure 16: Evolution of the radius $R_{\mathrm{PNS}}$, baryonic mass
$M_{\mathrm{PNS}}$, and magnetic energy $E_{\mathrm{PNS}}$ of the PNS for the
strong-field model (black) and the weak-field model (red). The PNS surface is
defined as the radius where the angle-averaged density reaches
$10^{11}\,\mathrm{g}\,\mathrm{cm}^{-3}$.

(a) (a) (b) (b)

Figure 17: Angle-averaged total magnetic field strength (solid lines, root-
mean-square averages), and the dipole component of the radial magnetic field
$B_{r}$ (dotted lines) as a function of density (a) for the strong-field model
and (b) for the weak-field model. The vertical dotted black lines indicate the
fiducial PNS surface density $\rho=10^{11}\mathrm{g\,cm^{-3}}$.  Figure 18:
Evolution of the kick and angular momentum of PNS for the strong-field model
(black/blue) and the weak-field model (red). For the strong-field model, we
show both the time integral of the entire angular momentum flux onto the PNS
(black), and the time integral of magnetic torques alone (blue). The advective
angular momentum flux and the magnetic torques are evaluated on a sphere of
radius $30\,\mathrm{km}$ around the origin. Different line styles are used for
the three components of the angular momentum vector (solid: $J_{x}$, dash-
dotted: $J_{y}$, dotted: $J_{z}$).

We next consider the kick velocity and angular momentum of the newly born
neutron star, which we plot in Figure 18. Following Scheck et al. (2008), we
evaluate the kick indirectly from the momentum of the ejecta assuming angular
momentum conservation since the innermost part of our grid is treated in
spherical symmetry and does not allow the PNS to move freely. Specifically, we
compute

$\mathbf{v}_{\mathrm{kick}}=-\frac{1}{M}\int_{r>R_{\mathrm{gain}}}\rho\mathbf{v}\,\mathrm{d}V,$
(26)

where $M$ is the total mass222Although the choice between the baryonic and
gravitational mass would only make a minor difference here, we note that it is
appropriate to use the baryonic PNS mass for evaluating the kick in a pseudo-
Newtonian simulation. Although the inertial mass of the neutron star is
eventually reduced by neutrino emission as it cools, near-isotropic emission
of neutrinos in the PNS rest frame will not change the PNS velocity, and
decrease the PNS momentum proportionately to the decrease in inertial mass.
inside the gain radius, and the integral covers the entire volume outside the
gain radius. black In both models, the kick is almost aligned with the
$z$-axis of the spherical polar grid, reflecting the explosion geometry shown
in Figure 15, and thus the forcing by initial perturbations and/or the
spherical polar grid geometry (cp. discussion in Section 4.2). As expected,
the kick points in the direction opposite to the most prominent plume, i.e.,
towards the North pole in the strong-field model and towards the South pole in
the weak-field model, at least by the end of the simulation. Since the
explosion in the weak-field model has not yet sufficiently advanced to
extrapolate the final value of the kick, which has only reached
$84\,\mathrm{km}\,\mathrm{s}^{-1}$ by the end of the simulation. black The
kick direction may also still be subject to change, since other non-grid
aligned plumes and downflows are still present in the weak-field model and may
impart momentum onto the PNS in other directions. In the strong-field model,
the kick reaches $326\,\mathrm{km}\,\mathrm{s}^{-1}$, and although it is still
growing, its rate of change is decreasing. Since accretion downflows have
already been quenched and the ejecta are expanding steadily, the final kick
will only grow further by a modest amount. We can estimate the final kick
velocity of these models using the extrapolation method demonstrated in Müller
et al. (2019b) to find final PNS kick velocities of
$345\,\mathrm{km}\,\mathrm{s}^{-1}$ and $110\,\mathrm{km}\,\mathrm{s}^{-1}$
for the strong and weak-field simulations, respectively. We thus expect the
final kick to lie within the typical range observed for young neutron stars
(Hobbs et al., 2005).

The second panel in Figure 18 shows the angular momentum $\mathbf{J}$ of the
PNS, which is obtained by integrating the magnetic torques and hydrodynamic
angular momentum fluxes onto the PNS in time and over the surface of a sphere
of radius $30\,\mathrm{km}$,

$\frac{\mathrm{d}\mathbf{J}}{\mathrm{d}t}=-\oint(\mathbf{r}\times\mathbf{v})\rho
v_{r}\,\mathrm{d}A+\frac{1}{4\pi}\oint(\mathbf{r}\times\mathbf{B})B_{r}\,\mathrm{d}A.$
(27)

For the strong-field model, the contribution from magnetic torques is also
shown separately. We note, however, that our choice of the PNS surface being
at $30\,\mathrm{km}$ is somewhat arbitrary; other works such as Wongwathanarat
et al. (2013) have used a considerably larger analysis region up to radii of
several hundred kilometres. The choice of a small radius is designed to better
capture the current PNS spin, but the final spin of the PNS depends rather
sensitively on how much of the material in its vicinity is eventually accreted
or blown out333The sensitivity of the PNS spin to rather small amounts of
accretion has been demonstrated most perspicuously in the extreme case of
fallback accretion at later times, which can easily spin up the PNS to
millisecond periods (Chan et al., 2018; Chan et al., 2020; Stockinger et al.,
2020).. The eventual multi-dimensional “mass cut” cannot be predicted at this
rather early stage of the explosion yet. This introduces a significant
uncertainty on the final spin, since the shells outside $30\,\mathrm{km}$
still carry considerable angular momentum and have not yet been brought into
rigid rotation with the PNS by the end of the simulation. Experiments with
different analysis radii (corresponding to different hypothetical mass cuts)
indicate that the PNS angular momentum could still change by a small factor,
with the direction of the spin appearing more uncertain than the spin
magnitude. By contrast, $\mathbf{v}_{\mathrm{kick}}$ is quite insensitive to
the precise choice of the analysis region. The results for the weak-field
model are only shown for the sake of completeness; as for the kick, the
angular momentum shows no sign of saturation in this case yet. Its absolute
value at the end of the simulation is $5.11\times
10^{45}\,\mathrm{g}\,\mathrm{cm}^{2}\,\mathrm{s}^{-1}$ and is still growing.

In the strong-field model, the angular momentum of the PNS has largely frozen
out, however. By the end of the simulation, it has reached about $1.36\times
10^{47}\,\mathrm{g}\,\mathrm{cm}^{2}\,\mathrm{s}^{-1}$. It is interesting to
note that the spin-up of the PNS in the strong-field model is driven mostly by
magnetic torques.

Using the analytic approximation for the neutron star moment of inertia from
Lattimer & Schutz (2005),

$I\approx
0.237M_{\mathrm{grav}}R^{2}\left[1+4.2\left(\frac{M_{\mathrm{grav}}}{M_{\odot}}\frac{\mathrm{km}}{R}\right)+90\left(\frac{M_{\mathrm{grav}}}{M_{\odot}}\frac{\mathrm{km}}{R}\right)\right],$
(28)

the angular momentum can be translated into a neutron star spin period
$P=J/I$. Assuming a final neutron star radius $R$ of $12\,\mathrm{km}$, we
estimate the moment of inertia as $\mathrm{1.52\times 10^{45}g\,cm^{3}}$ for
the strong-field model and $\mathrm{1.58\times 10^{45}g\,cm^{3}}$ for the
weak-field model, which corresponds to spin periods of about $73\,\mathrm{ms}$
and $1.86\,\mathrm{s}$, respectively. Prima facie, this puts the strong-field
model in the range of typical spin-period of young _pulsar_ (Faucher-Giguère &
Kaspi, 2006; Perna et al., 2008; Popov & Turolla, 2012; Igoshev & Popov,
2013), and not of magnetars, which are observed to rotate rather slowly
(Olausen & Kaspi, 2014; Kaspi & Beloborodov, 2017). However, highly magnetised
neutron stars will also spin-down rapidly; even the youngest observed
magnetars must already have been spun down considerably after the supernovae.
We therefore extrapolate the spin period using the neutron star spin-down
formula for the period derivative $\dot{P}$ in terms of $P$ and the surface
dipole field $B_{\mathrm{surf}}$ (cp. Lorimer & Kramer, 2004),

black

$\dot{P}=\frac{1}{P}\left(\frac{B_{\mathrm{surf}}}{3\times
10^{19}}\right)^{2}\,\mathrm{s^{-1}}.$ (29)

Since the final surface dipole field strength is still subject to
uncertainties, we calculate the spin-down for several choices of
$B_{\mathrm{surf}}$ in the range around the tentative dipole field strength of
$2\times 10^{14}\,\mathrm{G}$ at the end of our simulation. The extrapolated
spin period after several hundred years for the strong-field model (Figure 19)
reach several seconds, which is compatible with observations (Olausen & Kaspi,
2014), which show magnetar spin periods of $1\,\mathrm{s}<P<12\,\mathrm{s}$.
The results can, e.g., be compared to Swift J1818.0–1607 (Blumer & Safi-Harb,
2020), one of the youngest known magnetars, which is approximately 500 years
old and has a spin period of $1.36\,\mathrm{s}$. For a surface dipole field
strength of $10^{14}\,\mathrm{G}$, we obtain a spin period of
$\mathord{\sim}0.6\,\mathrm{s}$ after 500 years. While this is still faster
than Swift J1818.0–1607, it is of the same magnitude, and thus appears roughly
compatible with observed constraints on magnetar spins, especially considering
that stochastic variations among strongly magnetised neutron stars are
expected.

Finally, we consider the orientation of the kick and spin direction.
Observational evidence has suggested generic spin-kick alignment in neutron
stars (Johnston et al., 2005; Noutsos et al., 2012, 2013), but simulations
have yet to reveal a robust mechanism for spin-kick alignment; recent attempts
to connect spin-kick alignment to rotation (Powell & Müller, 2020) or fallback
(Janka et al., 2022) have not yielded wholly convincing results. At least
during the early explosion phase, magnetic fields in non-rotating supernova
models do not appear to provide a mechanism for spin-kick alignment either, as
can be seen from Figure 20, which shows the angle $\alpha$ between the spin
and kick direction for our models. We find a misalignment of about
$50^{\circ}$ and $35^{\circ}$ at the end of the simulation for the strong-and
weak-field models. With supernova explosion models now including magnetic
fields, it becomes increasingly difficult to account for spin-kick alignment
during the early explosion phase. Spin-kick alignment may indeed be a long-
term effect and related, e.g., to later fallback (Janka et al., 2022)
(although the mechanism proposed by Janka et al. 2022 likely does not work as
intended because of a problem with angular momentum non-conservation).

Figure 19: Rotation period of the neutron star as a function of age using the
spin-down formula (29). Figure 20: The angle $\alpha$ between the spin and
kick directions of the PNS based on an evaluation of the PNS angular momentum
using the hydrodynamic fluxes and magnetic torques at a radius of
$30\,\mathrm{km}$.

##  5 Conclusions

We investigated the role of strong magnetic fields in non-rotating progenitors
on core-collapse supernova explosions, specifically its effect on shock
revival, explosion dynamics, and the properties of the compact remnant. To
this end, we performed 3D simulations of the collapse and explosion of a
$16.9M_{\odot}$ star using a twisted-torus initial field and two different
initial field strengths ($10^{12}\,\mathrm{G}$ and $10^{6}\,\mathrm{G}$ at the
centre). It has been speculated that strong pre-supernova fields of order
$10^{12}\,\mathrm{G}$ might arise in massive stars as fossil fields from a
stellar merger (Ferrario & Wickramasinghe, 2005, 2006; Schneider et al., 2019;
Schneider et al., 2020).

Our simulations show that strong magnetic fields have a noticeable impact on
the neutrino-driven explosion mechanism, even in slowly/non-rotating
progenitors. The strong-field model undergoes shock revival already early on
about $0.1\,\mathrm{s}$ after bounce and reaches an explosion energy
$9.3\times 10^{50}\,\mathrm{erg}$ by the end of our simulation at
$0.3\,\mathrm{s}$ after bounce. Thus, the model reaches a typical supernova
explosion energy faster than in non-magnetised 3D models of core-collapse
supernovae from slowly rotating progenitors. While the corresponding weak-
field model also explodes somewhat later, it shows a much slower rise of the
explosion energy to only $1.2\times 10^{50}\,\mathrm{erg}$ by the end of the
simulation.

The principal role of magnetic fields in the precipitous and powerful
explosion of the strong-field model is to provide an extra boost to neutrino
heating and hydrodynamic turbulence to trigger the explosion early, while the
explosion is still primarily driven by neutrinos. This is similar to the
results of Müller & Varma (2020), who already found that field amplification
by a turbulent dynamo can aid shock revival by the neutrino-driven mechanism;
the key difference is that the effect of magnetic fields comes into play early
on without the need for an extended growth phase. The saturation field
strength in the gain region, corresponding to about 40% of kinetic
equipartition is reached almost immediately once neutrino-driven convection
develops. In the weak-field model, field amplification by the turbulent dynamo
is too slow to have a substantial impact before a neutrino-driven explosion
develops

As in Müller & Varma (2020), magnetic fields aid shock revival in the strong-
field model by increasing the convective efficiency in the gain region and
reducing the binding energy of the gain region. In addition, the strong fields
in the progenitor also create sizable perturbations of $\pm 10\%$ in the pre-
shock ram pressure, reminiscent of the effect of infalling convective eddies
from shell burning in perturbation-aided explosions (Müller & Janka, 2015b;
Müller et al., 2017; Abdikamalov & Foglizzo, 2020). The perturbations are
dominated by a large-scale $\ell=1$ dipole mode. Large-scale ram-pressure
perturbations of such amplitudes tend to be strongly conducive to shock
revival (Müller & Janka, 2015b; Abdikamalov et al., 2016; Huete et al., 2018).
In future, the evolution of strong magnetic seed fields in the progenitor and
their impact on shock revival need to be investigated using a more consistent
approach by basing supernova simulations directly on models of
magnetoconvection in massive stars (Varma et al., 2021) and complementing
simulations with analytic theory for the amplification and dynamical role of
seed perturbations, similar to the hydrodynamic case (Takahashi & Yamada,
2014; Müller et al., 2016; Huete et al., 2018; Abdikamalov & Foglizzo, 2020;
Abdikamalov et al., 2021).

The high explosion energy of the strong-field model is primarily a consequence
of early shock revival. Its energetics is consistent with the physics known
from other models of neutrino-driven explosions (Marek et al., 2006; Müller,
2015; Bruenn et al., 2016; Müller et al., 2017). Nucleon recombination is
primarily, but not exclusively, responsible for the delivery of the explosion
energy. Due to early shock revival, the volume-integrated neutrino heating
rate is high and the specific binding energy of the gain region is modest.
These factors allow for a high mass outflow rate, a very efficient conversion
of neutrino heating into recombination energy in the ejecta region, and hence
a rapid growth of the explosion energy. As in early explosions in 2D
simulations (Bruenn et al., 2016), mechanical $P\,\mathrm{d}V$ work and
(specific to our model) magnetic stresses also contribute significantly to the
energy budget of the ejecta.

The strong-field model has evolved sufficiently far to tentatively extrapolate
some properties of the neutron star born in the supernova. Our simulation
predicts a gravitational neutron star mass of $1.50M_{\odot}$, a kick velocity
of about $345\,\mathrm{km}\,\mathrm{s}^{-1}$, and a birth spin period of about
$73\,\mathrm{ms}$. The dipole component of the magnetic field at the surface
of the proto-neutron star is about $2\times 10^{14}\,\mathrm{G}$ at the end of
the simulation, but is decreasing significantly at this point. Such a field
could spin down the neutron star to periods of second within a few hundred
years, which would be compatible with the spin period of young magnetars. We
do not find spin-kick alignment, and the mechanism that might bring about such
an alignment remains opaque.

However, while the raw neutron star parameters from the simulation appear to
be quite consistent with observational constraints on young neutron stars and
magnetars in particular, our simulations rather prompt more questions about
neutron star formation channels and birth properties. It remains unclear
whether the proto-neutron star in the strong-field would actually evolve into
a magnetar. The surface dipole field strength is decreasing by the end of the
simulation, and we cannot follow phenomena that may shape its final magnetic
fields on much longer time scales, such as the ablation of surface material by
the neutrino-driven wind, field burial (Romani, 1990; Torres-Forné et al.,
2016) or breakout of fields from the proto-neutron star convection zone.
Regardless of these uncertainties, it is remarkable, though, that the weak-
field model actually develops a strong surface dipole field by the end of the
simulation than the strong-field model. As far as the proto-neutron star
surface field is concerned, memory of the initial conditions appears to be
erased within a few hundred milliseconds after collapse. This casts some doubt
on the viability of the fossil-field scenario for magnetar formation. Whether
amplification and reconnection processes in the supernova core indeed make the
neutron star magnetic field independent of the initial field in the long term
requires further investigation, however, both because of potential numerical
uncertainties (resolution, treatment of the induction equation, etc.) and
because the initial field in our simulation was still prescribed by hand.
Aside from the question of neutron star surface fields, it is also difficult
to square the high explosion energy of the strong-field model with the
characteristics of supernova remnants around magnetars, which point to have
normal or modest explosion energies (Vink & Kuiper, 2006).

Nonetheless, our simulations underscore that it is important to consider the
role of magnetic fields in supernova explosions – beyond the more spectacular
scenario of magnetorotational supernovae – in order to more fully understand
the explosion dynamics and the progenitor-remnant connection. As with
hydrodynamic simulations, a more integrated view of magnetic field effects
from the progenitor stage to the remnant phase is required for this purpose.

## Acknowledgements

BM acknowledges support by ARC Future Fellowship FT160100035. This work is
based on simulations performed within 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. This work has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement No. 945806).

## Data Availability

The data underlying this article will be shared on reasonable request to the
authors, subject to considerations of intellectual property law.

## References

  * Abdikamalov & Foglizzo (2020) Abdikamalov E., Foglizzo T., 2020, MNRAS, 493, 3496
  * Abdikamalov et al. (2016) Abdikamalov E., Zhaksylykov A., Radice D., Berdibek S., 2016, MNRAS, 461, 3864
  * Abdikamalov et al. (2021) Abdikamalov E., Foglizzo T., Mukazhanov O., 2021, MNRAS, 503, 3617
  * Aguilera-Dena et al. (2020) Aguilera-Dena D. R., Langer N., Antoniadis J., Müller B., 2020, ApJ, 901, 114
  * Akiyama et al. (2003) Akiyama S., Wheeler J. C., Meier D. L., Lichtenstadt I., 2003, ApJ, 584, 954
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Bisnovatyi-Kogan et al. (1976) Bisnovatyi-Kogan G. S., Popov I. P., Samokhin A. A., 1976, Ap&SS, 41, 287
  * Blumer & Safi-Harb (2020) Blumer H., Safi-Harb S., 2020, ApJ, 904, L19
  * Bollig et al. (2021) Bollig R., Yadav N., Kresse D., Janka H.-T., Müller B., Heger A., 2021, ApJ, 915, 28
  * Braithwaite & Nordlund (2006) Braithwaite J., Nordlund Å., 2006, A&A, 450, 1077
  * Braithwaite & Spruit (2004) Braithwaite J., Spruit H. C., 2004, Nature, 431, 819
  * Bruenn et al. (2016) Bruenn S. W., et al., 2016, ApJ, 818, 123
  * Bugli et al. (2020) Bugli M., Guilet J., Obergaulinger M., Cerdá-Durán P., Aloy M. A., 2020, MNRAS, 492, 58
  * Buras et al. (2006) Buras R., Janka H. T., Rampp M., Kifonidis K., 2006, A&A, 457, 281
  * Burrows et al. (2007) Burrows A., Dessart L., Livne E., Ott C. D., Murphy J., 2007, ApJ, 664, 416
  * Cantiello et al. (2007) Cantiello M., Yoon S. C., Langer N., Livio M., 2007, A&A, 465, L29
  * Chan et al. (2018) Chan C., Müller B., Heger A., Pakmor R., Springel V., 2018, ApJ, 852, L19
  * Chan et al. (2020) Chan C., Müller B., Heger A., 2020, MNRAS, 495, 3751
  * Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
  * Couch et al. (2015) Couch S. M., Chatzopoulos E., Arnett W. D., Timmes F. X., 2015, ApJ, 808, L21
  * Dedner et al. (2002) Dedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T., Wesenberg M., 2002, Journal of Computational Physics, 175, 645
  * Donati & Landstreet (2009) Donati J.-F., Landstreet J. D., 2009, ARA&A, 47, 333
  * Duncan & Thompson (1992) Duncan R. C., Thompson C., 1992, ApJ, 392, L9
  * Endeve et al. (2010) Endeve E., Cardall C. Y., Budiardja Mezzacappa A., 2010, ApJ, 713, 1219
  * 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
  * Ferrario & Wickramasinghe (2005) Ferrario L., Wickramasinghe D. T., 2005, MNRAS, 356, 615
  * Ferrario & Wickramasinghe (2006) Ferrario L., Wickramasinghe D., 2006, MNRAS, 367, 1323
  * Ferrario et al. (2016) Ferrario L., Melatos A., Zrake J., 2016, Magnetic Field Generation in Stars. Springer New York, New York, NY, pp 81–113, doi:10.1007/978-1-4939-3550-5_4, https://doi.org/10.1007/978-1-4939-3550-5_4
  * Georgy et al. (2009) Georgy C., Meynet G., Walder R., Folini D., Maeder A., 2009, A&A, 502, 611
  * Georgy et al. (2012) Georgy C., Ekström S., Meynet G., Massey P., Levesque E. M., Hirschi R., Eggenberger P., Maeder A., 2012, A&A, 542, A29
  * Grunhut et al. (2017) Grunhut J. H., et al., 2017, MNRAS, 465, 2432
  * Gurski (2001) Gurski K., 2001, in APS Meeting Abstracts. p. D2.007
  * Heger et al. (2005) Heger A., Woosley S. E., Spruit H. C., 2005, ApJ, 626, 350
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Huete et al. (2018) Huete C., Abdikamalov E., Radice D., 2018, MNRAS, 475, 3305
  * Igoshev & Popov (2013) Igoshev A. P., Popov S. B., 2013, MNRAS, 432, 967
  * Janka et al. (2008) Janka H. T., Müller B., Kitaura F. S., Buras R., 2008, A&A, 485, 199
  * Janka et al. (2022) Janka H.-T., Wongwathanarat A., Kramer M., 2022, ApJ, 926, 9
  * Jardine et al. (2022) Jardine R., Powell J., Müller B., 2022, MNRAS, 510, 5535
  * Johnston et al. (2005) Johnston S., Hobbs G., Vigeland S., Kramer M., Weisberg J. M., Lyne A. G., 2005, MNRAS, 364, 1397
  * Kamchatnov (1982) Kamchatnov A. M., 1982, Sov. Phys. JETP, 55, 69
  * Kasen & Woosley (2009) Kasen D., Woosley S. E., 2009, ApJ, 703, 2205
  * Kaspi & Beloborodov (2017) Kaspi V. M., Beloborodov A. M., 2017, ARA&A, 55, 261
  * Keane & Kramer (2008) Keane E. F., Kramer M., 2008, MNRAS, 391, 2009
  * Lai & Goldreich (2000) Lai D., Goldreich P., 2000, ApJ, 535, 402
  * Lattimer & Prakash (2001) Lattimer J. M., Prakash M., 2001, ApJ, 550, 426
  * Lattimer & Schutz (2005) Lattimer J. M., Schutz B. F., 2005, ApJ, 629, 979
  * Lattimer & Swesty (1991) Lattimer J. M., Swesty D. F., 1991, Nuclear Phys. A, 535, 331
  * Lattimer & Yahil (1989) Lattimer J. M., Yahil A., 1989, ApJ, 340, 426
  * LeBlanc & Wilson (1970) LeBlanc J. M., Wilson J. R., 1970, ApJ, 161, 541
  * Lorimer & Kramer (2004) Lorimer D. R., Kramer M., 2004, Handbook of Pulsar Astronomy. Cambridge Observing Handbooks for Research Astronomers Vol. 4, Cambridge Universityx Press, Cambridge (UK)
  * Mabanta & Murphy (2018) Mabanta Q. A., Murphy J. W., 2018, ApJ, 856, 22
  * Makarenko et al. (2021) Makarenko E. I., Igoshev A. P., Kholtygin A. F., 2021, MNRAS, 504, 5813
  * Marek et al. (2006) Marek A., Dimmelmeier H., Janka H. T., Müller E., Buras R., 2006, A&A, 445, 273
  * Masada et al. (2015) Masada Y., Takiwaki T., Kotake K., 2015, ApJ, 798, L22
  * Matsumoto et al. (2020) Matsumoto J., Takiwaki T., Kotake K., Asahina Y., Takahashi H. R., 2020, MNRAS, 499, 4174
  * Matsumoto et al. (2022) Matsumoto J., Asahina Y., Takiwaki T., Kotake K., Takahashi H. R., 2022, arXiv e-prints, p. arXiv:2202.07967
  * Mestel (1999) Mestel L., 1999, Stellar magnetism (International Series of Monographs on Physics; 99). Clarendon Press Oxford
  * Mestel (2001) Mestel L., 2001, in Mathys G., Solanki S. K., Wickramasinghe D. T., eds, Astronomical Society of the Pacific Conference Series Vol. 248, Magnetic Fields Across the Hertzsprung-Russell Diagram. p. 3
  * Meynet et al. (2011) Meynet G., Eggenberger P., Maeder A., 2011, A&A, 525, L11
  * Miyoshi & Kusano (2005) Miyoshi T., Kusano K., 2005, Journal of Computational Physics, 208, 315
  * Moss (2001) Moss D., 2001, in Mathys G., Solanki S. K., Wickramasinghe D. T., eds, Astronomical Society of the Pacific Conference Series Vol. 248, Magnetic Fields Across the Hertzsprung-Russell Diagram. p. 305
  * 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 & Hillebrandt (1979) Müller E., Hillebrandt W., 1979, A&A, 80, 147
  * Müller & Janka (2015a) Müller B., Janka H. T., 2015a, MNRAS, 448, 2141
  * Müller & Janka (2015b) Müller B., Janka H. T., 2015b, MNRAS, 448, 2141
  * Müller & Varma (2020) Müller B., Varma V., 2020, MNRAS, 498, L109
  * Müller et al. (2016) Müller B., Viallet M., Heger A., Janka H.-T., 2016, ApJ, 833, 124
  * Müller et al. (2017) Müller B., Melson T., Heger A., Janka H.-T., 2017, MNRAS, 472, 491
  * Müller et al. (2019a) Müller B., et al., 2019a, MNRAS, 484, 3307
  * Müller et al. (2019b) Müller B., et al., 2019b, MNRAS, 484, 3307
  * Murphy et al. (2019) Murphy J. W., Mabanta Q., Dolence J. C., 2019, MNRAS, 489, 641
  * Noutsos et al. (2012) Noutsos A., Kramer M., Carr P., Johnston S., 2012, MNRAS, 423, 2736
  * Noutsos et al. (2013) Noutsos A., Schnitzeler D. H. F. M., Keane E. F., Kramer M., Johnston S., 2013, MNRAS, 430, 2281
  * Obergaulinger & Aloy (2021) Obergaulinger M., Aloy M. Á., 2021, MNRAS, 503, 4942
  * Obergaulinger et al. (2009) Obergaulinger M., Cerdá-Durán P., Müller E., Aloy M. A., 2009, A&A, 498, 241
  * Obergaulinger et al. (2014) Obergaulinger M., Janka H. T., Aloy M. A., 2014, MNRAS, 445, 3169
  * Olausen & Kaspi (2014) Olausen S. A., Kaspi V. M., 2014, ApJS, 212, 6
  * Perna et al. (2008) Perna R., Soria R., Pooley D., Stella L., 2008, MNRAS, 384, 1638
  * Petrovic et al. (2005) Petrovic J., Langer N., Yoon S. C., Heger A., 2005, A&A, 435, 247
  * Popov & Turolla (2012) Popov S. B., Turolla R., 2012, Ap&SS, 341, 457
  * Powell & Müller (2020) Powell J., Müller B., 2020, MNRAS, 494, 4665
  * Quentin & Tout (2018) Quentin L. G., Tout C. A., 2018, MNRAS, 477, 2298
  * Rampp & Janka (2000) Rampp M., Janka H. T., 2000, ApJ, 539, L33
  * Raynaud et al. (2020) Raynaud R., Guilet J., Janka H.-T., Gastine T., 2020, Science Advances, 6, eaay2732
  * Romani (1990) Romani R. W., 1990, Nature, 347, 741
  * Sawai et al. (2013) Sawai H., Yamada S., Suzuki H., 2013, ApJ, 770, L19
  * Scheck et al. (2008) Scheck L., Janka H. T., Foglizzo T., Kifonidis K., 2008, A&A, 477, 931
  * Schneider et al. (2019) Schneider F. R. N., Ohlmann S. T., Podsiadlowski P., Röpke F. K., Balbus S. A., Pakmor R., Springel V., 2019, Nature, 574, 211
  * Schneider et al. (2020) Schneider F. R. N., Ohlmann S. T., Podsiadlowski P., Röpke F. K., Balbus S. A., Pakmor R., 2020, MNRAS, 495, 2796
  * Schneider et al. (2021) Schneider F. R. N., Podsiadlowski P., Müller B., 2021, A&A, 645, A5
  * Schöller et al. (2017) Schöller M., et al., 2017, A&A, 599, A66
  * Spruit (2009) Spruit H. C., 2009, in Strassmeier K. G., Kosovichev A. G., Beckman J. E., eds, IAU Symposium Vol. 259, Cosmic Magnetic Fields: From Planets, to Stars and Galaxies. pp 61–74, doi:10.1017/S1743921309030075
  * Stockinger et al. (2020) Stockinger G., et al., 2020, MNRAS, 496, 2039
  * Suijs et al. (2008) Suijs M. P. L., Langer N., Poelarends A. J., Yoon S. C., Heger A., Herwig F., 2008, A&A, 481, L87
  * Takahashi & Yamada (2014) Takahashi K., Yamada S., 2014, ApJ, 794, 162
  * Thompson & Duncan (1993) Thompson C., Duncan R. C., 1993, ApJ, 408, 194
  * Torres-Forné et al. (2016) Torres-Forné A., Cerdá-Durán P., Pons J. A., Font J. A., 2016, MNRAS, 456, 3813
  * Tricco et al. (2016) Tricco T. S., Price D. J., Bate M. R., 2016, Journal of Computational Physics, 322, 326
  * Varma & Müller (2021) Varma V., Müller B., 2021, MNRAS, 504, 636
  * Varma et al. (2021) Varma V., Müller B., Obergaulinger M., 2021, Monthly Notices of the Royal Astronomical Society, 508, 6033
  * Vink & Kuiper (2006) Vink J., Kuiper L., 2006, MNRAS, 370, L14
  * Wieczorek et al. (2016) Wieczorek M., MMesch de Andrade E. S., Oshchepkov I., heroxbd 2016, SHTOOLS/SHTOOLS: Version 4.0, doi:10.5281/zenodo.206114, https://doi.org/10.5281/zenodo.206114
  * Wongwathanarat et al. (2013) Wongwathanarat A., Janka H.-T., Müller E., 2013, A&A, 552, A126
  * Woods (2008) Woods P. M., 2008, in Bassa C., Wang Z., Cumming A., Kaspi V. M., eds, American Institute of Physics Conference Series Vol. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More. pp 227–233, doi:10.1063/1.2900149
  * Woosley & Bloom (2006) Woosley S. E., Bloom J. S., 2006, ARA&A, 44, 507
  * Woosley & Heger (2006) Woosley S. E., Heger A., 2006, ApJ, 637, 914