Three-dimensional Simulations of Magnetar-powered Superluminous Supernovae

We use the Larmor formula for the luminosity of the magnetar (Lyne & Graham-Smith 1990; Chen et al. 2016):

$$\begin{eqnarray}\begin{array}{rcl}{L}_{{\rm{m}}} & = & \displaystyle \frac{32{\pi }^{4}}{3{c}^{2}}{\left({\mathrm{BR}}_{\mathrm{ns}}^{3}\sin \alpha \right)}^{2}{P}^{-4}\\ & & \sim \,1.0\times {10}^{49}{B}_{15}^{2}{P}_{\mathrm{ms}}^{-4}\quad \mathrm{erg}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 1 }$$

where B₁₅ is the surface dipole field strength at the equator of the NS in units of 10¹⁵ G, P_ms is the initial spin period in milliseconds, R_ns is the radius of the NS and α is the inclination angle between the spin and magnetic axes. In our models, R_ns = 10⁶ cm and α = 30°. The spin energy is

$$\begin{eqnarray}&&E=\displaystyle \frac{1}{2}I{\omega }^{2}\sim 2\times {10}^{52}{P}_{\mathrm{ms}}^{-2}\quad \mathrm{erg}.\end{eqnarray} \tag{ 2 }$$

If the magnetic field is constant then the evolution of the luminosity, period, and spin energy of the magnetar over time can be approximated by

$$\begin{eqnarray}\begin{array}{rcl}{L}_{{\rm{M}}}(t) & \sim & {\left(1+t/{t}_{m}\right)}^{-2}{E}_{0}{t}_{m}^{-1}\ \mathrm{erg}\ {{\rm{s}}}^{-1},\\ P(t) & \sim & {\left(1+t/{t}_{m}\right)}^{1/2}{P}_{0}\ \mathrm{ms},\\ E(t) & \sim & {\left(1+t/{t}_{m}\right)}^{-1}{E}_{0}\ \mathrm{erg},\end{array}\end{eqnarray} \tag{ 3 }$$

where P₀ = P_ms(0), E₀ = E(P₀), and ${t}_{{\rm{m}}}\sim 2\times {10}^{3}{P}_{\mathrm{ms}}^{2}{B}_{15}^{-2}$ is the magnetar spin-down timescale.

The evolution of the magnetar is shown in Figure 1 for several magnetic fields and initial spins. The spin energies of the 1 and 10 ms magnetars are ∼2 × 10⁵² and 2 × 10⁵⁰ erg, respectively. Their initial luminosities can vary from ∼10⁴² to 10⁴⁹ erg s⁻¹ depending on period and magnetic field. The 1 ms magnetar with a magnetic field of 10¹⁴–10¹⁵ G can release 10⁵¹ erg in 10²–10⁴ s, with the remainder being emitted in just 10³–10⁷ s. Consequently, dipole radiation from the magnetar can efficiently convert its spin energy to luminosity.

Zoom In Zoom Out Reset image size

Which choice of magnetar parameters is most likely to produce a superluminous event? We show the magnetar luminosity over time in days in Figure 1(c). The average luminosity of the SN itself at early times is plotted as the horizontal dashed line at 10⁴³ erg s⁻¹. The magnetar produces a superluminous event only if its luminosity exceeds that of the SN by the time photons from the NS can diffuse out of the ejecta, typically on timescales of 50–100 days. How much of this energy is converted into kinetic energy of the ejecta or escapes as SN luminosity in any given event remains highly uncertain. However, if, in the most optimistic scenario, all of the energy of the magnetar is emitted as radiation, only the (1–4) × 10¹⁴ G, 1 ms magnetars can power luminous SN after 50 days. Those with higher fields or longer periods lose energy much earlier and are no longer bright by the time their photons exit the ejecta. We therefore adopt the 4 × 10¹⁴ G, 1 ms magnetar as the fiducial case in our study, which has a total spin energy of ∼2 × 10⁵² erg.

Although it is thought that the spin-down energy of the magnetar is converted into high-energy photons, neutrinos, cosmic rays, and winds, how it is partitioned among them is not well understood and the spectrum of the photons themselves is not well-constrained by models. Furthermore, the efficiency with which photons, particles, and winds dynamically couple to local ejecta is not completely clear. Most calculations of magnetar-powered SLSN light curves assume that 100% of the dipole radiation goes into luminosity when in reality some goes into PdV work on the ejecta that changes its dynamics. We deposit all the energy of the magnetar into the ejecta as heat for simplicity. This approximation yields an upper limit to the violence of any fluid instabilities that could affect the light curve of the explosion because none of the energy of the magnetar directly escapes the flow as radiation.

