On the Origin of Pulsar and Magnetar Magnetic Fields

Observations have revealed a dichotomy in the population of neutron stars (NSs): flaring, slowly rotating magnetars with strong inferred magnetic fields are distinct from faster-rotating radio pulsars with weaker fields. Explaining these two populations remains an open problem in the modeling of core-collapse supernovae and in subsequent NS evolution. While a large number of core-collapse simulations successfully match observed supernova properties and leave behind bound NS remnants (Couch & O'Connor 2014; Lentz et al. 2015; Burrows et al. 2019, 2020; Müller et al. 2019; Vartanyan et al. 2019; Vartanyan & Burrows 2020; Bollig et al. 2021), it is less clear how such modeling can account for the diversity of NS magnetic fields seen in nature.

To date only approximately 30 magnetars have been discovered (Olausen & Kaspi 2014), ³ though their transient nature implies that 10% or more of all young NSs could fall into this population (Kaspi & Beloborodov 2017). The cataloged spin periods range from 2 to 12 s, and their surface dipolar fields as calculated from the periods and period derivatives (see Igoshev et al. 2021, for a recent review) are 10¹³–10¹⁵ Gauss. Magnetars are young, with most having characteristic spin-down ages of less than 10⁴ yr. They are distributed closer to the galactic plane than radio pulsars (Kaspi & Beloborodov 2017), and eight are associated with known supernova remnants (Olausen & Kaspi 2014).

This population stands in contrast to that of the isolated radio pulsars, of which there are over 3000 cataloged members (Manchester et al. 2005). ⁴ The bulk of pulsar periods are between 0.1 and a few seconds, and young pulsars have periods of 0.5 s or less (Narayan 1987; Emmering & Chevalier 1989). Observations are well matched by initial periods being 300 ± 150 ms (Faucher-Giguère & Kaspi 2006). Pulsar surface dipole magnetic fields span 10¹¹–10¹³ Gauss, with the youngest objects falling in the 10¹²–10¹³ Gauss range (Kaspi & Kramer 2016).

The key to this dichotomy lies in the magnetic fields with which NSs are endowed during core collapse. While magnetars also have notably slower rotation rates, this is consistent with being born rotating as fast as (or possibly faster than) their less magnetized counterparts and subsequently rapidly losing angular momentum due to strong magnetic braking. Indeed, population modeling by Jawor & Tauris (2022) shows that magnetar initial periods must be less than 2 s.

One idea is that the field strengths seen in newly born NSs are simply the result of amplification due to flux freezing during collapse (Ferrario & Wickramasinghe 2006, 2008; Hu & Lou 2009), with the diversity of outcomes reflecting the range of main-sequence stellar fields. Spruit (2008) counters that the change in radius from that of the stellar core to that of the NS could not amplify reasonable fields to magnetar levels. On the other hand, Cantiello et al. (2016) note asteroseismic modeling indicates some massive stars may well be sufficiently magnetized to become magnetars through this simple mechanism. Still, flux freezing can be only a part of the explanation, since, as Spruit (2008) notes, a strong pre-existing field would couple the core to the envelope, reducing core angular momentum below what is seen in young pulsars.

Alternatively, field amplification can proceed more dynamically. Rotation is likely an important ingredient, with the interplay of rotation and magnetic fields being an important part of precollapse stellar evolution (Heger et al. 2005) as well as appearing in the dynamics of the protoneutron star (PNS) itself (Nagakura et al. 2020). The growth of field via compression, linear and exponential dynamos, and the magnetorotational instability (MRI; Velikhov 1959; Balbus & Hawley 1991) might produce jet-launching magnetars in the cores of long-duration gamma-ray bursts (Obergaulinger & Aloy 2020; Aloy & Obergaulinger 2021; Obergaulinger & Aloy 2021a, 2021b) and/or hypernovae (Burrows et al. 2007). These models generally rely on rapid rotation and, indeed, with such short rotational timescales even very small seed fields can be amplified to magnetar levels in less than 1 s (Mösta et al. 2015).

A related mechanism for amplification is that of a convective dynamo. Given the expectation of convective regions in PNSs shortly after bounce (Burrows 1987; Burrows & Lattimer 1988), Duncan & Thompson (1992) estimate that strong convective dynamo action during this time can greatly amplify magnetic fields, explaining magnetars in particular.

Within the convective dynamo framework, there happens to be a natural explanation for a dichotomy of outcomes. As shown by Christensen & Aubert (2006) and Olson & Christensen (2006) in the context of local planetary systems, there are simple scaling laws obeyed by these systems. In particular, the existence of a dipole field depends critically on the Rossby number relating the rotational and convective timescales. Dynamos that are sufficiently rapidly rotating yield a saturated field dominated by a dipole whose strength scales with the power provided by buoyancy, while slower rotation leads to a broader multipolarity.

Fortuitously for our purposes, these results are generally found to hold so long as the magnetic Reynolds number is sufficiently large, as it is in the PNS case, and they are otherwise largely independent of dimensionless numbers associated with diffusive effects, including the magnetic Reynolds, Ekman, Prandtl, and magnetic Prandtl numbers (Christensen & Aubert 2006). While momentum, thermal, and magnetic diffusivities can vary wildly across astrophysical dynamos, PNS Rossby numbers are similar to those of solar-system bodies. The relationship between convective flux and dipolar field strength has been successfully extended from planetary bodies to low-mass, rapidly rotating stars, including both T Tauri stars and old M dwarfs (Christensen et al. 2009). We seek in this paper to extend this connection further to encompass PNS magnetic-field generation.

Some recent work has focused on simulating field growth in the PNS context, with varying strategies for simplifying the problem. Franceschetti & Del Zanna (2020) see dynamo action with a prescribed, stationary velocity field and a polytropic PNS. Reboul-Salze et al. (2021) explore MRI effects with an incompressible (and, therefore, nonconvective) model. With an imposed heat flux and the anelastic equations, Raynaud et al. (2020) find strong field growth in rapidly rotating models. Masada et al. (2021) use a nuclear equation of state and drive convection by fixing the lepton fraction profile, whose gradient drives PNS convection in nature; they see large-scale eddies more akin to meridional flows than standard small-scale turbulence, which go to the center of the model and generate a dipolar field even with slow rotation. None of these models directly include neutrino transport, and the physics they include and conclusions they reach are varied.

We believe it is worthwhile to explore PNS dynamos from the heretofore neglected angle of Rossby number and convective flux scaling, applying the lessons learned from large suites of dynamo simulations to the PNS regime. In this light, we analyze a number of core-collapse simulations performed using the Fornax code (Skinner et al. 2019; Vartanyan et al. 2019). These include three-dimensional models with progenitors ranging in mass from 9 M_⊙ to 25 M_⊙ (Burrows et al. 2020), where one new 9 M_⊙ model (unpublished) features initial rotation. While these simulations do not include magnetic fields, their advanced treatment of the other relevant processes of core collapse, including neutrino transport, neutrino–matter coupling, and the equation of state, allow us to use them to confidently anchor our analysis to physically realistic and self-consistent conditions. All of our models attain vigorous convection in and around a PNS; Figure 1 shows a rendering of one snapshot with velocity tracer particles and a density isosurface delineating the PNS.

Zoom In Zoom Out Reset image size

In Section 2, we analyze our models in terms of convective scaling relations. Section 3 explores the convective region in more detail. We discuss the implications of this analysis in Section 4. Details about diffusivities and dimensionless numbers can be found in the Appendix.

The Rossby number Ro = v_conv/ΩR characterizes the rate of rotation relative to eddy turnover times, with faster rotation corresponding to lower Rossby numbers. According to Christensen & Aubert (2006) and Olson & Christensen (2006), convective dynamos with a modified Rossby number

$$\begin{eqnarray}&&{\mathrm{Ro}}_{{\ell }}=\displaystyle \frac{\bar{{\ell }}}{\pi }\cdot \displaystyle \frac{{v}_{\mathrm{conv}}}{{\rm{\Omega }}R}\lesssim 0.12\end{eqnarray} \tag{ 1 }$$

should produce strong dipolar fields. Here, we take v_conv to be the radial contribution to v_turb,

$$\begin{eqnarray}&&{v}_{\mathrm{conv}}=\langle {\left({v}^{r}-{v}_{\mathrm{ave}}^{r}\right)}^{2}{\rangle }_{{\rm{\Omega }}}^{1/2},\end{eqnarray} \tag{ 2 }$$

thus more unambiguously associating it with convection. The average angular scale $\bar{{\ell }}$ is defined in the Appendix. The evolution of radial profiles of Ro_ℓ is shown in Figure 2. The 9 M_⊙ rotating model is always below the critical threshold near 30 km, and for several hundred milliseconds is below the threshold exterior to this region.

Zoom In Zoom Out Reset image size

In the Ro_ℓ ≲ 0.12 regime, Christensen et al. (2009) provide a scaling relation for the strength of the saturated dipole field, expressed in our notation as

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\mathrm{dip}}^{2}}{8\pi }={{cf}}_{\mathrm{ohm}}{f}_{\mathrm{eff}}{(\langle \rho \rangle {F}_{\mathrm{out}}^{2})}^{1/3}.\end{eqnarray} \tag{ 3 }$$

