text:
# The Magnetorotational Instability in Core Collapse Supernova Explosions

## 1 Introduction

There is accumulating evidence that core collapse supernovae are distinctly and significantly asymmetric. A number of supernova remnants show intrinsic “bilateral” structure (Dubner et al. 2002). Jet and counter jet structures have been mapped for Cas A in the optical (Fesen & Gunderson 1996; Fesen 2001; and references therein). Chandra X-ray Observatory (CXO) and XMM Newton observations of Cas A show that the jet and counter jet and associated structure are observable in the X-ray and that the intermediate mass elements are ejected in a roughly toroidal configuration (Hughes et al. 2000; Hwang et al. 2000; Willingale et al. 2002). HST observations of the debris of SN 1987A show that the ejecta are asymmetric with an axis that roughly aligns with the small axis of the rings (Pun et al. 2001; Wang et al. 2002). Spectropolarimetry shows that substantial asymmetry is ubiquitous in core-collapse supernovae, and that a significant fraction of core-collapse supernovae have a bi-polar structure (Wang et al. 1996, 2001). The strength of the asymmetry observed with polarimetry is higher (several %) in supernovae of Type Ib and Ic that represent exploding bare non-degenerate cores (Wang et al. 2001). The degree of asymmetry also rises as a function of time for Type II supernovae (from

<

<

\mathrel{\vbox{\hbox{$<$}\nointerlineskip\hbox{$\sim$}}} 1% to

>

>

\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}} 1%) as the ejecta expand and the photosphere recedes (Wang et al. 2001; Leonard et al. 2000, 2001). Both of these trends suggest that it is the core collapse mechanism itself that is responsible for the asymmetry.

Two possibilities are being actively explored to account for the observed asymmetries. One is associated with the rotational effect on convection (Fryer & Heger 2000) and another is due to the effect of jets (Khokhlov et al. 1999; Wheeler et al. 2000; Wheeler, Meier & Wilson 2002). The time scale of neutrino emission is short compared to the dynamical time of the overlying stellar mantle and envelope. Neutrino asymmetries can yield short-lived, impulsive effects (Shimizu, et al. 1994; Burrows & Hayes 1996; Lai et al. 2001), but there are questions of whether expansion and transverse pressure gradients will eliminate transient asymmetries before homologous expansion is achieved (Chevalier & Soker 1989). Rayleigh–Taylor and Richtmyer–Meshkov instabilities might produce “finger” asymmetries that are preserved, but it is unclear that such finer scale perturbations can reproduce the common feature of a single symmetry axis that is substantially independent of space and time (Wang et al. 2001). Jet calculations have established that non-relativistic axial jets of energy of order $`10^{51}`$ erg originating within the collapsed core can initiate a bi-polar asymmetric supernova explosion that is consistent with the spectropolarimetry (Khokhlov et al. 1999; Khokhlov & Höflich 2001; Höflich et al. 2001). The result is that heavy elements (e.g. O, Ca) are characteristically ejected in tori along the equator. Iron, silicon and other heavy elements in Cas A are distributed in this way (Hwang et al. 2000), and there is some evidence for this distribution in SN 1987A (Wang et al. 2002). Radioactive matter ejected in the jets can alter the ionization structure and hence the shape of the photosphere of the envelope even if the density structure is spherically symmetric (Höflich et al. 2001). This will generate a finite polarization, even though the density distribution is spherical and the jets are stopped deep within the star and may account for the early polarization observed in Type II supernovae (Leonard et al. 2000; Wang et al. 2001). If one of the pair of axial jets is somewhat stronger than the other, jets can, in principle, also account for pulsar runaway velocities that are parallel to the spin axis (Helfand et al. 2001, and references therein). While a combination neutrino–induced/jet–induced explosion may prove necessary for complete understanding of core-collapse explosions, jets of the strength computed by Khokhlov et al. (1999) are sufficient for supernova explosions. In this paper, we will explore the possible conditions that could lead to the formation of buoyant bi-polar MHD outflow.

Immediately after the discovery of pulsars there were suggestions that rotation and magnetic fields could be a significant factor in the explosion mechanism (Ostriker & Gunn 1971; Bisnovatyi-Kogan 1971; Bisnovatyi-Kogan & Ruzmaikin 1976; Kundt 1976). Typical dipole fields of $`10^{12}`$ G and rotation periods of several to several tens of milliseconds yield electrodynamic power of $`10^{4445}`$ erg s-1 that is insufficient to produce a strong explosion. The evidence for asymmetries and the possibility that bi-polar flows or jets can account for the observations suggest that this issue must be revisited. The fact that pulsars like those in the Crab and Vela remnants have jet-like protrusions (Weisskopf et al. 2000; Helfand et al. 2001) also encourages this line of thought. The present-day jets in young pulsars may be vestiges of much more powerful MHD jets that occurred when the pulsar was born. The transient values of the magnetic field and rotation could have greatly exceeded those observed today. Tapping that energy to power the explosion could be the very mechanism that results in the modest values of rotation and field the pulsars display after the ejecta have dispersed.

In the current context, one would like to know not only whether or not the rotation and magnetic field of a nascent neutron star can power the supernova explosion, but, specifically, whether or not bi-polar flows or jets form and whether or not they are sufficiently energetic to drive the explosion. Possible physical mechanisms for inducing axial flows, asymmetric supernovae, and related phenomena driven by magnetorotational effects were considered by Wheeler et al. (2000). The means of amplifying magnetic fields by linear wrapping associated with differential rotation in the neutron star and possibly by an $`\alpha `$ – $`\mathrm{\Omega }`$ dynamo were discussed. Whereas Wheeler et al. (2000) focused on the effect of the resulting net dipole field, Wheeler, Meier, & Wilson (2002) recognized that the toroidal field would be the dominant component of the magnetic field and explored the capacity of the toroidal field to directly generate axial jets by analogy with magneto-centrifugal models of jets in active galactic nuclei (Koide et al. 2000; and references therein). Wheeler, Meier, & Wilson (2002) found that the production of a strong toroidal field, substantially stronger than the $`10^{12}`$ G field of a pulsar, and strong axial jets driven by that field are possible.

A caveat to consider is that Wheeler, Meier, & Wilson (2002) again considered only amplification of the field by “wrapping,” a process that only increases the field linearly, and hence rather slowly in time. While conditions of very rapid rotation might exist that lead to a sufficiently rapid growth of the magnetic field, that is not guaranteed. In addition, reconnection might limit the field before it can be wrapped the thousands of times necessary to be interesting. Here we consider the effects of magnetic shearing, the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998), on the strongly shearing environment that must exist in a nascent neutron star. This instability is expected to lead to the rapid exponential growth of the magnetic field with characteristic time scale of order the rotational period. While this instability has been widely explored in the context of accretion disks, this is the first time it has been applied to core collapse. We will argue that this instability must inevitably occur in core collapse, that it is likely to be the dominant mechanism for the production of magnetic flux in the context of core collapse, and that it has the capacity to produce fields that are sufficiently strong to affect, if not cause, the explosion.

In §2 we describe our assumptions by which we treat the angular momentum and field amplification in the context of a numerical collapse calculation. In §3, we present the results of these calculations. Section 4 presents a discussion and conclusions.

## 2 Method

The collapse of a model iron core of a 15 M progenitor is simulated with a one-dimensional flux-limited diffusion code (Myra et al., 1987) that computes radial profiles for various quantities. The core collapse code does not include rotation, convection, or magnetic fields. To treat rotation of the iron core, we attribute various initial rotational profiles to the initial density profile. The evolution of the angular velocity profile $`\mathrm{\Omega }(r)`$ is then computed using the radial density profiles produced by the core collapse code assuming that the specific angular momentum of a given shell is constant (§2.1). The magnetic field is obtained using the $`\rho (r)`$, $`\mathrm{\Omega }(r)`$, and $`d\mathrm{\Omega }/dr`$ profiles according to the theory of the MRI (§2.2). The dynamical feedback of the magnetic field through mechanisms such as magnetic pressure, torque, and hoop stress are neglected here. The MHD luminosity is estimated using the magnetic field calculated (§2.3).

In this calculation, the initial iron core of 1.5 M forms a collapsing homologous core of mass 0.69 M. The bounce occurs about 166 ms after the onset of collapse and the shock forms at the boundary of the initial homologous core at a radius of 13.8 km. After bounce, the outer boundary of the protoneutron star (PNS) extends to about 80 – 90 km. The bounce shock reaches 500 km at about 120 ms after bounce where it stalls and remains at about this radius for the duration of the calculation. Between the boundary of the PNS and the shock radius, there is a region of transition where shock–heated material settles onto the PNS. By 210 ms after bounce, the PNS develops a core in hydrostatic equilibrium (v $``$ 100 km s-1) with radius 9 km and mass 0.2 M. The hydrostatic core grows out to 24.4 km and 0.69 M by 350 ms and stays constant to the end of the simulation. The hydrostatic core mass of the PNS at this time encompasses all the initial homologous core mass, 0.69 M. Although the original homologous core and the hydrostatic core have the same mass, the size of the PNS core is larger because the density decreases after bounce.