It is not known for certain how much time is required for a proto-NS to evolve into a magnetar so we turn it on in KEPLER ∼100 s after the explosion, just after the shock has broken out of the surface of the star and the ejecta is expanding nearly homologously. This choice is reasonable, given that even if the magnetar had formed instantaneously it would have emitted at most 0.6% of its total spin energy in the first 100 s, which is negligible in comparison to the energy it later deposits in the flow or even the much lower energy of the original explosion.

We consider the same star as Chen et al. (2016), a 6 ${M}_{\odot }$ carbon–oxygen (CO) core that evolved from a 24 ${M}_{\odot }$ zero-age main-sequence star that loses its hydrogen envelope and part of its helium envelope prior to explosion (Sukhbold & Woosley 2014). The core has a radius of 2 × 10¹¹ cm and ${}^{12}{\rm{C}}$ and ${}^{16}{\rm{O}}$ mass fractions of 0.14 and 0.86, respectively. This star evolves through carbon, neon, oxygen, and silicon burning, building up a 1.45 ${M}_{\odot }$ iron core prior to collapse. We use a CO star because many SLSNe have been found to be Type I, implying that the progenitor has lost its outer layers to a stellar wind or binary interaction. It was evolved until its iron core reached collapse velocities of ∼1000 km s⁻¹, at which point the SN was initialized by the injection of linear momentum in its core.

CC SNe are not emergent features of KEPLER simulations because it lacks the physics and dimensionality required to capture core bounce and the heating of the atmosphere of the proto-NS by neutrinos, both of which are thought to launch the SN shock. We instead trigger the explosion by injecting linear momentum into the core at a predetermined mass coordinate with an energy corresponding to the energy chosen for the explosion. Here, the CO star explodes as a CC SN with an energy of 1.2 × 10⁵¹ erg, producing a 1.45 ${M}_{\odot }$ NS (magnetar) and ejecting 0.22 ${M}_{\odot }$ of ${}^{56}\mathrm{Ni}$. KEPLER includes the gravity due to the NS and the self-gravity of the ejecta.

Multidimensional simulations of the evolution of the cores of massive stars have shown that convective mixing of the oxygen and silicon layers occurs prior to the collapse of the iron core (Meakin & Arnett 2007). The mechanics of the explosion can be sensitive to this convective structure. Furthermore, numerous studies of core collapse and bounce have found that fluid instabilities and mixing occur even at these early stages of CC SNe (Burrows et al. 1995; Janka & Mueller 1996; Mezzacappa et al. 1998; Woosley & Janka 2005; Kifonidis et al. 2006). Even on short timescales of several ms, accretion shock instabilities along with neutrino heating produce strong mixing and break the spherical symmetry of the core. KEPLER cannot reproduce any of this mixing in 1D from first principles so we artificially mix some of the ${}^{28}\mathrm{Si}$ core using the mixing length theory in KEPLER.

As noted above, the magnetar is turned on 100 s after linear momentum is injected into the core of the star to trigger the SN. Its luminosity is uniformly deposited as heat in the innermost 10 Lagrangian zones of the mesh, which enclose the central 0.12 ${M}_{\odot }$ of the ejecta. These zones have an approximately constant energy generation rate per gram. The total energy deposited in the ejecta evolves over time as shown in the plot in Figure 1(c) corresponding to the magnetic field and period of our magnetar. Density and velocity profiles of the ejecta 110 s after the explosion are shown in Figure 2.

Zoom In Zoom Out Reset image size

The velocity of the forward shock is ∼2 × 10⁹ cm s⁻¹ and the innermost ejecta is expanding at ∼10⁸ cm s⁻¹. Density and velocity jumps appear at ∼10¹⁰ cm as luminosity from the magnetar heats the inner ejecta and it expands, plowing up gas and forming a dense shell that moves supersonically with respect to the ejecta. If allowed to evolve further in 1D, the shell would reach Δ_S ∼ 1000–2000 within a few thousand seconds after the explosion (see Figure 4 of Chen et al. 2016). We instead halt the simulation and port these profiles into CASTRO because all nuclear burning due to energy from the magnetar has ceased but the shell is not yet prone to hydrodynamical instabilities.