Here c = 0.63 is a constant of proportionality fit by Christensen et al. (2009); we assume f_ohm = 1 (all dissipation is ultimately ohmic). F_out is the total flux (convective and neutrino) at the outer edge of the convection zone. Given the density and temperature scale heights H_ρ = −ρ/(d ρ/dr) and H_T = −T/(dT/dr), as well as the convective flux F_conv, the efficiency factor is

$$\begin{eqnarray}&&{f}_{\mathrm{eff}}=\left\langle {\left(\displaystyle \frac{{H}_{\rho }}{{H}_{T}}\right)}^{2/3}{\left(\displaystyle \frac{\rho }{\left\langle \rho \right\rangle }\right)}^{1/3}{\left(\displaystyle \frac{{F}_{\mathrm{conv}}}{{F}_{\mathrm{out}}}\right)}^{2/3}\right\rangle .\end{eqnarray} \tag{ 4 }$$

All averages are volume averages over the convective region.

In order to apply this scaling, we need to measure the fluxes in our simulations. Because they are not static, we are careful to do this in the reference frame in which the material at a given radius and time has a vanishing average radial momentum. The convective enthalpy flux comes from terms involving the kinetic energy, internal energy, and pressure:

$$\begin{eqnarray}&&{F}_{\mathrm{conv}}={\left\langle \left(\displaystyle \frac{1}{2}\rho {v}_{\mathrm{turb}}^{2}+u+p\right){v}_{\mathrm{conv}}^{r}\right\rangle }_{{\rm{\Omega }}}.\end{eqnarray} \tag{ 5 }$$

For comparison, we also consider the neutrino flux, which is simply the angle-averaged radial component in this frame:

$$\begin{eqnarray}&&{F}_{\nu }={\left\langle {F}_{{\nu }_{e}}^{r}+{F}_{{\bar{\nu }}_{e}}^{r}+{F}_{{\nu }_{x}}^{r}\right\rangle }_{{\rm{\Omega }}},\end{eqnarray} \tag{ 6 }$$