### 2.1 Angular Velocity Profile

The evolution of massive stars can be significantly affected by their rotation (Eddington, 1925; Kippenhahn et al., 1970; Endal & Sofia, 1976, 1978; Kippenhahn & Weigert, 1989; Meynet & Maeder, 2000; Maeder & Meynet, 2000), metallicity (Meynet et al., 1994), and mass loss (Chiosi & Maeder, 1986). Massive stars that are the progenitors of core collapse supernovae are known to be rapid rotators as main sequence stars with mean $`v\mathrm{sin}i`$ ranging from 100 km s-1 to 300 km s-1 (Slettebak, 1956; Conti & Ebbets, 1977; Fukuda, 1982; Penny, 1996; Howarth et al., 1997). Maeder et al. (1999) showed that the number ratio of Be to B-type stars increases as the local metallicity decreases, suggesting that massive stars rotate faster in lower metallicity environments. According to the parameterization of mass loss rates employed by Heger et al. (2000), the mass loss rates are enhanced by faster rotation and higher metallicity. Evolved stars are also observed to have higher mass loss rates. These effects are interdependent. Although mass loss significantly influences the evolutionary tracks in the HR diagram, it does not affect the evolution of the central core since after the He burning phase the core evolves independently from the envelope (Chiosi & Maeder, 1986). Effects of the magnetic field on stellar evolution (Heger & Woosley, 2002) may influence the timing of the decoupling between the core and the envelope (Fryer & Heger, 2000), which could affect the final angular momentum of the core. It is important to understand the evolution in both the envelope and the core of the progenitor star up to the precollapse stage since the outcome of core collapse may depend significantly on the precollapse conditions. For our study, the structure of the precollapse iron core, including the angular momentum distribution, is of particular significance.

The first calculations of the redistribution of angular momentum in the evolution of rotating massive stars were carried by Endal & Sofia (1978), who obtained angular velocity profiles up to the ignition of carbon burning. The evolution of the angular momentum distribution of massive stars up to the onset of core collapse was recently simulated by Heger et al. (2000) for supernova progenitors (8$`M_{\mathrm{}}`$$`25M_{\mathrm{}}`$) with solar composition. This study included the effects of rotation and angular momentum transport with mass loss. They found that the final specific angular momentum distribution prior to core collapse had no dependence on the initial surface velocity and only a weak dependence on the initial mass; larger initial masses attain smaller specific angular momenta prior to collapse because higher mass loss rates carry away more angular momentum. Heger et al. (2000) also found that the $`\mathrm{\Omega }(r)`$ profile in the precollapse core is not a smooth function of radius, but consists of step functions at fossil convective shell boundaries. The $`15M_{\mathrm{}}`$ model of Heger et al. (2000) attains an angular velocity of 10 rad s-1 (see their Fig. 8) in the center of the iron core at the precollapse stage. These simulations did not include the effects of a magnetic field. It is possible that the iron core rotates more slowly if the effect of magnetic braking is included (Spruit & Phinney, 1998; Heger & Woosley, 2002). Fryer & Heger (2000) studied the rotational effects on pure hydrodynamic core collapse explosions with initial velocity profiles obtained by Heger et al. (2000) with a central rotational velocity of 4 rad s-1. We adopt the initial angular velocity profile, $`\mathrm{\Omega }_0(r)`$, of Fryer & Heger (2000) (hereafter called the FH profile) as one case to study.

An analytical form of the $`\mathrm{\Omega }_0(r)`$ profile had been adopted in some previous work (Mönchmeyer & Müller, 1989; Yamada & Sato, 1994; Fryer & Heger, 2000). One representation of the profile takes the form (Mönchmeyer & Müller, 1989):

$$\mathrm{\Omega }_0(r)=\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}\frac{R^2}{r^2+R^2},$$
(1)
where $`R=10^8`$ cm (hereafter called the MM profile). The FH and MM profiles with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=4.0`$ rad s-1 resemble each other (Fryer & Heger, 2000), except that the MM profile is a continuous function, and the FH profile is broken by steps. The pioneering work on core collapse supernova explosions of rotating stars by LeBlanc & Wilson (1970) employed an initial solid body rotation with 0.7 rad s-1.

We adopt the following three profiles for this study; solid body, MM, and FH profiles, as illustrated in Fig. 1. These profiles are scaled by the central value of $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$ that we choose to range from 0.2 to 10 rad s-1. As the collapse progresses, the PNS will rotate near breakup speed for sufficiently high $`\mathrm{\Omega }_0(r)`$. The solid body case is most susceptible to this problem. For this work, we monitor the rotation and omit cases that are not approximately self-consistent with our spherically-symmetric treatment (although not necessarily unphysical). The adopted profiles that we study in detail have small enough angular momentum that little departure from spherical geometry will occur.

Once an $`\mathrm{\Omega }_0(r)`$ profile is adopted, the $`\mathrm{\Omega }(r)`$ profile is computed using the radial density profiles produced by the core collapse code described above. We assume that the specific angular momentum is conserved between consecutive time steps as follows:

$`\mathrm{\Omega }^{}(r^{})`$ $`=`$ $`\left({\displaystyle \frac{r}{r^{}}}\right)^2\mathrm{\Omega }(r)`$ (2)
$`=`$ $`\left({\displaystyle \frac{\rho ^{}}{\rho }}\right)^{\frac{2}{3}}\mathrm{\Omega }(r),`$ (3)
where the primed quantities represent the future time step (the appropriate expressions for finite shells in difference form was used in the code as applied below). It is inevitable that the collapsing core spins up and generates strong differential rotation for very general choices of the $`\mathrm{\Omega }_0(r)`$ profile, since the inner regions collapse larger relative distances than the outer regions. A strong shear must form at the boundary of the PNS. At bounce, the original homologous core has a positive gradient in $`\mathrm{\Omega }(r)`$, but by 210 ms after bounce, this gradient has nearly disappeared and after that the gradient is monotonically negative outward (see below; Figs. 7, 11).

Ruderman et al. (2000) noted that the collapse of a white dwarf to a PNS gives a positive $`\mathrm{\Omega }(r)`$ gradient since the relativistic degenerate core of the white dwarf has a steeper density profile than the PNS. The PNS will thus be relatively more compact for a given central density. There are two important differences in the current calculations. The most critical is that we are not considering a core collapsing in isolation as for the accretion-induced collapse scenario. Rather, the PNS forms within the massive star collapse ambience, and the PNS must be strongly differentially rotating with respect to the still-infalling matter. This will generate a strong shear at the boundary of the PNS that would not pertain to a collapsing isolated white dwarf. Another, more subtle, difference is the equation of state. Here we are considering the collapse of a partially degenerate iron core to a PNS. The equation of state is not as different as for the highly relativistic white dwarf collapsing to non-relativistic neutron star. An examination of the relative density profiles shown in Fig. 2 reveals that the homologous core at bounce is somewhat less centrally condensed than the original iron core; this results in the temporary positive gradient in $`\mathrm{\Omega }`$ within the homologous core as shown below in Figs. 7, 11. About 50 ms after bounce, the density profile is nearly identical to that of the initial iron core, giving the nearly flat rotation profile shown below. After that, the density profile becomes somewhat more centrally condensed than the original iron core and the rotation profile decreases monotonically outward even deep within the PNS.

### 2.2 Magnetic Field

The MRI generates turbulence in a magnetized rotating fluid body that amplifies the magnetic field and transfers angular momentum. In this paper we attempt, for the first time, to apply the physics of the MRI to the core collapse supernova environment where the MRI criterion is broadly satisfied. The MRI should pertain in this environment and amplify the magnetic field exponentially and perhaps, in turn, power MHD bi-polar flow or jets. Key questions are the amplitude of the magnetic field and the effect on the dynamics. Those two questions are deeply related, but with the limitations of this proof-of-principle study in mind, we will estimate the strength of the magnetic field due to the MRI during the core collapse process and postpone a study of feedback effects to a later work.

Ignoring entropy gradients for the moment, the condition for the instability of the slow magnetosonic waves in a magnetized, differentially rotating plasma is (Balbus & Hawley, 1991, 1998):

$$\frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}+(𝐤𝐯_\mathrm{A})^2<0,$$
(4)
where

$$v_\mathrm{A}=\frac{B}{\sqrt{4\pi \rho }}$$
(5)
is the Alfvén velocity. When the magnetic field is very small, and/or the wavelength is very long, $`(𝐤𝐯_\mathrm{A})^2`$ is negligible, and the instability criterion for the MRI is simply that the angular velocity gradient be negative (Balbus & Hawley, 1991, 1998), i.e.:

$$\frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}<0.$$
(6)
The growth of the magnetic field associated with the MRI is exponential with characteristic time scale of order $`\mathrm{\Omega }^1`$ since $`\tau 2\pi /|\omega |2\pi |d\mathrm{\Omega }^2/d\mathrm{ln}r|^{1/2}`$. The time scale for the maximum growing mode is given by (Balbus & Hawley, 1998):

$$\tau _{\mathrm{max}}=4\pi \left|\frac{d\mathrm{\Omega }}{d\mathrm{ln}r}\right|^1.$$
(7)
We thus expect the MRI to dominate any process such as “wrapping of field lines” (see Wheeler et al. 2000 and references therein) that only grows linearly in time, even if on about the same time scale. The MRI will also operate under conditions of moderate rotation that are not sufficient to compete with the PNS convective time scales to drive the sort of $`\alpha `$ – $`\mathrm{\Omega }`$ dynamo invoked by, e.g., Duncan & Thompson (1992). The resulting unstable flow is expected to become non-linear, develop turbulence, and drive a dynamo that amplifies and sustains the field. The field will grow because of the MRI dynamo action until it reaches a saturation field limit. The details of the dynamo action due to MHD turbulence are beyond the scope of this paper, and we ignore the effect of the turbulence on the hydrodynamics.

#### 2.2.1 Saturation Fields

An order of magnitude estimate for the saturation field can be obtained by equating the shearing length scale $`l_{\mathrm{shear}}dr/d\mathrm{ln}\mathrm{\Omega }`$ to the characteristic mode scale $`ł_{\mathrm{mode}}v_\mathrm{A}(d\mathrm{\Omega }/d\mathrm{ln}r)^1`$. The resulting saturation magnetic field is given by:

$$B_{\mathrm{sat}}^24\pi \rho r^2\mathrm{\Omega }^2.$$
(8)
This is the same result as obtained by setting the Alfvén velocity equal to the local rotational velocity, $`v_\mathrm{A}=r\mathrm{\Omega }`$. In the subsequent discussion, we will adopt the field in eq. (8) as the fiducial saturation field. Note that this expression does not depend on the shear, although the field growth in this case will depend on the shear through the instability criterion. Another simple estimate of the saturation field is to set the mode length equal to the local radius, on the grounds that the mode cannot be larger in wavelength than the physical size of the region in which it grows. The resulting saturation field is:

$$B_{\mathrm{rad}}^24\pi \rho r^2\mathrm{\Omega }^2\left(\frac{d\mathrm{ln}\mathrm{\Omega }}{d\mathrm{ln}r}\right)^2.$$
(9)
This differs from the previous estimate only by the logarithmic gradient of $`\mathrm{\Omega }`$, a factor normally of order unity.

The empirical value of the saturation field obtained by the numerical simulation of Hawley et al. (1996) is:

$`B_{\mathrm{sim}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{\rho }{\pi }}}r\mathrm{\Omega }`$ (10)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}B_{\mathrm{sat}}.`$
This saturation field is achieved after turbulence is fully established, which takes about 20 rotations following the initial exponential growth (Hawley et al., 1996). The field of eq. (10) is another possible saturation limit. Note that for conditions of rotation at much less than Keplerian, to which we restrict the current analysis, these saturation fields and that to be described below, represent fields that are much less than the equipartition field, for which $`B^2/8\pi `$ is comparable to the ambient pressure, i.e. for the current calculations, $`c_\mathrm{s}r\mathrm{\Omega }v_\mathrm{A}`$ (Fig. 14).