We map 1D spherically symmetric KEPLER blast profiles onto 3D cartesian grids in CASTRO with a conservative scheme developed by Chen et al. (2013) that preserves the mass, momentum, and energy of the ejecta. CASTRO is a multidimensional Adaptive Mesh Refinement (AMR) hydrodynamics code (Almgren et al. 2010; Zhang et al. 2011) with an unsplit piecewise-parabolic method hydro scheme (Colella & Woodward 1984), multi-species advection, and several equations of state (EOSs). At early stages of the simulation we use the Helmholtz EOS, which includes electron and positron pairs of arbitrary degeneracy, ions, and radiation (Timmes & Swesty 2000), to model the partly degenerate CO core but later switch to an ideal gas EOS after the density of the expanding ejecta falls below 10⁻¹² ${\rm{g}}\,{\mathrm{cm}}^{-3}$. We evolve mass fractions for ${}^{1}{\rm{H}}$, ${}^{4}\mathrm{He}$, ${}^{12}{\rm{C}}$, ${}^{16}{\rm{O}}$, ${}^{24}\mathrm{Mg}$, ${}^{28}\mathrm{Si}$, ${}^{40}\mathrm{Ca}$, ${}^{44}\mathrm{Ti}$, ${}^{48}\mathrm{Cr}$, and ${}^{56}\mathrm{Ni}$ to study how they are mixed during the explosion but do not employ a reaction network because nuclear burning due to heating by the magnetar has ceased. We neglect heating due to radioactive decay of ${}^{56}\mathrm{Ni}$ (and its conversion into iron) because the total energy released in our models is less than 1% of that of the SN. The gravities due to the ejecta, its envelope, and the compact object are all included in our models. The self-gravity of the ejecta is calculated from the monopole approximation by constructing a spherically symmetric gravitational potential from the radial average of the density and then determining the corresponding gravitational force everywhere in the AMR hierarchy. This approximation is reasonable, given the globally spherical structure of the supernova ejecta. The magnetar is treated as a point source.

We add a circumstellar medium (CSM) to the outer boundary of the ejecta due to the loss of the H and He layers of the star prior to the explosion. The envelope can be ejected by gravity waves (Quataert & Shiode 2012; Shiode & Quataert 2014) or super-Eddington stellar winds (Owocki et al. 2017, 2019) before the death of the star. This measure also prevents superluminal velocities when the shock crashes out of the star into very diffuse densities because CASTRO lacks relativistic hydrodynamics. Shiode & Quataert (2014) found that such ejections can create wind profiles with total masses of 10⁻³–1 ${M}_{\odot }$ extending out to 10¹⁵–10¹⁶ cm. We therefore surround the CO core with a ${\rho }_{r}={\rho }_{0}{(r/{r}_{0})}^{-2}$ wind profile, where ρ₀ is the density at the outer boundary of the ejecta, which we set to 2.11 × 10⁻¹¹ ${\rm{g}}\,{\mathrm{cm}}^{-3}$ at 1.75 × 10¹¹ cm. ρ₀ is a factor of 10⁻⁴ lower than the density at the surface of the star and represents the abrupt drop in density between the star and the CSM. The total mass of the envelope is 0.81 ${M}_{\odot }$.

Our CASTRO root grid is 5 × 10¹² cm on a side with a resolution of 128³ and eight nested grids centered on the star, which has a radius of 1.75 × 10¹¹ cm. This hierarchy of grids has a maximum initial spatial resolution of 1.53 × 10⁸ cm to resolve the ejecta in which energy from the magnetar is deposited and fluid instabilities later form. After the launch of the simulation we allow up to four additional levels of refinement, which are triggered by gradients in density, velocity, and pressure. If the SN shock reaches the grid boundary we double the size of the mesh while holding the number of mesh points constant and conservatively map the flows onto this new grid using the approach by Chen et al. (2013). In our run the original box was doubled 12 times to a final size of 2.05 × 10¹⁶ cm. Our maximum spatial resolution is similar to the finest zone size in Chen et al. (2016), ∼1.2 × 10⁸ cm, so it is sufficient to capture fluid instabilities and mixing. Outflow boundary conditions were set on the grid. The simulations are evolved until the shock reaches ∼10¹⁶ cm about 200 days after the explosion.