where the three neutrino species followed are the ν_e , ${\bar{\nu }}_{e},$ and the sum of the ν_μ , ${\bar{\nu }}_{\mu }$, ν_τ , and ${\bar{\nu }}_{\tau }$ neutrinos aggregated into one species. Multiplying by 4π R² yields the corresponding luminosities L_conv and L_nu, which we compare for the 9 M_⊙ models in Figure 3. At all times, the convective flux is strongly peaked inside the PNS convective region, as characterized further in Section 3. The modestly rotating 9 M_⊙ model has convection going all the way to the center. The neutrino luminosities grow to be much larger, and they are relatively constant in radius outside 30 km.

Zoom In Zoom Out Reset image size

In Figure 4, in order to provide a sense of how strongly convection scales with progenitor mass, we plot the convective luminosities for all other three-dimensional models from Burrows et al. (2020) and the unpublished rotating 9 M_⊙ model. The peak luminosity grows with progenitor mass, reaching values approximately four times higher in the 25 M_⊙ case relative to the 9 M_⊙ case.

Zoom In Zoom Out Reset image size

Several features of PNS convective flux are more easily illustrated in mass coordinates. Figure 5 shows a mass–time diagram of L_conv for three models. For the 9 M_⊙ nonrotating model, there is appreciable convection in the center of the PNS only after ∼1.7 s. For the other nonrotating models, it generally takes as much time or longer to achieve this state. In contrast, the moderately rotating 9 M_⊙ model has convective flux throughout most of the mass as early as ∼0.5 s. ⁵ Furthermore, we consistently find a radiative region in the outer ∼0.1 M_⊙ of the PNS. This is independent of progenitor mass, PNS mass, and rotation. The solid lines in Figure 5 denote the PNS boundary, defined arbitrarily as the 10¹¹ g cm⁻³ surface. This ∼0.1 M_⊙ layer is important in that any field generated in the convective region must be transported through it in order to contribute to the more readily observable NS surface and exterior fields.

Zoom In Zoom Out Reset image size

Defining the convective region to end where the convective luminosity drops below 1% of its peak value, we can apply Equations (3), (4), and (5) to our models to find the predicted saturated dipole field assuming the Christensen et al. (2009) scaling relations hold. This can be done for each snapshot, and the results as functions of time are shown in Figure 6. As there is snapshot-to-snapshot noise in the lines, for clarity we smooth them with a Gaussian filter (10 ms standard deviation).

Zoom In Zoom Out Reset image size

For all our models, the scaling predicts dipole field strengths of at least 5 × 10¹⁵ Gauss, with slightly larger values at higher progenitor masses. That is, there is enough convective flux to drive the creation of magnetar-level fields. Note that while assuming the scaling relations from Christensen et al. (2009) hold in this regime of Ro_ℓ ≲ 0.12, they only give the saturated field strength. That is, Figure 6 shows the time evolution of the saturated value, not the growth of the field itself.

Given how long the saturated strength remains large (beyond the time limits of any of our simulations), it seems there is a good chance for saturation to be reached. For illustrative numbers, we consider the rotating 9 M_⊙ model 500 ms after bounce. In this snapshot, there are 1030 sound-crossing times (center to ρ = 10¹¹ g cm⁻³ surface) per second; the fastest azimuthal velocities correspond to 136 rotations per second (and thus there are 136 Ro_ℓ convective turnovers per second). Comparing these numbers to the 14 e-foldings needed to grow the field by a factor of 10⁶, we infer that 1–2 s of sustained dynamo action could well be sufficient to generate magnetar fields from weak progenitor fields. We return to the issue of saturation while employing a linear analysis in Section 3.

It is important to note that we are not predicting such strong fields are ubiquitous. The scaling relation only predicts field strengths in the limit that the modified Rossby number is sufficiently low. For slower rotation, weaker, multipolar fields are predicted, with dynamos near the critical value possibly showing polarity reversals akin to those found in Earth's magnetic field. Our single rotating model has regions of sufficiently low Rossby number (Figure 2). Much of the convection, including the peak, is found only interior to 20 km (Figure 3, upper right panel).

As we are particularly interested in the broader implications of rotation as it applies to dynamos, we take a more detailed look at the rotation profile of the rotating 9 M_⊙ model. Figure 7 shows the radial profile of the angular velocity in the midplane at various times, plotted against mass coordinate. Note that the ordinate indicates the angular velocity at each abscissa, not the average interior to each abscissa. Angular velocity shows a strong peak near 1.3 M_⊙ (∼20–40 km) that grows with time. This is the effect of material that has accreted later, conserving its angular momentum and, thus, increasing its angular velocity, while not being able to efficiently penetrate the outermost envelope of the PNS to spin up the core at this early phase. Since the more coarse gridding in the deep interior can lead to large uncertainties near the coordinate singularity, we exclude the inner 0.1 M_⊙ of material from the plot.

Zoom In Zoom Out Reset image size

This rotation is not uniformly distributed in latitude. Figure 8 renders the three-dimensional azimuthal velocity field inside R = 100 km, 300 ms after bounce. The inner 50 km have been cut away and replaced with the 1000 and 2500 km s⁻¹ isosurfaces.

Zoom In Zoom Out Reset image size

In order for field amplification to occur, there must be a source of free energy in the system, which can then drive appropriate motion in the fluid. In the inner PNS context, there is indeed turbulence across a large range of scales, driven by convection. Here convection is not driven by an entropy gradient; indeed, entropy increases outward in the cores of most of our models. Instead, it is driven by the lepton composition gradient, which is maintained by the process of deleptonization as neutrinos escape from the outside inward (Keil et al. 1996; Dessart et al. 2006; Nagakura et al. 2020). We note that Olson & Christensen (2006) find that models with heat generated within the convective zone (as would happen with radioactive decay in a planetary mantle) do not display the same low-Rossby-number asymptotics as those with heat sourced only at boundaries. As we do not have explicit heat sources, nor are we considering the gain region outside the PNS core, we believe the latter models from Olson & Christensen (2006) are more applicable. Still, lepton-gradient-driven convection is different enough that further study is warranted.