Another means of estimating the saturation field considers the MRI mode that grows at the maximum rate (eq. (7)). When a vertical field exists, the maximum unstable growing mode (Balbus & Hawley, 1998) implies a saturation field of:

$$B_{\mathrm{max}}^2=4\pi \rho \lambda _{\mathrm{max}}^2\mathrm{\Omega }^2\left[\frac{1}{8\pi ^2}\left(1+\frac{1}{8}\frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)\frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right],$$
(11)
where $`\lambda _{\mathrm{max}}`$ is the wavelength of the maximum growing mode. Once again, we argue that this fastest growing mode must be contained in the physical region in which growth occurs, so the wavelength is not allowed to exceed the local radius $`r`$. With $`\lambda _{\mathrm{max}}=r`$, eq. (11) becomes:

$$B_{\mathrm{max}}^2=B_{\mathrm{sat}}^2\left[\frac{1}{8\pi ^2}\left(1+\frac{1}{8}\frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)\frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right].$$
(12)
This expression for the saturation field depends on the shear explicitly as well as indirectly through the stability criterion. Yet another variation would be to set the wavelength of the maximum growing mode equal to the shear length, but all these variations just yield differences of order $`d\mathrm{ln}\mathrm{\Omega }/d\mathrm{ln}r`$. Note that for the maximum growing mode the expression for $`B_{\mathrm{max}}^2`$ acquires a negative value when

$`{\displaystyle \frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}}`$ $`<`$ $`8\mathrm{\Omega }^2\text{ or },`$
$`\kappa ^2`$ $`<`$ $`4\mathrm{\Omega }^2<0,`$ (13)
where $`\kappa `$ is the epicyclic frequency:

$$\kappa ^2=\frac{1}{r^3}\frac{d(r^4\mathrm{\Omega }^2)}{dr}=4\mathrm{\Omega }^2+\frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}.$$
(14)
When eq. (2.2.1) is true, the epicyclic motion dominates over the MRI and prevents growth of the perturbation. In practice, the gradient of $`\mathrm{\Omega }`$ may be reduced by mixing due to the epicyclic motion, and the MRI may eventually be active in a region in which it was at first suppressed by a strong negative gradient of $`\mathrm{\Omega }(r)`$. In this work, we simply turned off field amplification when this condition arose, thus, perhaps, minimizing the actual effect of the MRI.

In many cases, including our collapsing core, the background gas is unstable to convection due to an entropy gradient. The previous discussion can be modified to take this effect into account. When differentially rotating gas is unstable to the MRI and stable to convection, a perturbed fluid element experiences a stabilizing force due to the convective stability. Therefore, convective stability tends to stabilize the MRI. The modified instability criteria for the MRI with an entropy gradient are (Balbus & Hawley, 1991, 1998):

$`N^2+{\displaystyle \frac{\mathrm{\Omega }^2}{\mathrm{ln}\varpi }}<0,`$ (15)
$`\left({\displaystyle \frac{P}{z}}\right)\left({\displaystyle \frac{\mathrm{\Omega }^2}{\omega }}{\displaystyle \frac{s}{z}}{\displaystyle \frac{\mathrm{\Omega }^2}{z}}{\displaystyle \frac{s}{\omega }}\right)<0,`$ (16)
where $`\varpi `$ and $`z`$ are in cylindrical coordinates, s is entropy, and N is the Brunt-Väisälä frequency:

$$N^2=\frac{g\delta }{C_p}s$$
(17)
where $`\delta =\left(\frac{\mathrm{ln}\rho }{lnT}\right)_P`$, $`C_P`$ is the specific heat at constant pressure.

For our spherical model, the LHS of eq. (16) is zero. Note that the pure hydrodynamic counterparts to this instability, the Høiland criteria, can be obtained by simply replacing angular velocity with angular momentum, in eq. (15) and eq. (16), but the Høiland criteria cannot be obtained by taking the limit of zero magnetic field since eq. (15) and eq. (16) do not depend on field strength (Balbus & Hawley, 1998). The spherically symmetric collapse models used in this work develop regions that are unstable to convection, corresponding to a negative value of $`N^2`$, but are not able to compute the development of convection. The breakout of convection, although complex and time-dependent, would tend to lower the “superadiabatic” gradient and hence to reduce the value of $`|N^2|`$ compared to the values we compute here from eq. (17). For this reason, the values of $`N^2`$ we compute and employ in eq. (15) and other expressions may be exaggerated in convectively unstable regions, with implications for our calculations that will be described below.

For a rotating magnetized object with axial symmetry, the maximum growing mode in the presence of entropy gradients corresponding to $`\lambda _{\mathrm{max}}=r`$ is:

$$B_{\mathrm{max},\mathrm{en}}^2=B_{\mathrm{sat}}^2\left[\frac{\eta ^21}{4\pi ^2}+\frac{1}{8\pi ^2}\left(\eta +\frac{1}{8}\left(\xi \frac{N^2}{\mathrm{\Omega }^2}+\eta \frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)\right)\left(\xi \frac{N^2}{\mathrm{\Omega }^2}+\eta \frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)\right],$$
(18)
and the growth time scale is given by:

$$\tau _{\mathrm{max},\mathrm{en}}=2\pi \left|(\eta ^22\eta +1)\mathrm{\Omega }^2+\frac{\eta 1}{2}\left(\xi N^2+\eta \frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)+\frac{1}{16\mathrm{\Omega }^2}\left(\xi N^2+\eta \frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)^2\right|^{\frac{1}{2}},$$
(19)
where:

$`\xi ^2`$ $`=`$ $`(1\mathrm{sin}2\theta )^2,`$ (20)
$`\eta ^2`$ $`=`$ $`\mathrm{sin}^2\theta (1\mathrm{sin}2\theta ).`$ (21)
When the polar angle $`\theta =\pi /2`$, we obtain the instability criterion in the equatorial plane:

$$N^2+\frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}<0,$$
(22)
and the saturation field of the maximum growing mode becomes:

$$B_{\mathrm{max},\mathrm{en}}^2=B_{\mathrm{sat}}^2\left[\frac{1}{8\pi ^2}\left(1+\frac{1}{8}\left(\frac{N^2}{\mathrm{\Omega }^2}+\frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)\right)\left(\frac{N^2}{\mathrm{\Omega }^2}+\frac{d\mathrm{ln}\mathrm{\Omega }^2}{d\mathrm{ln}r}\right)\right],$$
(23)
which grows exponentially with the time scale

$$\tau _{\mathrm{max},\mathrm{en}}=4\pi \left|\frac{N^2}{2\mathrm{\Omega }}+\frac{d\mathrm{\Omega }}{d\mathrm{ln}r}\right|^1.$$
(24)
The field growth for the case of the saturation field in eq. (23) depends explicitly on the shear and Brunt-Väisälä frequency and depends indirectly on those quantities through the instability criterion and growth time. For $`N^2=0`$, the saturation field of the maximum growing mode and its growth time scale reduce to those of the no–entropy–gradient case (eq. (7), (12)).

From eq. (23), $`B_{\mathrm{max},\mathrm{en}}^2`$ becomes negative due to competition with the epicyclic frequency when

$$N^2+\frac{d\mathrm{\Omega }^2}{d\mathrm{ln}r}<8\mathrm{\Omega }^2.$$
(25)
Convective instability ($`N^2<0`$) assists the MRI to destabilize the flow according to the instability criterion in eq. (22). For the maximum growing mode, however, even when the instability criterion is met, the growth of the mode is inhibited when the magnitude of $`N^2`$ is negative and large because of the condition in eq. (25). In our calculations, the field with saturation limit $`B_{\mathrm{sat}}`$ grows whenever the instability criterion is satisfied. In this case, the MRI is inhibited in regions of convective stability and promoted in regions of convective instability. The field with saturation limit $`B_{\mathrm{max},\mathrm{en}}`$ grows when the instability criterion is satisfied, but epicyclic motions do not dominate, i. e., when $`8\mathrm{\Omega }^2<N^2+d\mathrm{\Omega }^2/d\mathrm{ln}r<0`$. Note that our estimates of $`N^2`$ in convectively unstable regions may exaggerate the effects of convection in eqs. (22), (23), (24), and (25). In particular, an exessively large negative value of $`N^2`$ in eq. (22) could promote the instability to the MRI artificially. This could affect either prescription for saturation fields. On the other hand, an excessively large negative value for $`N^2`$ in eq. (25) will inhibit the MRI from operating in conditions where it might actually do so by exaggerating the epicyclic effects. In this sense, our estimates of the activity of the MRI may be conservative.

In the numerical work to be presented below, we considered two prescriptions for the saturation field in the equatorial plane, the fiducial field $`B_{\mathrm{sat}}`$ corresponding to $`v_\mathrm{A}=r\mathrm{\Omega }`$ (eq. (8)), and $`B_{\mathrm{max},\mathrm{en}}`$, which corresponds to the maximum growing mode with the effects of the entropy gradient included (eq. (23)). We apply the instability criterion, eq. (22), and the growth time of eq. (24) to both saturation field prescriptions and the epicyclic limit, eq. (25), to the prescription of eq. (23). Cases corresponding to the numerical saturation limit (eq. (10)) can be obtained by simply dividing the resulting field profile in the fiducial saturation field case by 2$`\pi `$.

Any initial vertical field is turned into a radial and an azimuthal field by the MRI, leading to channel flows in 2D simulations and spiral streaming in 3D simulations (Balbus & Hawley, 1998). Other initial field configurations are also known to be unstable (Vishniac & Diamond, 1992; Balbus & Hawley, 1998). The flow induced by the MRI quickly develops turbulence. Hawley et al. (1996) found in their study of dynamo action due to the MRI that the saturation magnetic energy has a very weak dependency on the initial field strength and configuration as long as there is some vertical component to the original seed field. Both the field strength and the configuration of the magnetic field in the precollapse iron core are unknown. Here, we assume that in the precollapse iron core the magnetic field is much weaker than the equipartition field so that the MRI is strongly unstable, resulting in exponential growth of the field, and that there is some vertical component to the field so that the saturation fields we adopt are representative of the strongest limit expected from the MRI. These preliminary calculations represent a reasonable approximation of the effect of the MRI on growth of the magnetic field in the rotating and collapsing cores that precede the supernovae explosions.

#### 2.2.2 Time Dependent Calculation

We applied the time scale, $`\tau `$, of eq. (24) so that the magnetic field amplified by the MRI achieves the asymptotic saturation value $`B_{\mathrm{asym}}`$ only on this time scale. To account for this finite growth time to reach saturation, we take the growth of the magnetic field $`\mathrm{\Delta }B`$ during interval $`\mathrm{\Delta }t`$ to be:

$$\mathrm{\Delta }B=(1e^{\frac{\mathrm{\Delta }t}{\tau }})(B_{\mathrm{asym}}B_{\mathrm{old}}),$$
(26)
where Basym is calculated using eq. (8) or eq. (23) and $`B_{\mathrm{old}}`$ is the magnetic field achieved in the previous time step. Note that the growth time depends on $`N^2`$ and may be artificially short in strongly convectively unstable regions. The magnetic field for the current time step is then computed as:

$$B=B_{\mathrm{old}}+\mathrm{\Delta }B,$$
(27)
where $`\mathrm{\Delta }`$B goes to zero as the field asymptotically approaches the saturation limit.

### 2.3 MHD Power

MHD jets are common in systems with a central object that is accreting matter with angular momentum and a magnetic field (Meier et al., 2001). Examples of such systems are the engines of active galactic nuclei, micro quasars, gamma–ray bursts (collapsars), and young stellar objects. In the collapse of a rotating iron core, the environment for MHD jets to form is satisfied by a central PNS with infalling matter with angular momentum and a magnetic field. Previous studies have shown that buoyant MHD outflows form in the collapse of a magnetized rotating iron core (LeBlanc & Wilson, 1970; Symbalisty, 1984). Wheeler, Meier, & Wilson (2002) discussed jet formation in the context of magnetorotational core collapse supernovae, but based only on the linear wrapping of field lines. Although jet formation is beyond the scope of this paper, we can estimate the power produced by the magnetic field generated by the MRI that could drive bi-polar flows or jets.

The characteristic power of non-relativistic MHD outflow is given by Blandford & Payne (1982); see also (Meier, 1999; Wheeler, Meier, & Wilson, 2002):

$$L_{\mathrm{MHD}}=\frac{B^2r^3\mathrm{\Omega }}{2}.$$
(28)
The outflow carries energy, angular momentum, and mass. The ejected mass is contained in a cone with a fractional solid angle $`f_{\mathrm{\Delta }\mathrm{\Omega }}=\mathrm{\Delta }\mathrm{\Omega }/4\pi `$, where $`\mathrm{\Delta }\mathrm{\Omega }`$ is the solid angle of the cone about the rotational axis of the PNS. The critical power required for the plasma in the cone to escape from the gravitational potential during a free fall time is given by Meier (1999) and Wheeler, Meier, & Wilson (2002),

$$L_{\mathrm{crit}}=\frac{E_{\mathrm{esc}}}{\tau _{\mathrm{ff}}}=\frac{4\pi }{3}\rho r^2f_{\mathrm{\Delta }\mathrm{\Omega }}\left(\frac{GM}{r}\right)^{\frac{3}{2}}.$$
(29)
When $`L_{\mathrm{MHD}}`$ exceeds $`L_{\mathrm{crit}}`$, the plasma can be accelerated to escape velocity or higher. The ratio of $`L_{\mathrm{MHD}}`$ to $`L_{\mathrm{crit}}`$ is expressed as a dimensionless parameter,