The magnetar luminosity that was turned on in KEPLER is continued in CASTRO by depositing heat in a central sphere that has a radius of ∼2% that of the expanding ejecta. This sphere is always resolved by ∼100 zones. We sinusoidally vary the energy we deposit in the sphere in angle by an amplitude of 1% to approximate anisotropies in emission from the NS. If densities in the cavity formed by the expansion of the newly heated ejecta become too small, energy deposition can lead to superluminal velocities because CASTRO is not a relativistic code. To prevent this we add small amounts of mass with the energy that total ∼2.92 × 10⁻¹⁰ ${M}_{\odot }$ by the end of the run, which is negligible compared with the mass of ejecta. This mass injection also limits the maximum velocity of the wind to ∼30% speed of the light.

The sphere into which the magnetar initially deposits its energy has a radius of 3 × 10⁹ cm. At a luminosity of ∼10⁴⁸ erg s⁻¹ it deposits more than 1 × 10⁵¹ erg into the flows in the first 1000 s. This rapid injection of energy causes the ejecta around the magnetar to quickly expand, forming a hot bubble. Figure 3 shows densities in the bubble and surrounding ejecta 720 s after the formation of the magnetar. The radii of the bubble and forward shock marked by the white dashed circles are found by spherically averaging the densities and velocity flows and identifying the appropriate peaks, as in Figure 2. The forward shock of the SN has grown to ∼10¹² cm. Rayleigh–Taylor (RT) fluid instabilities have formed on two spatial scales, just behind the forward shock and at the boundary of the magnetar bubble. The first are driven by the forward shock and the second are driven by wind from the magnetar. Only the former are still visible because the latter have devolved into turbulent hot gas in the bubble, as seen in the eddy-like structures within the inner dashed circle of Figure 3. Both instabilities efficiently mix the ejecta.

Zoom In Zoom Out Reset image size

We plot 1D angle-averaged profiles of density and velocity in Figure 4. The shell driven outwards by the magnetar bubble is visible as the density bump at r ∼ 3 × 10¹¹ cm, and it moves at a velocity of 5 × 10⁸ cm s⁻¹. The anisotropies in the flows driven by luminosity from the magnetar at early times seen in the right panel of Figure 3 lead to elongation in the bubble at much later times, as we discuss below.

Zoom In Zoom Out Reset image size

We show densities and velocities 30 h after magnetar formation in Figure 5, after it has injected about 90% of its spin-down energy into the bubble and its luminosity has fallen from 10⁴⁸ to 10⁴⁶ erg s⁻¹. Because of the rapid release of nearly 20 times the energy of the original explosion, the bubble, which is at 1.25 × 10¹⁴ cm, has almost overtaken the forward shock of the SN, which is at 1.6 × 10¹⁴ cm. Remnants of the RT fingers due to the forward shock at 720 s are now visible as the diffuse structures enclosed by the bubble at r ∼ (0.5–1.1) × 10¹⁴ cm. The fingers at r ∼ (1.3–1.6) × 10¹⁴ cm are new RT instabilities due to the deceleration of the forward shock as it plows up the CSM. Their amplitudes and velocities have been strongly enhanced by pressure from the hot bubble behind them, which has accelerated iron-group elements in them to speeds of (1–2) × 10⁹ cm s⁻¹, comparable to those of the forward shock. This process could be the origin of the high-velocity ${}^{56}\mathrm{Ni}$ and ${}^{56}\mathrm{Fe}$ observed at unexpectedly early times in some CC SNe, such as SN 1987A.

Zoom In Zoom Out Reset image size

Past simulations of SN 1987A by Fryxell et al. (1991), Kifonidis et al. (2006), Joggerst et al. (2010), and Hammer et al. (2010) showed that significant mixing can occur during the explosion that partially explained its observed spectra but could not produce its high-velocity ${}^{56}\mathrm{Ni}$ and ${}^{56}\mathrm{Fe}$ lines. More recent work by Mao et al. (2015) and Wongwathanarat et al. (2015) produced high-velocity ${}^{56}\mathrm{Ni}$ features but only under special conditions that depended on the choice of progenitor star and explosion energy. Ono et al. (2020) proposed a binary progenitor model and jet-like explosions for the origin of such lines that were quite different than the progenitor inferred for SN 1987A.