Figure 9 highlights turbulent regions with the radial and tangential contributions to v_turb in the two 9 M_⊙ models, nonrotating and rotating. The nonrotating case shows the same qualitative features as found with the other progenitors: a well-defined convective region at the outer layers of the PNS (initially around 20–30 km), which slowly moves inward as the core cools and contracts, and another region of turbulent velocities exterior to this, going out beyond 150 km. The latter region is not convective at its base (interior to 60 km), as can be seen by the lack of turbulent radial velocities. This is consistent with the lack of convective flux here, as shown in Figure 3. Still, the material in this region, which has passed through the shock, maintains large velocities for extended periods of time, and thus may be relevant for magnetic-field generation.

Zoom In Zoom Out Reset image size

The rotating case displays several qualitative differences. The PNS convection zone can still be seen in radial velocity, though by around 500 ms after bounce its inner boundary extends all the way to the center. Furthermore, there is a distinct zone of turbulent motion in the inner 10 km of the PNS at early times. In general there are no parts of the region we are examining (interior to 150 km, and interior to the shock, which can be seen as the boundary between zero and nonzero velocities in the upper left of the lower right panel in Figure 9) that do not have significant turbulent velocities. We emphasize that bulk radial and azimuthal motion has been subtracted, leaving to show only smaller-scale motions.

A three-dimensional visualization of radial velocity in a single snapshot is provided in Figure 10. This shows the rotating 9 M_⊙ model 300 ms after bounce, out to R = 100 km. The large velocities from 50 to 100 km are visible, and one can clearly see the convective zone near R = 20 km. In this latter region, the convective cells' heights are a substantial fraction of the thickness of the turbulent layer itself.

Zoom In Zoom Out Reset image size

An important aspect of dynamos is the (kinetic) helicity in the small-scale velocity field v ', given by

$$\begin{eqnarray}&&h={\boldsymbol{v}}^{\prime} \cdot {\rm{\nabla }}\times {\boldsymbol{v}}^{\prime} .\end{eqnarray} \tag{ 7 }$$

Figure 11 shows the helicity on spherical surfaces (R = 20 km) 300 ms after bounce in the two 9 M_⊙ models, where the rotating model has azimuthally elongated features near the midplane. The three-dimensional rendering in Figure 12 shows the same quantity for R < 100 km, with a wedge at R = 20 km, for the rotating model. The nonconvective outer layers of the PNS, with essentially no helicity, separate the turbulent regions with large helicities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, field growth should be aided by larger magnitudes of helicity. This can be seen by writing the induction equation for the mean magnetic field as

$$\begin{eqnarray}&&\displaystyle \frac{\partial \left\langle {\boldsymbol{B}}\right\rangle }{\partial t}={\rm{\nabla }}\times (\left\langle {\boldsymbol{v}}\right\rangle \times \left\langle {\boldsymbol{B}}\right\rangle +\alpha \left\langle {\boldsymbol{B}}\right\rangle -(\eta +\beta ){\rm{\nabla }}\times \left\langle {\boldsymbol{B}}\right\rangle )\end{eqnarray} \tag{ 8 }$$

(Charbonneau 2014), where η is the magnetic diffusivity, $\alpha =-{\tau }_{{\rm{c}}}\left\langle h\right\rangle /3$ drives mean-field amplification on small scales, $\beta ={\tau }_{{\rm{c}}}\left\langle {({\boldsymbol{v}}^{\prime} )}^{2}\right\rangle /3$ incorporates the fact that small-scale turbulent reconnection acts as an effective diffusivity, and τ_c is a correlation timescale for the turbulence. ⁶ As shown in Figure 13, regions of large helicity broadly trace those of large tangential turbulent velocity (see Figure 9).

Zoom In Zoom Out Reset image size

Rotating models already have extra kinetic energy from their bulk motions, compared to equal-progenitor-mass nonrotating models. Further, they can supply free energy sourced from their differential rotation. The above discussion shows that, even beyond these points, rotating models experience more turbulent and helical motion, which makes them more favorable for dynamo action.

The amenability for dynamo growth can be seen more quantitatively in a linear-mode analysis. As analyzed by Shibata et al. (2021) in an NS context (see Brandenburg & Subramanian 2005 for a general treatment), a mode with frequency ω (positive imaginary parts unstable), wavenumber k, and wavenumber k_∥ parallel to the rotation axis has the dispersion relation

$$\begin{eqnarray}&&i\,\omega =\eta {k}^{2}\pm \sqrt{{\alpha }^{2}{k}^{2}-i\alpha {k}_{\parallel }{S}_{{\rm{\Omega }}}},\end{eqnarray} \tag{ 9 }$$

where S_Ω = r(dΩ/dr) measures the shear in cylindrical radius r. The first term under the radical corresponds to the α² dynamo, while the second term gives rise to the α–Ω dynamo. Shibata et al. (2021) note that when a large-scale dynamo dominates (as we would expect to be associated with dipolar-dominated fields), the α² term becomes negligible compared to the α–Ω term. An unstable mode (with k_∥ = k) will then exist for

$$\begin{eqnarray}&&k\lt {k}_{\mathrm{crit}}={\left(\displaystyle \frac{| \alpha {S}_{{\rm{\Omega }}}| }{2{\eta }^{2}}\right)}^{1/3}.\end{eqnarray} \tag{ 10 }$$

We can estimate this critical wavenumber in our rotating 9 M_⊙ model, using radius r ∼ 20 km, eddy size L ∼ 5 km (Figure 10), eddy velocity v ∼ 10³ km s⁻¹ (Figure 9), helicity h ∼ 10¹¹ cms⁻² (Figure 13), and shear dΩ/dr ∼ 10⁻³ s⁻¹ cm⁻¹. Then α = τ h/3, with τ = L/v. We approximate magnetic diffusivity to be η ∼ β = τ v²/3 ∼ 10¹³ cm² s⁻¹, if we assume turbulent diffusion should be substituted for that due to microscopic conductivity (Figure 15 in the Appendix). The result is instability for wavelengths greater than ∼7 km.