$$\nu \left(\frac{L_{\mathrm{MHD}}}{L_{\mathrm{crit}}}\right)^{\frac{1}{2}}=f_{\mathrm{\Delta }\mathrm{\Omega }}^{1/2}\left(\frac{3\mathrm{\Omega }}{\mathrm{\Omega }_{\mathrm{kep}}}\right)^{\frac{1}{2}}\left(\frac{v_\mathrm{A}}{v_{\mathrm{esc}}}\right).$$
(30)
Numerical calculations have shown that when $`\nu `$ is less than 1, the dynamics correspond to a Blandford-Payne-type outflow (Meier et al., 1997; Meier, 1999), but when $`\nu `$ is greater than 1, the character of the flow changes to that of a faster and more tightly collimated jet. For our case, we expect a non-relativistic Blandford-Payne-type outflow, which has a broad opening at the base of the jet. For simplicity, we assume $`f_{\mathrm{\Delta }\mathrm{\Omega }}=1`$.

## 3 Results

Because the code used for these calculations employs spherical geometry and cannot include rotational effects explicitly, our study is limited to the “slow rotation” case, where we can assume spherical geometry to good approximation. The initial rotational velocities in the iron core are so small that a significant departure from spherical geometry is not expected except for the solid rotation cases with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}>2.0`$ rad s-1. As the core collapses and density increases, however, it is inevitable that the rotational velocity increases (eq. (2), (3)) while generating strong differential rotation, even for initially slow rotation. Fig. 3 shows the time evolution of the maximum value of the ratio of $`\mathrm{\Omega }(r)`$ to the local Keplerian velocity as a function of $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$. For all initial rotational profiles, the ratio is linear with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$ at a given time. At 46 ms after bounce, some material is already rotating above break-up speed if Ω0_c

>

>

subscriptΩ0_cabsent\Omega_{\rm{0\_c}}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}} 5.5 rad s-1 for the MM and FH initial profiles and if Ω0_c

>

>

subscriptΩ0_cabsent\Omega_{\rm{0\_c}}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}} 3.0 rad s-1 for the solid initial profile. As time elapses, the solid profile case acquires break-up velocity for lower $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$; $`\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{kep}}`$ 1 for $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$ 1.2 at 387 ms after bounce. Therefore for Ω0_c

>

1.2

>

subscriptΩ0_c1.2\Omega_{\rm{0\_c}}\mathrel{\vbox{\hbox{$>$}\nointerlineskip\hbox{$\sim$}}}1.2, the envelope is expected to be severely distorted, possibly creating a disk-like structure. The initial differential rotation cases with the MM and FH profiles behave differently from the solid profile. The maximum value of the ratio $`\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{kep}}`$ decreases with time. To be self-consistent, we will analyze only the initial solid profile with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}0.2`$ rad s-1, and the MM and FH profiles with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}1.0`$ rad s-1. We stress that more rapid rotation in the initial models is not necessarily unphysical, only that we cannot treat such cases here. By selecting slow rotation cases, the structure is assumed to be stable to nonaxisymmetric instabilities.

As expected, the shear is the strongest at the boundary of the initial homologous core as seen in Fig. 4. At bounce the shear is positive inside of the initial homologous core, and the region is stable against the MRI though convective instability can destabilize the structure. The solid body profile possesses a similar shear profile to the FH and MM profiles, but the magnitude of the shear is smaller since this case started with smaller $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$. The spikes in the FH profile are features due to the discontinuities in $`\mathrm{\Omega }_0`$. By the end of the calculations, the shear is negative in most regions, and the magnitude of the shear is diminished for all cases.

Fig. 5 gives the evolution of the angular momentum per unit mass, $`j`$, for representative initial models. The structure is everywhere stable to the Rayleigh criterion except at the composition discontinuities in the FH model. The solid body model with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=0.2`$ rad s-1 has smaller $`j`$ in the initial model and in the post-collapse inner layers than the MM or FH cases with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=1.0`$ rad s-1 but has comparable values of $`j`$ in the outer regions at later times.

The collapsing core has convectively stable regions with $`N^2>0`$ and unstable regions with negative entropy gradient or $`N^2<0`$ (Fig. 6). In general, the PNS (less than 80 – 90 km) is convectively stable, and the region between the boundary of the PNS and the stalled shock is convectively unstable. After bounce, the convective stability in the original homologous core is sufficient to damp the MRI within the PNS in accord with eq. (22). For this reason, the field within 10 km presented below was amplified before bounce and then frozen in. When ambient gas is unstable to convection, the MRI is enhanced (eq. (22), (23), (24)), unless the epicyclic frequency is too high (eq. (25)) for the case with $`B_{\mathrm{max},\mathrm{en}}`$. Note that this one-dimensional code does not actually include convection, so the Brunt-Väisälä frequency we compute here is a formal quantity and likely to be higher than the eddy overturn frequency of fully developed convection. For this reason, the MRI may be artificially activated in convectively unstable regions with positive gradient of $`\mathrm{\Omega }`$. In practice, this does not appear to be a significant factor since the structure is broadly unstable to the MRI anyway. The growth time in convectively unstable regions could be artificially short. The most significant effect of our prescription for the Brunt-Väisälä frequency may be to exaggerate the epicyclic effects through eq. (25), especially for regions with small rotational velocity. The result is to shut off the MRI in regions with large negative $`N^2`$ and small $`\mathrm{\Omega }`$. This has a substantial effect on the case with saturation field corresponding to the maximum growing mode for which eq. (25) applies, but does not affect the case with the fiducial saturation field.

### 3.1 Initial Differential Rotation Case

Both of our initial differential rotation profiles, FH and MM, produce strong differential rotation (Fig. 7) due to collapse of the core with the strongest shear at the boundary of the initial homologous core (Fig. 4) which later becomes the hydrostatic PNS core. At bounce, the initial homologous core has positive $`\mathrm{\Omega }(r)`$ gradient, but the bounce of the initial homologous core eventually flattens the profile and makes the gradient negative. The differences between the FH and MM profiles are that the former preserves discontinuities present in the initial profile in subsequent rotational profiles, as clearly seen in the plot of $`\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{kep}}`$ in Fig. 7, and that the peak value of $`\mathrm{\Omega }(r)`$ at bounce in the FH profile is higher than that of the MM profile even though they have the same $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=1.0`$ rad s-1. This higher maximum value of $`\mathrm{\Omega }`$ at bounce is because the first step in the FH profile, which extends to 5.3 $`\times 10^7`$ cm in the initial iron core, gives a larger $`\mathrm{\Omega }_0`$ at finite radius than the MM profile (Fig. 1). Since both have the same $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$, the central rotational velocities stay the same: 1180 rad s-1 at bounce and 890 rad s-1 387 ms after bounce, which gives periods of 5.3 ms and 7.1 ms respectively. Even a relatively modest value of $`\mathrm{\Omega }_0`$ gives a very rapidly rotating PNS and hence strong differential rotation with respect to the infalling matter.

Because we have selected “slow rotator” conditions, the rotation remains sub–Keplerian (Fig. 7). At bounce, the peak of $`\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{kep}}`$ is at the boundary of the initial homologous core for both profiles, but the peak for the FH case has a higher value reflecting faster rotation for the reason mentioned above. At later times, however, the peaks move to the second hump which is located inside the stalled shock. This hump is at the same location as a maximum in entropy (Fig. 6) which is caused by shocked material with higher density.

For the given initial rotational profiles, the fiducial field $`B_{\mathrm{sat}}`$ and that corresponding to the maximum growing mode with entropy gradient $`B_{\mathrm{max},\mathrm{en}}`$ (eq. (23)) are amplified exponentially according to eq. (26) and (27). The resulting magnetic field and the ratio $`\beta ^1P_{\mathrm{mag}}/P_{\mathrm{gas}}`$ (where $`\beta `$ is the conventional $`\beta `$ in plasma physics) for $`B_{\mathrm{sat}}`$ are presented in Fig. 8 and for $`B_{\mathrm{max},\mathrm{en}}`$ in Fig. 9. For both magnetic field prescriptions, there is no significant difference between the MM and FH profiles, except that FH achieves a slightly higher magnetic field and thus larger $`\beta ^1`$ because of the higher rotational velocity at the boundary of the initial homologous core where the shear is the strongest.