These models failed to reproduce high-velocity ${}^{56}\mathrm{Ni}$ and ${}^{56}\mathrm{Fe}$ lines because mixing in them is mainly driven by the reverse shock before the forward shock breaks out the star. The violence of the mixing depends on the structure of the progenitor and is greatest in red supergiants, whose extended envelopes allow instabilities to grow for longer times, dredging material deep in the ejecta up to higher velocities. However, observations indicate that the progenitor of 1987A was probably a blue supergiant in which mixing is expected to be weak because the star is more compact. Furthermore, the velocities observed for ${}^{56}\mathrm{Ni}$-rich ejecta are usually ∼4000–6000 km s⁻¹, which cannot be achieved by mixing due to the reverse shock (which decelerates rather than accelerates over time). Energy from a magnetar or pulsar inside a CC SN can naturally explain the high-velocity ${}^{56}\mathrm{Ni}$ and ${}^{56}\mathrm{Fe}$ found in some explosions such as SN 1987A.

By 200 days the magnetar bubble and forward shock are at ∼3 × 10¹⁵ cm and 10¹⁶ cm, respectively. Approximately 99% of the energy of the magnetar has been released and the ejecta is now expanding homologously. As shown in Figure 6, RT instabilities appear on the largest scales just behind the forward shock while the interior of the magnetar bubble has become elongated. This asymmetry was seeded by turbulence in the bubble at early times, when hot gas along some lines of sight preferentially expanded through channels of lower density in the ejecta. A somewhat boxlike structure enclosing the wind bubble is visible (and more prominent in the right panel of Figure 6). It arises because there is less AMR refinement in the low densities surrounding the bubble because of the relatively small density gradients there. This structure becomes more distorted as the grid is doubled and is therefore ultimately an artifact of the mesh. The morphology of the magnetar bubble begins to freeze in mass coordinate at this time while the RT instabilities in the outer regions can continue to grow, depending on the CMS.

Zoom In Zoom Out Reset image size

We show how some of the elements have become mixed at this time in Figure 7. The Fe represents any remaining ${}^{56}\mathrm{Ni}$ that was synthesized during the explosion and the ${}^{56}\mathrm{Co}$ and ${}^{56}\mathrm{Fe}$ into which the rest decayed. Most of the mixing has occurred between ${}^{12}{\rm{C}}$ and ${}^{16}{\rm{O}}$ in the outer regions and ${}^{24}\mathrm{Mg}$, ${}^{28}\mathrm{Si}$, ${}^{40}\mathrm{Ca}$, and Fe in the bubble, but some Fe has been dredged up to ∼8 × 10¹⁵ cm. 1D angle-averaged profiles of density, velocity, and elemental abundances are shown in Figure 8. The flat mass fractions at small radii suggest that the inner ejecta is approaching chemical homogeneity. Fe at mass fractions of 10⁻³–10⁻² has high velocities, ∼10⁹ cm s⁻¹. The velocity profile reveals two free expansions, that of the central bubble (which continues to be pumped with energy at low levels by the magnetar) and that of the surrounding ejecta. The ripples at the outer edge of the profile are due to the RT fingers. The radius of the inner free expansion (∼1.5 × 10¹⁴ cm) is much smaller than that of the bubble because it is the gas most directly exposed to radiation from the now much dimmer magnetar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Chen et al. (2016) modeled the same magnetar studied here except in 2D and only at early times in its evolution. They omit the CSM of the star because they only evolved the explosion for short times after shock breakout. They found that mixing in the shell of the magnetar bubble is more efficient in 2D, with many thin filamentary structures forming in contrast to the smoother clumpy structures found in 3D. This primarily is due to differences in turbulent cascades in 2D and 3D.

In the absence of a CSM, when radiation can escape the ejecta at early times, mixing in the magnetar bubble could determine the spectrum and light curve of the event because it sets the temperatures and structures of the clumps emitting high-energy photons from the explosion. However, most of these events are expected to occur in a dense envelope so radiation breakout would happen much later. As in our 3D models, iron-group elements deep in the ejecta were accelerated to high velocities at early times by energy from the magnetar in Chen et al. (2016), reaching speeds of 1.2 × 10⁹ cm s⁻¹.