As shown in Shibata et al. (2021), the fastest growing mode has a wavenumber k_fast that is ∼0.4 times the critical wavenumber, corresponding to a wavelength of ∼20 km in our case. The corresponding growth rate is

$$\begin{eqnarray}&&{\mathfrak{I}}\left({\omega }_{\mathrm{fast}}\right)=\displaystyle \frac{3}{8}{\left(\displaystyle \frac{{\alpha }^{2}{S}_{{\rm{\Omega }}}^{2}}{2\eta }\right)}^{1/3}\sim 600\ {{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 11 }$$

Even considering the approximate nature of these estimates and the simple, linear nature of this analysis, ⁷ there are plenty of e-foldings to reach saturation in 1–2 s. In the highly conductive, convective PNS context, even modest rotation will lead to dynamo growth of field; the only question is what the nonlinear saturation strength is and whether rotation is strong enough to produce a dipole-dominated field, both of which we have discussed in Section 2.

We can make a simple estimate of the small-scale magnetic-field strengths we expect based on energy considerations. Using only the turbulent velocity field, we can ask what equivalent field has 10% of the fluid's kinetic energy on a given spherical shell at a given time. Radial profiles of this estimate, calculated for all models 450 ms after bounce, are shown in the upper panel of Figure 14. The three models that fail to explode are marked with dotted lines; the 9 M_⊙ rotating model is indicated with a dashed line. All models have a peak equivalent field strength within a factor of a few of 10¹⁵ Gauss corresponding to the PNS convection zone at around 15–20 km, with the rotating model maintaining this strength all the way to the center. In the range 30–150 km, both of the 9 M_⊙ models have much less turbulent kinetic energy than the other successfully exploding models, while those others show little variation with progenitor mass.

Zoom In Zoom Out Reset image size

Taking the peak equivalent field strength from each snapshot, we construct time series of maximum magnetic-field strengths, shown in the lower panel of Figure 14. All our models sustain a region where 10% of the turbulent kinetic energy corresponds to at least 5 × 10¹⁴ Gauss throughout the simulations, which go beyond 1 s after bounce in some cases. Note we are not suggesting that dipole fields would necessarily arise with such magnitudes; estimating the dipole strength requires the Rossby number and convective flux analysis of Section 2. Equipartition merely shows that there is sufficient free energy in the system to power magnetar fields, even if much of it only contributes to toroidal fields or to structures that are incoherent on large scales.

The observed NS population displays a wide range of magnetic-field strengths and spins. While to some degree evolution after the birth of the NS can lead to different outcomes (most notably the millisecond pulsars that have undergone accretion from a companion, as described in Bhattacharya & van den Heuvel 1991 or Phinney & Kulkarni 1994), much of the diversity must be attributed to the range of conditions before and during core collapse. In particular, the dichotomy between magnetars and the standard radio pulsars demands explanation.

There are a number of mechanisms that can amplify the magnetic fields in the cores of massive stars during collapse, resulting in pulsar or possibly magnetar strengths. Flux freezing, the MRI, and dynamos have all been considered, and it is likely that each plays some role in nature. Flux freezing has difficulty attaining magnetar fields, however, and both it and the MRI do not naturally result in a dichotomy of outcomes given a monomodal distribution of initial conditions.

On the other hand, convective dynamos can produce fields with a natural dichotomy of outcomes. As first realized in the context of planetary bodies, systems with modified Rossby numbers Ro_ℓ below a threshold of ∼0.12 acquire predominantly dipolar fields, while systems with relatively slower rotation (or faster convective motions) become multipolar, with relatively weaker dipoles (Christensen & Aubert 2006; Olson & Christensen 2006). The behavior near the transition is less clear; ⁸ this does not qualitatively affect our arguments, though should be considered in future work that quantifies the distribution of PNS Rossby numbers.

In the low-Rossby-number regime, the resulting dipole strength has a simple scaling with convective flux that is independent of the classic dimensionless numbers, and this has been extended to and tested against T Tauri stars and old M dwarfs (Christensen et al. 2009). The Rossby-number dichotomy motivates us to analyze core-collapse models in the same framework, applying the Christensen et al. (2009) scalings to a new physical regime, that of a PNS at the center of a core-collapse supernova.

It has already been shown that vigorous PNS convection is generic in core collapse (Dessart et al. 2006; Nagakura et al. 2020), and so there will always be interaction between this motion and the magnetic fields present in the progenitor. Moreover, even a moderately rotating progenitor, such as our 9 M_⊙ model, can attain Ro_ℓ < 0.12 (Figure 2) in and around the PNS. That is, the extreme rotation at an appreciable fraction of breakup sometimes employed to maximize magnetic growth may be unnecessary for producing large fields.

We note Masada et al. (2021) find that even very slow rotation (periods longer than 160 ms) can produce large dipolar fields rapidly (∼6–10 ms per e-folding) given large eddies (and thus long turnover times) in deep-core convection. All of their simulations show larger convective cells than we see, and it is unclear how generic this behavior is, or whether it might be the result of, e.g., the way in which they impose Y_e profiles to capture neutrino effects. Their results illustrate how important the details of PNS convection can be for generating NS magnetic fields.

A somewhat different perspective is given by Raynaud et al. (2020), whose simulations differ from those of Masada et al. (2021) by (1) being anelastic; (2) using heat flux rather than composition to drive convection; (3) using generally larger momentum, thermal, and magnetic diffusivities; and (4) employing very rapid rotation (with periods approaching breakup at 2 ms). Raynaud et al. (2020) find magnetic energy densities exceeding kinetic energy densities, with saturated field strengths scaling with rotation rate, from which they infer magnetostrophic rather than turbulent force balance dominates (see Brun et al. 2015). Their fields grow at a rate of ∼70 ms per e-folding. Rapid rotation might be sufficient to generate magnetar fields without the aid of convective flux, but it need not be necessary. Indeed, requiring all magnetars to be rotating near breakup at birth would lead to observable effects in supernovae at very early (Burrows et al. 2007) or slightly later (Kasen & Bildsten 2010) times, which would be in tension with the low observed rates of particularly luminous transients.

Our own suite of radiation-hydrodynamic simulations, though lacking magnetic fields, employs full neutrino transport and a nuclear equation of state, and therefore should develop particularly realistic convection. With progenitors ranging from 9 M_⊙ to 25 M_⊙, they show sufficient convective flux (Figure 4) to support even magnetar fields according to the Christensen et al. (2009) scalings (Figure 6). Predicted field strength increases slightly with progenitor mass, consistent with more massive progenitors leaving behind a more massive PNS, which subsequently has more convective flux (Nagakura et al.2020).

That there is enough kinetic energy in PNS convection to correspond to strong magnetic fields in equipartition (Figure 14) is not new (Endeve et al. 2010, 2012; Obergaulinger et al. 2014). This fact suggests that on small scales magnetar-magnitude fields are naturally generated in every NS at birth. However, the dynamo studies of Christensen & Aubert (2006), Olson & Christensen (2006), and Christensen et al. (2009) provide reason to believe that when sufficient rotation is present, large-scale, dipolar fields of magnetar magnitude will in fact arise. This connection is strengthened by the presence of large values of helicity in the regions of interest (Figures 12 and 13), which drives large-scale field growth in mean-field theory (Charbonneau 2014).

The scaling results we employ only predict the state of the magnetic field once it has reached nonlinear saturation. In order for convective dynamos to be relevant to NS field amplification, there must be sufficient time for field growth while convection lasts (i.e., before the core has fully deleptonized and cooled). Rough estimates of relevant timescales (see Section 2), as well as more detailed linear-mode analysis (Section 3), indicate that there is ample time for exponential growth of at least a factor of 10⁶ in field strength. The 9 M_⊙ models we analyze are run to ∼2 s after bounce, and, while convection may begin to diminish by this time (Figure 3), conditions are still amenable to the generation of strong dipolar fields in the low-Rossby-number regime.

Our focus has been on using dynamo scaling relations to inform PNS field generation. We note, however, that data for PNS fields can inform the general study of dynamos. The objects we consider have very different diffusion physics compared to planetary interiors and low-mass stars, with particularly high conductivities and with neutrino-mediated momentum and thermal diffusivities. While some dimensionless characterizations of our models (e.g., Rossby, Reynolds, and Prandtl numbers) are roughly commensurate with studies in different regimes, others (e.g., magnetic Reynolds and magnetic Prandtl numbers, ${\mathrm{Re}}_{{\rm{m}}}$ and Pr_m) are much larger. We refer the reader to the Appendix for more details about the diffusivities and dimensionless numbers in our models.

These differences in dimensionless numbers may be relevant according to Brandenburg (2014), who finds that the ratio of kinetic to magnetic dissipation in simulations is roughly ${\Pr }_{{\rm{m}}}^{1/3}$. Extrapolating this result from the largest value in that study, Pr_m = 20, to Pr_m ∼ 10¹⁴ applicable to PNS convection, we might expect f_ohm ∼ 10⁻⁵ in Equation (3), meaning the expected dipole might be reduced by a factor of ∼10² from what is plotted in Figure 6. We emphasize, however, that this extrapolation is uncertain, and in fact Brandenburg (2014) mentions that it may be that at large ${\mathrm{Re}}_{{\rm{m}}}$ (which is much larger in a PNS than in any simulation), the dissipation ratio may lose its dependence on Pr_m. Moreover, as noted by Thompson & Duncan (1993), the nominal neutrino viscous damping scale is comparable to the neutrino mean free path, so any turbulent cascade can continue down to much smaller scales set by electron viscosity, where the magnetic Prandtl number is much closer to unity.

Though explicit numerical tests with realistic magnetic diffusivities are intractable, we assume, and there is reason to believe (Christensen et al. 2009), that for sufficiently large dimensionless numbers the behavior of dynamos should enter an asymptotic regime and lose its dependence on these numbers. This assumption, and indeed the range of validity for dynamo scaling, will be tested over time as more PNS simulations are compared to NS observations.

In conclusion, we find it plausible that PNS convection plays a large role in determining salient properties of NSs at birth. Convective motions can further augment fields that have already been amplified via flux freezing. If sufficient rotation exists—and we emphasize that "sufficient" need not be close to breakup—dynamo action will lead to a dipolar field with a strength scaling with convective flux. This channel can explain the magnetars, as the strong dipole field will naturally remove angular momentum on short timescales, leaving a slowly rotating, highly magnetized NS to be observed. Without sufficient rotation, the convective regions still exist, but no magnetar-like dipole will be created, and the resulting NS will take longer to spin down, explaining standard radio pulsars.

There is important work to be done to be sure this model holds together. While we can be sure that strong convection is common in PNSs, it is still unclear what distribution of rotational profiles and seed magnetic fields we should expect in the late stages of massive-star evolution. In addition, future PNS modeling would benefit from a more comprehensive treatment of the physics. The models we study here excel at hydrodynamics and neutrino transport, including a sophisticated equation of state and detailed neutrino–matter interactions, but they do not include the (generally energetically subdominant) magnetic fields we know must be present. The most recent PNS/magnetic-field studies by various groups (e.g., Franceschetti & Del Zanna 2020; Raynaud et al. 2020; Masada et al. 2021; Reboul-Salze et al. 2021) have begun to address magnetic-field generation through dynamos or the MRI, though at the expense of approximating the physics that drives PNS convection. What is needed in the next theoretical phase is the merger of state-of-the-art three-dimensional supernova codes incorporating detailed neutrino transfer, equations of state, and initial models with a sophisticated treatment of magnetic fields. The resulting radiation-magnetohydrodynamic code will enable evolutionary simulations of the three-dimensional PNS cores birthed by supernovae and carried to late times. The ultimate goal is to determine the mapping between the initial massive-star progenitors and the final masses, spins, magnetic dipoles, and magnetic multipolarity structures of the residual NSs. Such a theory should usefully connect the radio pulsar and magnetar communities with the supernova theory community studying the creation of the NSs upon which the former focus. This is a long-awaited goal that we suggest may soon be within reach.

We thank Eric Blackman, Romain Teyssier, Tianshu Wang, and Hiroki Nagakura for fruitful discussions and Joe Insley of Argonne National Laboratory for the graphical rendering of Figure 1. We acknowledge support from the U. S. Department of Energy Office of Science and the Office of Advanced Scientific Computing Research via the Scientific Discovery through Advanced Computing (SciDAC4) program and grant No. DE-SC0018297 (subaward 00009650) and support from the U. S. National Science Foundation (NSF) under grant Nos. AST-1714267 and PHY-1804048 (the latter via the Max-Planck/Princeton Center (MPPC) for Plasma Physics). The bulk of the computations presented in this work were performed on Blue Waters under the sustained-petascale computing project, which was supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters was a joint effort of the University of Illinois at Urbana–Champaign and its National Center for Supercomputing Applications. The two 9 M_⊙ models were simulated on the Frontera cluster (under awards AST20020 and AST21003), and this research is part of the Frontera computing project at the Texas Advanced Computing Center (Stanzione et al. 2020). Frontera is made possible by NSF award OAC-1818253. Additionally, a generous award of computer time was provided by the INCITE program, enabling this research to use resources of the Argonne Leadership Computing Facility, a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Finally, the authors acknowledge computational resources provided by the high-performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Princeton University Office of Information Technology, and our continuing allocation at the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U. S. Department of Energy under contract DE-AC03-76SF00098.

Software: Fornax (Skinner et al. 2019; Vartanyan et al. 2019), VisIt.

In a turbulent, convective, conducting fluid, we generally expect the character of the flow to depend on three diffusivities: those of momentum, ν (i.e., the kinematic viscosity); thermal energy, κ; and magnetic field, η. Viscosity determines how far down in scale turbulence will be driven, and, in combination with thermal diffusivity, determines in many contexts whether convection drives turbulence. Some magnetic diffusivity is required for successful dynamos, allowing reconnection to change the field topology. Too much, however, prevents magnetic fields from being sustained.

Burrows & Lattimer (1988) note that ν (see van Weert & van den Horn 1983) and κ (see van den Horn & van Weert 1984) are dominated in the optically thick regions by the effects of neutrinos diffusing through the fluid. Thompson & Duncan (1993) approximate the scalings

$$\begin{eqnarray}&&{\nu }_{\mathrm{diff}}=2\times {10}^{9}\ {\rho }_{14}^{-4/3}{T}_{30}f{({Y}_{{\rm{e}}})}^{-1}\ {\mathrm{cm}}^{2}\ {{\rm{s}}}^{-1}\end{eqnarray} \tag{ A1 }$$

and

$$\begin{eqnarray}&&{\kappa }_{\mathrm{diff}}=1.5\times {10}^{10}\ {\rho }_{14}^{-2/3}{T}_{30}^{-1}\ {\mathrm{cm}}^{2}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ A2 }$$