Considerable differences occur between the $`B_{\mathrm{sat}}`$ and $`B_{\mathrm{max},\mathrm{en}}`$ cases. These differences are due to the fact that the $`B_{\mathrm{sat}}`$ case depends directly only on the azimuthal angular velocity and not the shear. This case depends on the radial gradients reflected in the shear and the Brunt-Väisälä frequency only indirectly through the instability criterion, eq. (22), and the growth time, eq. (24). The case with $`B_{\mathrm{max},\mathrm{en}}`$ depends both explicitly and implicitly on the angular velocity, the shear, and the Brunt-Väisälä frequency through eqs. (22), (23), (24), and (25). It is difficult to tell what the appropriate saturation field should be. Given the limitation of the current calculations, we can only argue that these fields are roughly representative of what one would expect during core collapse. We return to this issue in §4.

In the inner region (less than 7 km), for the $`B_{\mathrm{sat}}`$ case the field was amplified before bounce and then remained constant after bounce since the inner portions are sufficiently stable to convection to damp the MRI there. Near the boundary of the PNS, convective instability promotes the MRI giving the strong flat peak at $`10^6`$ cm in Fig. 8. The sharp decline on the outer boundary of the peak is due to damping by convective stability that diminishes in strength at larger radii, yielding a gradual radial increase in field strength at later time. For the $`B_{\mathrm{max},\mathrm{en}}`$ case, the instability criterion is fulfilled before bounce, but at the same time the condition in eq. (25) can be met: the epicyclic motion prevents growth of the MRI, and in our calculation the field growth is simply turned off while this condition holds. The result is smaller fields in the inner region for this case. This suppression of the fields for this case may be exaggerated by our prescription for $`N^2`$ in the absence of actual convective motion. In the outer region near the shock front, the angular velocity gradient is mostly negative and the shear is small. When ambient gas is convectively unstable the MRI amplifies the magnetic field for the $`B_{\mathrm{sat}}`$ case. For the $`B_{\mathrm{max},\mathrm{en}}`$ case, however, the condition in eq. (25) is met, again resulting in no amplification of the magnetic field. This interplay with convective instability results in the “spikiness” of the resulting field distribution. Once again, this suppression of the MRI may be artificial in these calculations. Some of the spikiness is exacerbated by evaluating gradients simply from zone to zone rather than attempting any smoothing. We have done calculations in which we have formally set $`N=0`$. The result is a substantially smooth distribution of fields at late times that resembles the upper envelope of the field distribution given in Fig. 8.

The $`\beta ^1`$ profiles have multiple peaks for the $`B_{\mathrm{sat}}`$ cases but has single peak for the $`B_{\mathrm{max},\mathrm{en}}`$ case. For the $`B_{\mathrm{sat}}`$ case, MHD activity is expected not only around the PNS boundary, but also the region between the boundary and the stalled shock. Note that at later times $`B_{\mathrm{sat}}`$ is larger than $`B_{\mathrm{max},\mathrm{en}}`$ for both the FH and MM profiles, and consequently $`\beta ^1`$ is larger for the $`B_{\mathrm{sat}}`$ case. The parameter $`\beta ^1`$ is especially sensitive to the amplitude of the magnetic field because it is proportional to the square of the magnetic field.

The time evolution of the field with saturation limit $`B_{\mathrm{sat}}`$ for the MM and FH profiles at 0.69 $`M_{\mathrm{}}`$ is plotted in Fig. 10. The initial solid body rotation case (§3.2) shows similar behavior. This mass shell contains the initial homologous core at bounce and later the hydrostatic core of the PNS. The radius of this mass shell is not constant. For this mass shell, the magnetic field is amplified shortly before bounce, and the amplitude quickly jumps to a few times $`10^{15}`$ G. The saturation field is reached in about 9 ms. In general, the quantity $`\mathrm{\Delta }B`$ in eq. (26) becomes smaller with time and becomes zero when the field reaches the saturation limit. In practice, if the shear drops so that the saturation field decreases, the field can decline as it approaches the saturation limit. There is no field dissipation in our treatment, but we do see the field peak then relax to a lower level in some portions of the PNS.

For the cases with initial differential rotation, the peak values of the magnetic field at the end of our calculation at 387 ms after bounce are $`B_{\mathrm{sat}}=3.5\times 10^{16}`$ G and $`B_{\mathrm{max},\mathrm{en}}=4.5\times 10^{15}`$ G for the FH profile, and $`B_{\mathrm{sat}}=3.3\times 10^{16}`$ G and $`B_{\mathrm{max},\mathrm{en}}=3.9\times 10^{15}`$ G for the MM profile. The amplitude of the magnetic field is remarkably high and above the QED limit ($`B_\mathrm{Q}=4.4\times 10^{13}`$ G), but remains less than equipartition. For the case of $`B_{\mathrm{sat}}`$, $`P_{\mathrm{mag}}`$ is $``$ 10% of $`P_{\mathrm{gas}}`$, and magnetic buoyancy may limit growth of the magnetic field (Wheeler, Meier, & Wilson, 2002).

### 3.2 Initial Solid Body Rotation Case

Even for an initial solid body rotating iron core, collapse generates a differential rotation profile similar to that of the FH and MM profiles (Figs. 7, 11). The value of $`\mathrm{\Omega }(r)`$ depends on $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$: for the solid body case, since we limit our study to $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}0.2`$ rad s-1, the angular velocity of the PNS remains only of order 100 rad s-1. For the solid body case, the ratio $`\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{kep}}`$ is expected to peak in the outer regions, since for a given $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$ the rotational velocity is much larger at larger radius for solid body rotation than for differential rotation. Since in this outer region all three profiles have similar specific angular momentum, the ratio $`\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{kep}}`$ is also similar among these three for our choices of $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$.

Despite the smaller initial angular velocity in the initial solid rotation case, the magnetic field is amplified to a large value for the $`B_{\mathrm{sat}}`$ case (Fig. 12). At 387 ms after bounce, the peak value of $`B_{\mathrm{sat}}`$ is $`9.0\times 10^{15}`$ G and $`2.3\times 10^{14}`$ G for $`B_{\mathrm{max},\mathrm{en}}`$.

The ratio $`\mathrm{\Omega }/\mathrm{\Omega }_{\mathrm{kep}}`$ peaks just inside of the stalled shock, but the resulting magnetic field is not maximized at this location. The maximum magnetic field is obtained at the boundary of the PNS where the shear is the strongest. This is a good illustration of how the shear, rather than angular momentum, powers the MRI . Since the magnitude of the shear is smaller than in the cases of the MM and FH profiles with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=1.0`$ rad s-1 (Fig. 4), the magnetic field strength is also smaller.

### 3.3 MHD Outflow Power

We expect that the magnetic field generated by the MRI will power MHD bi-polar outflow. Employing the characteristic power of a Blandford & Payne type MHD outflow (eq. (28)), the outflow luminosity $`L_{\mathrm{MHD}}`$ and the dimensionless parameter $`\nu `$ are calculated for the three initial rotational profiles, MM, FH, and solid body. The profiles of this MHD luminosity mimic those of the magnetic field. The parameter $`\nu `$ is always less than 0.11 for all cases, verifying our assumption that the resulting flows would be of the non-relativistic Blandford-Payne type. Fig. 13 gives the MHD luminosity and the parameter $`\nu `$ at different times for the MM and FH profiles using the saturation prescription $`B_{\mathrm{max},\mathrm{en}}`$. For the cases with initial differential profiles, the peaks of the MHD luminosity are at the boundary of the PNS which coincides with peaks in $`\nu `$ for the MM and FH profiles. The parameter $`\nu `$ peaks inside the stalled shock for the initial solid rotation case due to the relatively higher MHD luminosity and lower critical luminosity in that region.

Our calculations are limited to sub–Keplerian rotation and sub–equipartition fields, and yet they potentially produce significant MHD luminosity: for the saturation field $`B_{\mathrm{sat}}`$, the maximum values 387 ms after bounce are 6.6 $`\times 10^{53}`$ erg s-1 for FH, 4.2 $`\times 10^{53}`$ erg s-1 for MM, and 6.5 $`\times 10^{51}`$ erg s-1 for the initial solid rotation. For the saturation field $`B_{\mathrm{max},\mathrm{en}}`$, the peak values of MHD luminosity 387 ms after bounce are 1.4 $`\times 10^{52}`$ erg s-1 for FH, 9.0 $`\times 10^{52}`$ erg s-1 for MM, and 2.7 $`\times 10^{48}`$ erg s-1 for initial solid rotation. The investigation of how the MHD luminosity can be turned into a bi-polar flow is left for future work, although we outline some possibilities in the discussion below.

## 4 Discussion and Conclusions

No one would doubt that the progenitors of core collapse supernovae rotate and possess some magnetic field. Pulsars as remnants of core collapse are manifestly rotating and magnetized. The question has always been whether rotation and magnetic fields would be incidental perturbations or a critical factor in understanding the explosion. We have shown here that with plausible rotation from contemporary stellar evolution calculations and any finite seed field with a component parallel to the rotation axis, the magnetorotational instability can lead to the rapid exponential growth of the magnetic field to substantial values on times of a fraction of a second, comparable to the core collapse time. Even a relatively modest value of initial rotation gives a very rapidly rotating PNS and hence strong differential rotation with respect to the infalling matter.

This result seems reasonably robust. The reason is that the instability condition for the MRI is basically only that the gradient in angular velocity be negative. This condition is broadly satisfied in core collapse environments. Rotation can weaken supernova explosions without magnetic field (Fryer & Heger, 2000); on the other hand, rotational energy can be converted to magnetic energy that can power MHD bi-polar flow that may promote supernova explosions. The implication is that rotation and magnetic fields cannot be ignored in the core collapse context.

Even artificially limiting the post-collapse rotation to sub–Keplerian values as done here, we find fields in excess of $`10^{15}`$ G near the boundary of the neutron star. While this field strength is sub–equipartition, the implied MHD luminosities we have derived are of order $`10^{52}`$ erg s-1. This is a substantial luminosity and could, alone, power a supernova explosion if sustained for a sufficiently long time, a fraction of a second. As pointed out by Wheeler, Meier & Wilson (2002), the fields do not have to be comparable to equipartition to be important because they can catalyze the conversion of the large reservoir of rotational energy into buoyant, bi-polar MHD flow. Higher rates of initial rotation that are within the bounds of the evolutionary calculations could lead to even larger post-collapse rotation and even larger magnetic fields. If the initial rotation of the iron core proves to be substantially lower than we have explored here, then the MRI would be of little consequence to the explosion. The MHD luminosities derived here are comparable to the typical neutrino luminosities derived from core collapse, $`10^{52}`$ erg s-1. One important difference is that the matter beyond the PNS is increasingly transparent to this neutrino luminosity, whereas the MHD power is deposited locally in the plasma. Another difference is that the neutrino luminosity is basically radial so it resists the inward fall of the collapse, the very source of the neutrino luminosity itself. In contrast, hoop stresses associated with the magnetic field (see below) will tend to pull inward and force matter selectively up the rotation axis.

We note that for complete self-consistency, one should apply the MRI to the evolution of rotating stars where even a weak field renders the Høiland dynamical stability criterion “all but useless” in the words of Balbus & Hawley (1998). Some steps to explore magnetic viscosity in stellar evolution have been taken. Spruit & Phinney (1998) have invoked a magnetic viscosity that efficiently couples the inner core to the outer expanding giant envelope until the core contraction is more rapid than its rotation. They concluded that if the decoupling is as late as the phase of carbon depletion, the iron core might be rotating with an extremely low angular velocity, $`10^4`$ rad s-1, yielding an initial PNS rotation period of 100 s or an angular velocity of $``$ 0.06 rad s-1. This condition requires that the field be of order the rotation energy of the core, but if this condition holds then field amplification by the MRI in the collapse would be negligible. On the other hand, recent calculations by Heger & Woosley (2002) based on a prescription for magnetic viscosity by Spruit (2002) yield much more rapidly rotating iron cores. Heger & Woosley (2002) find PNS rotation rates of 4 to 8 ms, consistent with the values we have explored here. Clearly, much more must be done to understand the magnetorotational evolution of supernova progenitors.