where ρ₁₄ = ρ/(10¹⁴ g cm⁻³), T₃₀ = k_B T/(30 MeV), and $f{({Y}_{{\rm{e}}})=({Y}_{{\rm{e}}}^{1/3}+{(1-{Y}_{{\rm{e}}})}^{1/3})}^{-1}$. Outside the neutrinospheres, however, the neutrino mean free paths begin to grow larger than the length scales of interest, and so the effective diffusivities must be modified. As is done in Guilet et al. (2015) in the context of MRI growth, for the purposes of our estimations we define

$$\begin{eqnarray}&&\nu =\displaystyle \frac{{\nu }_{\mathrm{diff}}}{1+{(\lambda /{H}_{p})}^{2}}\end{eqnarray} \tag{ A3 }$$

and

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{{\kappa }_{\mathrm{diff}}}{1+{(\lambda /{H}_{p})}^{2}},\end{eqnarray} \tag{ A4 }$$

where

$$\begin{eqnarray}&&\lambda =200\ {\rho }_{14}^{-1/3}{T}_{30}^{-3}f({Y}_{{\rm{e}}})\ \mathrm{cm}\end{eqnarray} \tag{ A5 }$$

(Thompson & Duncan 1993) and H_p = − p/(dp/dr) is the pressure scale height at a given radius.