As expected, the shear and hence the saturation fields are often highest at the boundary of the PNS where strong MHD activity is anticipated. Unexpectedly, with our fiducial saturation field, the field can also be large within the standing shock compared to the local pressure since the shock compression there naturally leads to shear in a rotating environment and since gas pressure is relatively low. The strength of this secondary peak is about the same for the initial solid body profile with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=0.2`$ rad s-1 and for the MM and FH profiles with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=1.0`$ rad s-1. We have had to limit $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$ to smaller values in the case of initially constant angular velocity in order to not violate the sub–Keplerian condition at larger radii after collapse. The maximum magnetic fields achieved are generally about the same strength within a factor of 10 for all three initial angular velocity prescriptions we have explored, although the value of $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}`$ was chosen to be substantially smaller in the initially solid rotation case compared to those where there is a gradient in $`\mathrm{\Omega }`$.

The configuration of the magnetic field in a precollapse iron core is not well understood. In this calculation we have assumed there exists a seed vertical field to calculate the growth of the field due to the MRI; however, the MRI can amplify other components of the magnetic field. The final configuration of the magnetic field after collapse may be less uncertain since the system has a strongly preferred direction due to rotation. In the context of accretion disks, Hawley et al. (1996) show that their fiducial run with initial random magnetic field configuration results in 9% radial, 88 % toroidal, and 3 % vertical components. Most of the shear is in the radial direction, so the radial component is greatly amplified by the MRI and turned into toroidal field due to differential rotation (Balbus & Hawley, 1998). The dominant component is most likely to be the toroidal field.

Another uncertainty is the rotational profile. Although we have formally considered the MRI in a geometry with arbitrary pressure gradients, in practice our assumption of conservation of angular momentum in spherical shells effectively restricts our analysis of the resulting magnetic fields to the equatorial plane. It is not clear what profile to use in the PNS, since, we note, even the rotational profile of the Sun is not well understood. A full understanding of the rotational state of a PNS remains a large challenge.

We have assumed various prescriptions for the saturation field. All are variations on the theme that, within factors of order 2$`\pi `$, the saturation field will be given by the condition $`v_\mathrm{A}r\mathrm{\Omega }`$. Fig. 14 compares the velocities of sound ($`c_\mathrm{s}`$), Keplerian rotation ($`v_{\mathrm{kep}}`$) and model rotation ($`r\mathrm{\Omega }`$), and the Alfvén ($`v_\mathrm{A}`$) velocity corresponding to $`B_{\mathrm{sat}}`$ for the solid, MM and FH profiles we adopted. For the sub–Keplerian rotational velocity profiles we employed, the magnetic field saturates at $`v_\mathrm{A}r\mathrm{\Omega }`$. For the case of $`B_{\mathrm{max},\mathrm{en}}`$, $`v_\mathrm{A}<r\mathrm{\Omega }`$ since $`B_{\mathrm{max},\mathrm{en}}<B_{\mathrm{sat}}`$ for the MM and FH profiles, and $`B_{\mathrm{max},\mathrm{en}}`$ is especially small for the initial solid rotation profile. In all cases, the saturation field is sub–equipartition with Alfvén velocity less than Keplerian. In the numerical disk simulations, about 20 rotations are required to reach saturation. The region of maximum shear in these calculations, around 15 km, typically has an angular velocity of 500 rad s-1 or a period of about 0.013 s. That means that by the end of the current calculations at 0.387 s, there have been about 30 rotations. Although the prescriptions for the growth and saturation fields we use here are heuristic, this aspect of our results is certainly commensurate with the numerical simulations of the MRI.

The issues of the saturation field and the nature of astrophysical dynamos are still vigorously explored. There are concerns that even weak fields may limit turbulent cascade at the smallest scales and hence suppress dynamo processes. Vainshtein, Parker & Rosner (1993) point out that this may be avoided by the formation of Coriolis-twisted loops that develop in shearing strong toroidal fields. Papaloizou and Szuszkiewicz (1992) presented a global stability analysis for conditions where the angular velocity is constant along magnetic lines of force that generally agreed with the local MRI stability analysis of Balbus & Hawley (1998) for weak fields. For strong fields, however, their analysis suggested that fields corresponding to the saturation condition $`v_\mathrm{A}\sqrt{R\mathrm{\Omega }c_\mathrm{s}}`$ would be unstable. For the sub–Keplerian conditions we consider here, this saturation criterion would yield a ratio of magnetic to total pressure that is larger by a factor of $`c_\mathrm{s}/r\mathrm{\Omega }`$ than the criteria we used (see Fig. 14). Vishniac & Cho (2001) conclude that the MRI has the required properties for a dynamo, anisotropic turbulence in a shearing flow, to generate both disordered and ordered fields of large strength. The saturation limits we have adopted here are consistent with those found in numerical calculations of the MRI saturation, but this topic clearly deserves more study.

We have made one major assumption that must be checked more carefully in future work. That is that we can apply the MRI linear stability analysis derived for no background radial flow to a dynamic situation where radial inflow and convective motions are rapid. The character of the MRI has never been investigated in this regime with finite background flow. One rationale for our assumption is that at a given radius one can transfer to a co-moving frame where the local radial velocity is zero. This is formally correct, but since we then make assumptions about modes with wavelength comparable to the size of the system, this formality may be inapplicable to the problem we want to solve, the nature of the saturation fields. An ex-post factor justification for our results is that after collapse when the PNS settles nearly into hydrostatic equilibrium, the radial motions are, in fact, small and the traditional MRI is probably valid as illustrated in Fig. 15. The action of the MRI in the vicinity of the standing shock and the fields we derive there may, however, be questionable on this basis. This topic deserves more study as a general point of physics, not merely for the application to core collapse. We also note that the MRI could be altered in the context of large radiation pressure. The radiation pressure is not a major component of the total pressure in the current calculations, so we do not believe this to be a factor here, but this is worth considering in principle.

We have also neglected several feed-back processes. Among these are the build up of magnetic pressure and hoop stresses, the viscous coupling of shells that will tend to suppress the differential rotation, the effects of neutrino viscosity on turbulence, the effects of the magnetic field on neutrino transport, and the effects of centrifugal forces. We will describe some of these issues briefly.

Both the magnetic pressure and the magnetic viscosity are small for the sub–Keplerian conditions explored here. For most cases $`\beta ^1`$ is less than 0.1 for the conditions we have assumed, (the $`B_{\mathrm{sat}}`$ case for the FH profile with $`\mathrm{\Omega }_{0\mathrm{\_}\mathrm{c}}=1.0s^1`$ pushes this limit), so the direct dynamical effect of the magnetic field is expected to be small. The viscous time scale is $`\tau _{\mathrm{vis}}(\alpha \mathrm{\Omega })^1(\mathrm{r}/\mathrm{h})^2`$, where $`\alpha `$ is the viscosity parameter and h is the vertical scale height, with h $``$ r for our case. For a magnetically-dominated viscosity,

$$\alpha \frac{B_\mathrm{r}B_\varphi }{4\pi P}=\left(\frac{B_\mathrm{r}}{B_\varphi }\right)\frac{B_\varphi ^2}{4\pi P}2\left(\frac{B_\mathrm{r}}{B_\varphi }\right)\beta ^1.$$
(31)
With this expression for $`\alpha `$, the viscous time becomes:

$$\tau _{\mathrm{vis}}\frac{1}{2}\left(\frac{B_\varphi }{B_\mathrm{r}}\right)\left(\frac{1}{\beta ^1\mathrm{\Omega }}\right)\mathrm{\Omega }^1,$$
(32)
where the final inequality follows from $`B_\varphi >B_\mathrm{r}`$ and $`\beta ^1<1`$. This prescription for viscosity is reasonable in the absence of convection. In the portions of the structure that are convective, the viscosity could be enhanced significantly.

In addition to dynamic effects, we note that the fields generated here are well above the QED limit ($`B_\mathrm{Q}=4.4\times 10^{13}`$ G). In this exotic regime, such a strong field has radiative and thermodynamic effects (Duncan, 2000), although it is not clear that these have profound effects on the dynamics in core collapse supernovae.

We have not discussed the role of neutrinos here, although the processes of neutrino loss and de-leptonization are included in our calculation of the cooling PNS. It is possible that the neutrino flux affects the magnetic buoyancy (Thompson & Murray 2002) and that the magnetic fields affect the neutrino emissivity (Thompson & Duncan, 1996) and interactions with the plasma (Laming 1999). The time scale for shear viscosity due to neutrino diffusion is much longer than the times of interest here, although magnetic fields and turbulence can make it shorter (Goussard et al., 1998). The MRI provides magnetic field and turbulence, so this issue deserves further study. In addition to affecting the shear, the neutrino viscosity might also affect the turbulence needed to make the MRI work. Using the expression for the neutrino viscosity at sub-nuclear densities from Goussard et al. (1998),

$$\eta _\nu 2\times 10^{23}\left(\frac{T}{10\mathrm{MeV}}\right)^2\rho _{13}^1\mathrm{g}\mathrm{cm}^1\mathrm{s}^1,$$
(33)
gives a Reynolds number

$$Re\frac{vR\rho }{\eta _\nu }\frac{\mathrm{\Omega }R^2\rho }{\eta _\nu }10^5$$
(34)
near the PNS boundary at $`R10^6`$ cm with $`\mathrm{\Omega }1000`$ rad s-1. This could decrease to Re $`100`$ near the standing shock where the density is lower. The former is probably large enough to sustain turbulence in the strong shear, but the latter might not be.

When viscosity dominates the dissipation in a collisional plasma, the growth condition that the field growth time be less than the dissipation time yields a constraint on the magnetic field (Balbus & Hawley, 1998):

$$B\left(\frac{15\pi }{8}\nu \rho \mathrm{\Omega }\right)^{1/2}=5\times 10^7\mathrm{G}\mathrm{\Omega }_3^{1/2}\left(\frac{\mathrm{T}}{10\mathrm{MeV}}\right)^{5/4}$$
(35)
This equation implies that B must exceed $`10^8`$ G near the boundary of the PNS where T $`10`$ MeV and smaller values in cooler portions at larger radii. Even compression of moderate fields in the iron core should exceed this threshold. We note that if one were to use the expression for the neutrino viscosity from eq. (33), the corresponding limit on the field would exceed $`10^{13}`$ G for similar parameters. This limit may not be relevant, however, since in order to use the fluid equations, rather than a kinetic theory, to describe the instability and its damping, the physical lengthscale associated with the viscosity, the mean free path of the dissipative particles, must be substantially less than the damping scale of the turbulence. This is not likely to be the case for neutrinos in the vicinity of the region of maximum shear in this problem. The damping length where Re $``$ 1 would be about $`10^5`$ of $`10^6`$ cm or about 10 cm from eq. (34), whereas the mean free path of the neutrinos (neglecting blocking and other complications) is about $`10^4`$ cm at the density of $`10^{13}`$ g cm-3 and temperature of about 10 MeV that characterize the region of maximum shear (Arnett, 1996). There may be regions deep in the neutron star where the neutrino mean free path is short enough that this becomes an issue, but we do not think that the neutrinos can affect the eddy turnover and instability of the MRI in regions that are significant for maximum field growth.

An obvious imperative is to now understand the behavior of the strong magnetic fields we believe are likely to be attendant to any core collapse situation. The fields will generate strong pressure tensor anisotropies that can lead to dynamic response even when the magnetic pressure is small compared to the isotropic ambient gas pressure. As argued in Wheeler, Meier & Wilson (2002), a dominant toroidal component is a natural condition to form a collimated magnetorotational wind, and hence polar flow. A first example of driving a polar flow with the MRI is given by Hawley & Balbus (2002). A key ingredient to force flow up the axis and to collimate it is the hoop stress from the resulting field. We have examined the acceleration implied by the hoop stresses of the fields we have derived here, $`a_{\mathrm{hoop}}=B_\varphi ^2/4\pi \rho r`$. We find that the hoop stresses corresponding to the peak saturation fields can be competitive with, and even exceed, the net acceleration of the pressure gradient and gravity. The large scale toroidal field is thus likely to affect the dynamics by accelerating matter inward along cylindrical radii. The flow, thus compressed, is likely to be channeled up the rotation axes to begin the bi-polar flow that will be further accelerated by hoop and torsional stresses from the field, the “spring and fling” outlined in Wheeler et al. (2002). The MRI is expected to yield a combination of large scale and small scale magnetic fields. Even the small scale fields in the turbulent magnetized medium may act like a viscoelastic fluid that would tend to drive circulation in along the equator and up the rotation axis (Williams, 2002) where the same small scale fields could collimate and stabilize the flow even in the absence of large scale toroidal fields (Li, 2002). On the other hand, small scale fields can result in dissipation by reconnection, an issue that we have not treated here, relying implicitly on the numerical simulations that show a growth in the field that is not eliminated by such effects.

The dynamics of these jets may depart substantially from pure hydrodynamical jets, since they will tend to preserve the flux in the Poynting flow and reconnection can accelerate the matter (see Spruit, Daigne & Drenkhahn 2001 and references therein). The poloidal component of the field can be another contributor to plasma waves (Wheeler et al. 2000).

For a complete understanding of the physics in a core collapse supernova explosion, a combination of neutrino–induced and jet–induced explosion may be required. Understanding the myriad implications of this statement will be a rich exploration.

We thank Steve Balbus, Roger Blandford, Stirling Colgate, Peter Höflich, Julian Krolik, Pawan Kumar, John Scalo, Jim Stone, Chris Thompson, Ethan Vishniac, Peter Williams and Ellen Zweibel for valuable discussions. We also thank the anonymous referee for several valuable comments. JCW is grateful to the hospitality of the Aspen Center for Physics where this work germinated. This work was supported in part by NASA Grant NAG59302 and NSF Grant AST-0098644.