Angle-averaged radial profiles of these diffusivities, calculated from our 9 M_⊙ model at various times, are shown in Figure 15. In the core of the PNS (approximately the inner 30 km), ν ranges from 10⁸ to 10¹² cm² s⁻¹, while κ ranges from 10¹⁰ to 10¹² cm² s⁻¹. They both range from 10¹¹ to 10¹³ cm² s⁻¹ in the region exterior to this (inside 150 km).

Zoom In Zoom Out Reset image size

Magnetic diffusivity arises from electrical resistivity, which is taken to be dominated by electron–proton scattering. Following Thompson & Duncan (1993), we define

$$\begin{eqnarray}&&\eta =\displaystyle \frac{1}{3\times {10}^{4}}\ {\rho }_{14}^{-1/3}{(5{Y}_{{\rm{e}}})}^{-1/3}\ {\mathrm{cm}}^{2}\ {{\rm{s}}}^{-1}.\end{eqnarray} \tag{ A6 }$$

Radial profiles of η are also shown in Figure 15, where it varies from 10⁻⁵ to 10⁻³ cm² s⁻¹ inside the PNS, increasing up to 10⁻² in the exterior region. If we let v_turb be the rms velocity of turbulence, subtracting off bulk radial motion at each radius, and, in the case of rotating models, net rotation at each radius and latitude, and if we define ${v}_{\mathrm{ave}}^{r}={\left\langle \rho v\right\rangle }_{{\rm{\Omega }}}/{\left\langle \rho \right\rangle }_{{\rm{\Omega }}}$, ${v}_{\mathrm{ave}}^{\theta }=0$, and ${v}_{\mathrm{ave}}^{\phi }=0$ (nonrotating models) or ${v}_{\mathrm{ave}}^{\phi }={\left\langle \rho v\right\rangle }_{\phi }/{\left\langle \rho \right\rangle }_{\phi }$, then

$$\begin{eqnarray}&&{v}_{\mathrm{turb}}={\left\langle {({\boldsymbol{v}}^{\prime} )}^{2}\right\rangle }_{{\rm{\Omega }}}^{1/2},\end{eqnarray} \tag{ A7 }$$

where ${\boldsymbol{v}}^{\prime} ={\boldsymbol{v}}-{{\boldsymbol{v}}}_{\mathrm{ave}}$. This provides a characteristic velocity scale.

Next, we define a characteristic turbulent length scale as follows. Decompose each component of the turbulent velocity field via real spherical harmonics Y_{ℓ m} , calculating the power in each harmonic degree ℓ up to a practical maximum of ${{\ell }}_{\max }=32$, and finding the average of ℓ weighted by this power:

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{\ell }m}=\oint ({\boldsymbol{v}}-{{\boldsymbol{v}}}_{\mathrm{ave}}){Y}_{{\ell }m}\,d{\rm{\Omega }},\end{eqnarray} \tag{ A8 }$$

with

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{\ell }}=\sum _{m=-{\ell }}^{{\ell }}{{\boldsymbol{v}}}_{{\ell }m}{Y}_{{\ell }m}\end{eqnarray} \tag{ A9 }$$

and

$$\begin{eqnarray}&&\bar{{\ell }}=\displaystyle \frac{{\sum }_{{\ell }=1}^{{{\ell }}_{\max }}{\ell }{\left\langle {v}_{{\ell }}^{2}\right\rangle }_{{\rm{\Omega }}}}{{\sum }_{{\ell }=1}^{{{\ell }}_{\max }}{\left\langle {v}_{{\ell }}^{2}\right\rangle }_{{\rm{\Omega }}}}.\end{eqnarray} \tag{ A10 }$$

Take the length scale to be

$$\begin{eqnarray}&&{L}_{\mathrm{turb}}=\displaystyle \frac{\pi R}{\bar{{\ell }}},\end{eqnarray} \tag{ A11 }$$

where R is the spherical radius.

Dimensionless numbers constructed from diffusivities provide important characterizations of the fluid dynamics. The Reynolds number

$$\begin{eqnarray}&&\mathrm{Re}\,=\,\displaystyle \frac{{v}_{\mathrm{turb}}{L}_{\mathrm{turb}}}{\nu },\end{eqnarray} \tag{ A12 }$$

averaged over angle in a given snapshot, is plotted in the first panel of Figure 16 for the suite of models. It lies between 10³ and 10⁵ in all models within the PNS core, delineated by a sharp change in density, and is perhaps an order of magnitude lower than this immediately exterior to the core. Curiously, $\mathrm{Re}$ is not so large as to prohibit resolution of the viscous scale in future simulations. The three models that do not explode (those with 13, 14, and 15 M_⊙ progenitors; Burrows et al. 2020), have clearly distinct $\mathrm{Re}$ profiles (and, indeed, profiles of other dimensionless numbers). Otherwise, there is little variation of $\mathrm{Re}$ with progenitor mass.

Zoom In Zoom Out Reset image size

The magnetic Reynolds number can be defined similarly:

$$\begin{eqnarray}&&{\mathrm{Re}}_{{\rm{m}}}\,=\,\displaystyle \frac{{v}_{\mathrm{turb}}{L}_{\mathrm{turb}}}{\eta }.\end{eqnarray} \tag{ A13 }$$

As shown in Figure 16, ${\mathrm{Re}}_{{\rm{m}}}$ is always above 10¹⁶, indicating a highly conducting plasma where magnetic diffusion at characteristic turbulent scales is quite small. Excluding models that do not explode, only the 9 M_⊙ models are outliers in ${\mathrm{Re}}_{{\rm{m}}}$, with the rotating model having larger values interior to 30 km and the nonrotating model having lower values exterior to this. We do not, however, expect the differences to result in qualitative differences in terms of turbulence and dynamo generation, given how large ${\mathrm{Re}}_{{\rm{m}}}$ is in all cases.

We can also construct the Prandtl number

$$\begin{eqnarray}&&\Pr =\displaystyle \frac{\nu }{\kappa }\end{eqnarray} \tag{ A14 }$$

and magnetic Prandtl number

$$\begin{eqnarray}&&{\Pr }_{{\rm{m}}}=\displaystyle \frac{\nu }{\eta },\end{eqnarray} \tag{ A15 }$$

shown in the bottom panels of Figure 16. While the former is close to unity over the region of interest, the large conductivity of the fluid results in values of Pr_m ranging from 10¹³ to 10¹⁶.