Magnetar-Powered Supernovae in Two Dimensions. II. Broad-Line Supernovae Ic

TL;DR

This study uses two-dimensional simulations to explore how magnetar-powered supernovae, especially those associated with broad-line Type Ic supernovae, are influenced by the central engine, affecting nucleosynthesis, explosion morphology, and luminosity.

Contribution

It provides new insights into the effects of different central engine types on supernova explosion characteristics and nucleosynthetic yields through detailed hydrodynamical simulations.

Findings

Less than 0.05 solar masses of Ni are produced, regardless of engine type.

Explosion morphology can diagnose properties of the central engine.

Without circumstellar medium, these supernovae are relatively dim with peak magnitude around -16.5.

Abstract

Nascent neutron stars with millisecond periods and magnetic fields in excess of $10^{16}$ Gauss can drive highly energetic and asymmetric explosions known as magnetar-powered supernovae. These exotic explosions are one theoretical interpretation for supernovae Ic-BL which are sometimes associated with long gamma-ray bursts. Twisted magnetic field lines extract the rotational energy of the neutron star and release it as a disk wind or a jet with energies greater than 10$^{52}$ erg over $\sim 20$ sec. What fractions of the energy of the central engine go into the wind and the jet remain unclear. We have performed two-dimensional hydrodynamical simulations of magnetar-powered supernovae (SNe) driven by disk winds and jets with the CASTRO code to investigate the effect of the central engine on nucleosynthetic yields, mixing, and light curves. We find that these explosions synthesize less than 0.05 Msun of Ni and that this mass is not very sensitive to central engine type. The morphology of the explosion can provide a powerful diagnostic of the properties of the central engine. In the absence of a circumstellar medium these events are not very luminous, with peak bolometric magnitudes $M_b \sim -16.5 $ due to low Ni production.

Tables

Table 1. Table 1 Magnetar-powered Explosions

Model	Type	Engine Type	${}_{}{}^{56}\mathrm{Ni}$ Mass	Remnant	$M_{ej}$
			$M_{\odot }$		$M_{\odot }$
A	J	J(${2.5}^{\circ }$, $ϵ$)	$0.016$	$4.54M_{\odot }$ BH	5.46
B	W	W($ϵ$)	0.038	$1.66M_{\odot }$ NS	8.34
C	J+W	J(${2.5}^{\circ }$, 0.1$ϵ$), W(0.9$ϵ$)	0.036	$1.66M_{\odot }$ NS	8.34
D	J+W	J(${10}^{\circ }$, 0.1$ϵ$), W(0.9$ϵ$)	0.037	$1.66M_{\odot }$ NS	8.34

Table 2. Table 2 Ni 56 superscript Ni 56 {}^{56}\mathrm{Ni}  masses

Energy Injection Zone	Mode A	Mode B	Mode C	Mode D
	$M_{\odot }$	$M_{\odot }$	$M_{\odot }$	$M_{\odot }$
2000 - 3000 km	0.024	0.036	0.037	0.036
2000 - 4000 km	0.015	0.040	0.039	0.040
2000 - 5000 km	0.016	0.038	0.036	0.037

Equations

$$E=\frac{1}{2}I{\omega }^2\approx 2\times 10^{52}{P}_{\mathrm{ms}}^{-2}\mathrm{erg},$$

$$\begin{array}{rl}\frac{dE}{dt}& =-\frac{32{\pi }^4}{3c^3}(B{R}_n^3\sin \alpha {)}^2{P}^{-4}\\ & \approx 10^{49}{B}_{15}^2{P}_{\mathrm{ms}}^{-4}\mathrm{erg}{\mathrm{s}}^{-1},\end{array}$$

$$\begin{gathered} {\tau }_d\approx E/|\frac{dE}{dt}| \\ \approx 2000B_{15}^{-2}{P}_{\mathrm{ms}}^2\sec . \end{gathered}$$

$$t_{dyn}=\frac{446}{\sqrt{{\rho }_{\mathrm{c}}}}\sec ,$$

$$T_{\mathrm{s}}(r)=1.33\times 10^{10}{\left(\frac{E_{\exp }}{10^{51}\mathrm{erg}}\right)}^{1/4}{\left(\frac{r}{10^8\mathrm{cm}}\right)}^{-3/4}\mathrm{K},$$

$$t_{\mathrm{diff}}\sim 14{\kappa }^{\frac{1}{2}}{M}_{\mathrm{ej}}^{\frac{3}{4}}{E}^{-\frac{1}{4}}\mathrm{days}$$

Full text

**affiliationtext: EACOA Fellow, email: ken.chen@nao.ac.jp

1. Introduction

Most stars from 30 - 80 ${M}_{\odot}$ eventually collapse to black holes because the energy released by core collapse cannot drive a shock that is powerful enough to overcome the ram pressure of infall, so core bounce fails to produce an explosion (e.g. Sukhbold et al., 2015; Ertl et al., 2016). But this picture can change with rapidly rotating stars, in which a neutron star (NS) with a period of a few milliseconds may be born. Rotation can amplify the magnetic field of the NS above $10^{15}$ G, creating a magnetar. The magnetar might spin down quickly by magnetic braking and release its rotational energy in the form of a radiatively-dominated disk wind (Duncan & Thompson, 1992; Thompson & Duncan, 1993). During a brief, early phase of braking, the radiation can be in the form of x-rays and soft gamma-rays (Kouveliotou et al., 1998; Gaensler et al., 2005; Woosley & Bloom, 2006; Maeda et al., 2007; Mereghetti, 2008; Esposito et al., 2009). In some cases, if magnetorotational instabilities arise they can launch a collimated jet that pierces the outer layers of the star and produces a gamma-ray burst (GRB) (e..g. LeBlanc & Wilson, 1970; Burrows et al., 2007; Uzdensky & MacFadyen, 2006; Mösta et al., 2014).

If both a disk wind and jet are present, a highly asymmetric SN explosion may accompany the burst (e..g. Metzger et al., 2011; Soker, 2016). Such events release energies of up to $\sim 10^{52}{\mathrm{erg}}$, 10 times those of conventional core-collapse (CC) SNe (Paczyński, 1998; Iwamoto et al., 1998). These magnetar-powered SNe are likely observed as broad-line Type Ic SNe (SNe Ic-BL), which are often referred to as hypernovae, and they are among the most energetic explosions in the Universe (Smidt et al., 2014). SNe Ic-BL have a very broad absorption lines of oxygen and iron but lacking of helium and hydrogen in their spectra. Their light curves (LCs) peak at absolute magnitude $\approx$ -18 $-$ -20 mag and the shape of LC is different from that of SNe Ic, which shows a broad peak and a slow tail (Iwamoto et al., 2000). The ejecta of SNe Ic-BL expand at velocities about $30,000$ km/s which is much faster than that of normal SNe (Vink et al., 2015). About $0.05-0.3{{M}_{\odot}}$ ${}^{56}\mathrm{Ni}$  is estimated to form in SNe Ic-BL (Cano, 2013; Prentice et al., 2016). Since their explosion engines are closely tied to their central remnants, SNe Ic-BL are promising candidates for studying the physics of compact objects because they may account for the SNe associated with some GRBs (Kulkarni et al., 1998; Patat et al., 2001; Chornock et al., 2010). Less extreme ($10^{14}-10^{15}$ G, $5-10$ ms) magnetars may explain superluminous SNe (SLSNe Woosley, 2010; Kasen & Bildsten, 2010; Chen et al., 2016), because a substantial fraction of the total rotational energy of the neutron star is emitted as light at late time. Chen et al. (2016) used 2D simulations to study the mangetar-powered SNe by neutron stars of a constant magnetic field strength of $4\times 10^{14}$ G, with initial rotational periods of 1 ms and 5 ms. They found that fluid instabilities cause a strong mixing and fracture shells of ejecta into filamentary structures. The observational signatures of the resulting supernova could be very different from those predicted by the 1D models.

Predicting the observational signatures of magnetar-powered SNe is a key to properly identifying them as more are discovered by the new SN factories such as the Palomar Transient Factory (PTF; Law et al., 2009), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Kaiser et al., 2002) and the Large Synoptic Survey Telescope (LSST; Ivezic et al., 2008). Multidimensional simulations that bridge a large range of spatial scales are needed to model the light curves of magnetar-powered SNe because of their inherent asymmetry. Previous studies have mainly focused on just the physics of the central engines of these explosions (e.g., Burrows et al., 2007) but larger-scale simulations of mixing in jet-powered SNe (Couch et al., 2011; Papish & Soker, 2014a, b) and neutrino-wind driven SNe (Joggerst et al., 2010; Nordhaus et al., 2010; Wongwathanarat et al., 2015) have now been done. In this paper we study magnetar-powered SNe Ic-BL driven by central engines that are combinations of jets and disk winds. Our two-dimensional models include nucleosynthesis but not radiation transport, and they evolve the explosion from early internal mixing to breakout and then homologous expansion in order to calculate their observational signatures. In Section 2 we describe our progenitor models, explosion simulations and light curve calculations. The evolution of the explosions, including mixing, is examined in Section 3 and we discuss our results and conclude in Section 4.

2. Numerical Setup

We take as our initial conditions the collapsing carbon-oxygen core of a massive star that has been evolved from central carbon ignition to the onset of core collapse. At this stage, profiles for the star are mapped into the CASTRO code and exploded with a variety of central engines. Blast profiles from CASTRO are then evolved with the radiation hydrodynamics code STELLA to obtain light curves for these events.

2.1. KEPLER Progenitor Model

The progenitor star in all our models is a non-rotating 10 ${M}_{\odot}$ carbon-oxygen core with initial ${}^{12}\mathrm{C}$ and ${}^{16}\mathrm{O}$ mass fractions of 0.172 and 0.828, respectively. This progenitor, taken from another simulation campaign (Sukhbold & Woosley, 2014), roughly corresponds to the CO core of a non-rotating, solar-metallicity 35 ${M}_{\odot}$ star prior to the loss of its H and He envelope. It is evolved from central carbon ignition until the beginning of collapse with the 1D stellar evolution code KEPLER (Weaver et al., 1978). The actual mass of the zero-age main sequence star can vary depending on mass loss, rotation, magnetic field, and other physical parameters since more than one combination of progenitor mass and physical parameters can lead to the same CO core mass.

We consider a stripped CO core because the SNe Ic-BL observed to date are likely the explosions of massive stars that have shed their helium envelopes (Iwamoto et al., 1998; Nakamura et al., 2001). Such SNe must either be born with rapid rotation rates (Yoon et al., 2006) or be spun up at late times in a common envelope phase with another star (tidal locking) (Fryer et al., 1999; Izzard et al., 2004) or compact object (Zhang & Fryer, 2001). Either case usually results in the expulsion of the hydrogen and helium envelopes in some type of outburst (a luminous blue variable star, or LBV; Baraffe et al., 2001) or the ejection of a dense shell just before the death of the star. However, a large mass ejected from a massive star may reduce its angular momentum and prevent it from forming a rapidly rotating core, which is required for magnetar formation. An alternative path is quasi-chemically homogeneous evolution, which results in a large CO core that retains its angular momentum because there is no massive ejection. Tidal locking in the case of close binaries could then provide additional angular momentum to create a rapidly rotating core.

Our model was evolved for 14,100 years to form a 1.51 ${M}_{\odot}$ iron core with a radius of 1570 km. The 1D stellar evolution model was halted when any region of the star began to collapse faster than 1000 km s*-1*. At this point it was mapped into a 2D grid in CASTRO. The composition and velocity profiles of the pre-SN core are shown in Fig. 1.

2.2. 2D CASTRO Models

CASTRO is a multidimensional adaptive mesh refinement (AMR) radiation hydrodynamics code (Almgren et al., 2010; Zhang et al., 2011) with an unsplit piecewise parabolic method hydro scheme (Colella & Woodward, 1984). CASTRO has a Helmholtz equation of state based on Timmes & Swesty (2000), which has relativistic electron-positron pairs of arbitrary degeneracy, ions (which are treated as an ideal gas) and photons. We advect 17 isotopes, ${}^{1}\mathrm{H}$ , ${}^{3}\mathrm{He}$, ${}^{4}\mathrm{He}$, ${}^{12}\mathrm{C}$, ${}^{14}\mathrm{N}$, ${}^{16}\mathrm{O}$, ${}^{20}\mathrm{Ne}$, ${}^{24}\mathrm{Mg}$, ${}^{28}\mathrm{Si}$, ${}^{32}\mathrm{S}$, ${}^{36}\mathrm{Ar}$, ${}^{40}\mathrm{Ca}$, ${}^{44}\mathrm{Ti}$, ${}^{48}\mathrm{Cr}$, ${}^{52}\mathrm{Fe}$, ${}^{54}\mathrm{Fe}$, ${}^{56}\mathrm{Ni}$, and use a simplified prescription for nuclear burning that takes all elements to be burned to ${}^{56}\mathrm{Ni}$ when the gas temperature and density exceeds $4.3\times 10^{9}$ K and $1\times 10^{6}\mathrm{g}\,\mathrm{cm}^{-3}$. Densities, velocities, temperatures and mass fractions from KEPLER are mapped onto a 2D cylindrical grid in CASTRO with the conservative scheme of Chen et al. (2013), which preserves energies and masses over a large range of spatial scales. The mapping is done just before the formation of the NS, when infall velocities near the core reach $\sim 1000$ km s*-1*.

We only simulate a quadrant of the star in 2D, so the mesh is $2\times 10^{12}$ cm in both $r$ and $z$, or about twenty times the radius of the star, $r_{*}\sim 1.1\times 10^{11}$ cm. The star is shrouded by a low-density envelope $\rho\propto r^{-3.1}$ that prevents mixing as the forward shock plows up the circumstellar medium (CSM) after breakout. The root grid has 2562 zones and up to eight levels of refinement for an additional factor of up to 256 (28) in spatial resolution. The grid is refined on gradients in density, velocity, and pressure. This approach provides an effective simulation domain of $65,536\times 65,536$ zones.

The explosion energy from the magnetar is injected by hand. We center eight nested grids on the site of energy injection, each of which has twice the resolution of the grid above it for a maximum resolution equal to that of the lowest level of the AMR hierarchy. Reflecting and outflow boundary conditions are imposed on the inner and outer boundaries in both $r$ and $z$, respectively. We use a monopole approximation for self-gravity, in which the gravitational potential is constructed from the radial average of the density and used to calculate gravity forces everywhere in the AMR hierarchy by linear interpolation. This approximation is efficient and valid because the star is nearly spherically symmetric. Point source gravity due to the compact remnant is also included in our models.

We estimate the explosion energy and its timescale as follows. The radius of the NS is assumed to be $R_{n}\sim 10^{6}$ cm and the moment of inertia for a typical NS is $I\sim 10^{45}$ g cm2 so its initial rotational energy is

[TABLE]

where $P_{ms}$ is the period of the magnetar in milliseconds. This energy can be released through dipole radiation,

[TABLE]

where $B_{15}=B/10^{15}$G and assume $\sin\alpha=1$ for simplicity. The spin-down time scale can be approximated as

[TABLE]

If a magnetar forms with a rotation period of 1 ms and magnetic field stress of $10^{16}$ G a total energy of $\sim 2\times 10^{52}$ erg can be deposited into the surrounding core in 20 sec in the form of a very energetic disk wind or a collimated jet after the formation of the NS. Since the specifics of the central engine are not known, we consider four cases by varying the fraction of the energy that goes into the wind (isotropic) or the jet (anisotropic). Both observational evidence (Gaensler & Slane, 2006; Younes et al., 2016) and theoretical studies (Proga et al., 1998; Proga, 2000) of the disk wind from the accretion disk suggest that such wind can be highly inhomogeneous and anisotropic. Because our simulation cannot resolve the central disk, yet do not include MHD. Therefore, the mixing from our spherical disk wind model can be the lower limit of actual mixing.

Engine A is a purely jet-driven explosion, with all the energy going into a jet with a half opening angle $\phi=2.5^{\circ}$. This model is a practical example mentioned in (Gilkis et al., 2016) before. Engine B is a wind driven explosion in which all the energy of the explosion is deposited isotropically into the surrounding core. Engine C is both a jet and a wind in which 90% of the energy is isotropic and 10% goes into a jet with $\phi=2.5^{\circ}$. Engine D is the same as Engine C but with $\phi=10^{\circ}$. Engines A and B are the two extreme cases while Engine C is motivated by the conventional GRB SN and Engine D represents the case of a wobbling jet. The precession is caused by kink instabilities which may often occur in core-collapse jets in (Bromberg & Tchekhovskoy, 2016). Such instabilities disperse the energy in a larger opening angle, so the jet may die before reaching the stellar surface without producing a GRB. Gilkis (2016) also suggested that a wide jet can form if there are non-axisymmetric patterns in the core-collapse. We summarize the four engines in Table 1. These models are associated with the jet-feedback mechanism (JFM) of CCSNe (Papish & Soker, 2011; López-Cámara et al., 2013). In JFM scenario, the jets launched from the central compact object must be fast and narrow, they would deposit their energy inside the star through shock waves, then forming two hot bubbles to push out infalling gas. Therefore, the jet becomes slow, massive, and wide. If the jet feedback is effective, accretion would halt on early and result in regular supernova explosions. Otherwise, the accretion lasts longer and it supply more energy to the jet and eventually creates much energetic explosions (Soker, 2016).

We assume that the wind and jet deposit their energy uniformly in a region 2,000 - 5,000 km from the center of the star. In the wind, all of the energy is deposited as internal energy of the gas. In the jet, $50\%$ of the energy goes into internal energy and the rest goes into kinetic energy by injecting a highly relativistic momentum flux with a speed of $1.5\times 10^{10}$ cm s*-1* along the polar axis. In some cases the jet can blow out so much gas that extremely low densities result, which can cause numerical difficulties in the runs. To prevent blowout from creating a complete vacuum in the energy injection region we added a total mass of $5\times 10^{-8}$ ${M}_{\odot}$ to these zones, which is completely negligible in comparison to the mass of star.

Gas falling into the central 2,000 km of the grid is assumed to accrete onto the NS or black hole (BH). 364 grids at the deepest level of refinement resolve the injection region during the simulation. In each case, the shock is evolved until it reaches the outer grid boundary at $\sim$ 20 $r_{*}$, when the ejecta are expanding homologously. Due to the large number of levels of refinement, the four models took about 800,000 CPU hours on $\it Hopper$ at the National Energy Research Scientific Computing Center (NERSC).

2.3. STELLA

We calculate light curves for our explosions with the 1D multigroup radiation hydrodynamics code STELLA (Blinnikov & Bartunov, 1993; Blinnikov et al., 2006; Blinnikov & Tolstov, 2011). STELLA can calculate spectral energy distributions (SEDs) for the blast profiles at each time step. Multicolor LCs can be obtained by convolving filter functions with the SEDs. All our SN light curves are calculated with 100 frequency bins from 1 to $5\times 10^{4}$ Å on a log scale. STELLA implicitly evolves time-dependent equations of the angular moments of the intensity averaged over a frequency bin. Local thermodynamic equilibrium is assumed in determining the ionization states of materials. STELLA has been extensively used for modeling SN light curves (Blinnikov et al., 2000; Chugai et al., 2004; Woosley et al., 2007; Tominaga et al., 2011; Moriya et al., 2011).

3. Magnetar Powered Explosions

The central engine injects $10^{51}$ erg/s into a shell of $r=2000-5000$ km for 20 sec, heating the gas and producing a strong shock that reverses the collapse. The energy is evenly distributed through the shell, and in models B, C, and D the heated gas reverses infall in the shell in under a second. Most of the ${}^{56}\mathrm{Ni}$ forms in the first 0.5 sec of collapse, when the gas is heated to $4.3\times 10^{9}$ K and compressed to $\rho>10^{6}$ $\mathrm{g}\,\mathrm{cm}^{-3}$ by energy from the magnetar below and by infall from above. In models B - D, most of the ${}^{56}\mathrm{Ni}$ is formed by compression heating and later expelled by the magnetar wind. 0.026 - 0.04 ${M}_{\odot}$ of ${}^{56}\mathrm{Ni}$ are ejected by the explosion, which drives a strong shock from the core at 1 - 4 $\times 10^{9}$ cm s*-1*. Energy from the magnetar primarily governs the dynamics of the explosion, not energy from nuclear burning.

${}^{56}\mathrm{Ni}$ masses are listed in Table 1 for all four models. The isotropic models create more ${}^{56}\mathrm{Ni}$ than the jet model, in which material is burned in a much smaller solid angle. Nevertheless, the jet produces more ${}^{56}\mathrm{Ni}$ than its small solid angle alone might suggest because it dredges material up from greater depths that would otherwise have fallen onto the compact remnant. Models B, C and D produce about the same ${}^{56}\mathrm{Ni}$ mass, 0.037 ${M}_{\odot}$, which is significantly less than would be expected for a CC SN for reasons we discuss in greater detail below.

3.1. ${}^{56}\mathrm{Ni}$ Production

Distributing the energy of the magnetar evenly from radii of 2000 - 5000 km is reasonable but ad hoc. To determine how sensitive ${}^{56}\mathrm{Ni}$ production is to the thickness of the injection site we repeated models A - D at higher resolution on much smaller meshes. The grid was reduced to $10^{10}$ cm along both axes, with $512^{2}$ zones and no AMR. We considered three injection regions: 2000 - 3000 km, 2000 - 4000 km and 2000 - 5000 km. Although energy from the magnetar is deposited for 20 sec, densities at the injection site fall very quickly so no ${}^{56}\mathrm{Ni}$ will be made after 0.5 sec after the explosion. We therefore only evolve the SN for the first 5 sec in these models. The outcomes of all twelve runs are summarized in Table 2.

${}^{56}\mathrm{Ni}$ masses are sensitive to the structure of the progenitor for two reasons. First, the structure determines how much mass is burned to ${}^{56}\mathrm{Ni}$ by the magnetar wind and by compression heating from collapse. It also determines how much of this ${}^{56}\mathrm{Ni}$ then falls back onto the compact remnant. We show the structure of the presupernova core in Fig. 2. Gas below $r\sim 1.57\times 10^{8}$ cm (the radius of the iron core) will collapse directly to a proto-neutron star. Most of the ${}^{56}\mathrm{Ni}$ forms at radii of $\sim 1.57-2.38\times 10^{8}$ cm, where densities and temperatures are driven to $5.7-25\times 10^{6}$ $\mathrm{g}\,\mathrm{cm}^{-3}$ and $3.07-4.64\times 10^{9}$ K by the wind and by collapse. If the magnetar forms promptly, the gas in this shell would be burned to 0.21 ${M}_{\odot}$ of ${}^{56}\mathrm{Ni}$, about what would be expected from a CC SN.

But in most scenarios the magnetar does not turn on immediately after proto-neutron star formation. If the delay is just 0.2 sec, ${}^{56}\mathrm{Ni}$ below $r\sim 1.85\times 10^{8}$ cm would fall back onto the NS and only 0.05 - 0.12 ${M}_{\odot}$ would be ejected. If the delay is 0.5 sec, ${}^{56}\mathrm{Ni}$ at $r\leq 2.0\times 10^{8}$ cm will fall back and even less ${}^{56}\mathrm{Ni}$ will be ejected, 0.01 - 0.07 ${M}_{\odot}$. Our choice of energy deposition at $r\geq 2\times 10^{8}$ cm implies a 0.5 sec delay in the explosion, which is why on average only 0.037 ${M}_{\odot}$ of ${}^{56}\mathrm{Ni}$ is produced in models B

• D.

The reason why ${}^{56}\mathrm{Ni}$ production is not sensitive to the thickness of the injection region is that no ${}^{56}\mathrm{Ni}$ is formed at $r>$ 2300 km in our models. We again find that the jet in model A dredges up more ${}^{56}\mathrm{Ni}$ from lower radii than it forms on its own. The energy was injected in a half-opening angle of only $2.5^{\circ}$. If all the mass within this solid angle was burned to ${}^{56}\mathrm{Ni}$ in the $r=2000-3000/4000$ km tests, it would only make $1.65\times 10^{-3}$ ${M}_{\odot}$, not the 0.015 - 0.016 ${M}_{\odot}$ it actually dredges up.

The extra material dredged up by the jet in model A is due to lateral pressure forces that drive gas sideways out of the jet when energy is suddenly injected into its solid angle. The gas that is driven sideways also experiences shear forces upward from gas deeper within the solid angle that is mostly expanding outward along the axis of the jet, and the combination of the two motions is what brings up the extra material. In the Model A runs, the simulation with magnetar energy injection at $r=2000-3000$ km produces 0.024 ${M}_{\odot}$ of ${}^{56}\mathrm{Ni}$, nearly twice that of the other two sites, because the lateral pressure forces are much stronger when the injected energy is concentrated in the thinner shell, and they drive more vigorous sideways mixing.

Our choice of injection site is consistent with simple analytic estimates of the radius out to which material can be burned to ${}^{56}\mathrm{Ni}$ in the core. The hydrodynamic time of the core is

[TABLE]

where $\rho_{\rm{c}}$ is its mean density. If we take $\rho_{\rm{c}}=10^{6}$ $\mathrm{g}\,\mathrm{cm}^{-3}$, then $t_{dyn}\sim 0.5$ sec, which is when most explosive burning and ${}^{56}\mathrm{Ni}$ production would happen. Over this time the magnetar injects $0.5\times 10^{51}$ $\mathrm{erg}$ into the outer layers of the core, and most of this energy goes into explosive burning. Except at small radii near the origin of the shock, the peak temperature at radius $r$ can be obtained by setting $(4/3)r^{3}aT_{\rm{s}}^{4}\sim E_{\rm{exp}}$, where $T_{\rm{s}}$ is the shock temperature and $a$ is the radiation density constant. This equation assumes that the heat capacity of the material behind the shock is dominated by the radiation field, and that expansion and pressure waves behind the shock are capable of maintaining nearly isothermal conditions there. The shock temperature at radius $r$ can then be expressed as (Woosley et al., 2002)

[TABLE]

In our simulation, $E_{\rm{exp}}=0.5\times 10^{51}$, so a spherical volume with $r=2.91\times 10^{8}{\mathrm{cm}}$ will reach $T_{s}\sim 4.5\times 10^{9}$ K for ${}^{56}\mathrm{Ni}$ synthesis, which is consistent with our CASTRO models.

3.2. Accretion and Breakout

The jet in model A gradually broadens as it propagates through the star at $\sim$ 15% of the speed of light, breaking through its surface at 19 sec. The shock driven by the wind in model D breaks out at 73 sec and the shock driven by both a wind and jet breaks out of the star at 38 and 57 sec in models C and D. In model A, a small fraction of the 0.016 ${M}_{\odot}$ of the ${}^{56}\mathrm{Ni}$ ejected by the jet is present in the bow shock. Since the jet propagates mostly along the axis of the star, accretion continues inward along the equatorial plane. We show accretion rates and blowout for model A in Fig. 3 and Fig. 4. We calculate the mass accretion by assuming the mass flux across a surface at radius of $r=2\times 10^{8}$ cm. However, this accretion is subject to fluid instabilities. At $t=50$ sec, two data points are missing because accretion flux becomes an outflow due to the strong fluid instabilities before the accretion shuts down. Accretion persists for 52 sec in model A until it is completely shut down by the cocoon of hot gas from the jet. The rate varies from 0.01 - 0.3 ${M}_{\odot}$ s*-1* and the total accreted mass is 3.03 ${M}_{\odot}$, enough to collapse the NS to a 4.54 ${M}_{\odot}$ black hole (BH). From Fig. 3, we estimate the BH forms about 10 sec after accretion starts when the mass of central compact object exceeds 2 ${M}_{\odot}$. Energy injection by the jet may become less efficient after the BH forms, but exactly how the NS collapses to a BH is beyond the scope of our study. While we just assume that the mass that falls back onto the NS turns it into a BH, angular momentum transport can delay BH formation and radiation from the BH may then reduce accretion. In Fig. 4 the core of the star has been completely blown out along the axis of the jet but remains intact along the equatorial plane, and likely falls back onto the BH at later times.

3.3. Mixing

Although all four models have the same injected energy, the geometry of injection leads to a range of mixing during the explosion. We show the four SNe at shock breakout in Fig. 5. In model A the jet produces a strong collimated outflow in the polar direction. The flow exhibits knots and kinks due to jet instabilities and some shear instabilities, but the latter do not grow to sufficient amplitudes to produce much mixing before breakout. Some heavy elements deep in the star are dredged up by the jet. Model B exhibits much more mixing, on par with what would be expected for a CC SN. The collision between the wind and the collapsing outer layers of the star drives the dynamical instabilities. The formation of a reverse shock also drives Rayleigh-Taylor (RT) insatiabilities and contributes to mixing significantly. In models C and D, twisting in the jet and the disk wind leads to both shear and RT instabilities, as seen in the elongated structures in the ejecta. There are traces of ${}^{56}\mathrm{Ni}$ in the jet, so gamma-ray emission from ${}^{56}\mathrm{Ni}$ is possible at early phases of the explosion. In all four cases the star becomes unbound by 80 sec after the explosion, but mixing continues after breakout.

The CSM around massive stars is usually diffuse because their strong winds and large ionizing UV fluxes drive away gas in their vicinity. The shock therefore accelerates rapidly when it breaks out of the star. The dramatic drop in density between the surface of the star and the CSM sometimes crashed the hydro solver in our simulations (a problem that has been reported by others; see, e.g., Whalen et al., 2013a, b). To prevent numerical difficulties we take the CSM to fall off from the surface of the star as $\rho\propto r^{-3.1}$. This diffuse envelope also prevents the formation of reverse shocks in the ejecta that can cause additional mixing after breakout so we can study just the mixing that is intrinsic to the explosion itself. In model A the velocity of the jet at breakout is $\sim$ 7 $\times$ 109 cm s*-1*, or $\sim$ 25% of the speed of light, so it is mildly relativistic. The breakout velocity is sensitive to the CSM density so if it falls off dramatically the shock can accelerate to nearly the speed of light and produce a GRB, although this does not happen in our models.

We run our simulations until the shock or jet reaches $\sim$ 20 $r_{*}$ and mixing ceases. At this stage the ejecta are expanding homologously and their energy is almost entirely kinetic (their internal energy is $\sim$ 10% of the total energy). We show snapshots of the ejecta for all four models at this stage in Fig. 6. The jet continues to broaden well after breakout. The original opening angle of $5^{\circ}$ has grown to $46^{\circ}$ as it partially thermalizes with gas in its path, forming a shock that is perpendicular to the jet. Models B and D have similar morphologies. Jets with large opening angles soon become hard to distinguish from explosions driven just by disk winds because the energy of the jet is rapidly spread throughout the ejecta.

We summarize mixing in all four models in Fig. 7. The metals mostly trace the outline of the jet in model A, although there is clearly some mixing. A ring pattern created by the cocoon of the jet is also visible. There is very little mixing along the equatorial plane of the star. There is much more mixing in models B, C, and D. The shells of elements that build up inside the star over its life are completely disrupted by the shock. Models B and D, which are somewhat difficult to distinguish in density, are more easily differentiated by isotope distribution, which can provide a diagnostic of central engine type which can be examined in SN spectra.

4. Light Curves

To calculate light curves for these explosions we average the flow variables at each radius over a sample of 10 angles on the grid. We then map these angle-averaged radial profiles onto 1D spherical grids in the STELLA code. The blast profiles are taken at the onset of homologous expansion. We only calculate light curves for models B, C, and D because the jet in model A is highly directional and relativistic so a different technique would be required to compute its light curve. We evolve all three explosions out to 60 days.

Bolometric light curves for the three models are shown in Fig. 8. Each exhibits a brief, extremely luminous pulse due to shock breakout that lasts for about an hour. The flux then dims and later rebrightens as photons from the radioactive decay of ${}^{56}\mathrm{Ni}$ begin to diffuse out of the ejecta, which occurs on a timescale $t_{\mathrm{diff}}$ (Arnett, 1996), where

[TABLE]

and $\kappa$ is the opacity of the ejecta in units of the Thompson scattering opacity of free electrons in hydrogen-free gas (0.2 cm2 g*-1*), $M_{\mathrm{ej}}$ is the mass of the ejecta through which photons must diffuse to reach the surface in units of solar mass and $E$ is the energy of the explosion in units of 1051 erg. Taking $\kappa=$ 1, $M_{\mathrm{ej}}=$ 4 because of the partial dredging up of ${}^{56}\mathrm{Ni}$, and $E=$ 20, we find that $t_{\mathrm{diff}}\sim$ 18 days, which corresponds to the time at which ${}^{56}\mathrm{Ni}$ rebrightening peaks in all three magnetar-powered SNe. ${}^{56}\mathrm{Co}$ decay follows, and the light curves gradually fade on timescales comparable to its half-life, $\sim$ 77 days.

Because less than 0.05 ${M}_{\odot}$ of ${}^{56}\mathrm{Ni}$ is made, bolometric luminosities peak at $M_{b}\approx-16.5$, which is not much brighter than normal CC SNe. Mixing may dredge ${}^{56}\mathrm{Ni}$ up to nearly the surface of the ejecta in some models. We show ${}^{56}\mathrm{Ni}$ mass fractions and velocities at the end of all four runs in Fig. 9. Except in model B, which is driven by only a wind, a clump of ${}^{56}\mathrm{Ni}$ has broken through the ejecta along the axis of the jet. The ${}^{56}\mathrm{Ni}$ rich ejecta have velocities of $1-3\times 10^{9}$ cm s*-1* that are large enough to Doppler shift their spectra. Note the roughly evenly spaced velocity contours, which indicate homologous expansion of the ejecta.

Model A could produce a GRB without a luminous SN component. If the jet is sustained by the accretion of gas falling in from the equator, it could drive a GRB from 20 s to 100 sec. This model could explain GRB 060614 (Gal-Yam et al., 2006; Della Valle et al., 2006), which synthesized very little ${}^{56}\mathrm{Ni}$, $5\times 10^{-4}$ ${M}_{\odot}$, lasted $\sim 100$ s, and did not produce an optically luminous SN ($M_{v}>-12.3$ mag).

5. Summary and Conclusions

We have performed 2D simulations of magnetar-powered SNe with the CASTRO AMR code. These transients are presumed to be observed as SNe Ic-BL. We examine central engines in the form of disk winds, jets, and combinations thereof by varying the morphology of energy injection in the core of the star. These engines lead to a variety of mixing in the explosion before and after shock breakout. We find that infall from the equatorial plane can create a BH if there is a collimated jet; otherwise, a NS forms.

Although magnetar-powered SNe can be ten times more energetic than normal CC SNe, as is observed with SNe Ic-BL, they do not necessarily make superluminous transients. Light curves for our SNe Ic-BL manifest a sharp, very luminous transient due to shock breakout followed by dimming and then mild ${}^{56}\mathrm{Ni}$ rebrightening at 10 - 20 days. Our SNe Ic-BL are not very luminous because they only make 0.02 - 0.05 ${M}_{\odot}$ of ${}^{56}\mathrm{Ni}$, which is consistent with models by Nishimura et al. (2015) and Suwa & Tominaga (2015). The ${}^{56}\mathrm{Ni}$ production is small because most of it is made by compression heating during the core collapse instead of energy injection from the magnetar, and much of it falls onto the nascent NS before it can be expelled by the magnetar wind. Unless energy is deposited at a very high rate, the magnetar is unlikely to produce enough ${}^{56}\mathrm{Ni}$ to rebrighten its light curve. Suwa & Tominaga (2015) suggests that magnetars with magnetic fields $>2\times 10^{16}$ G and rotation periods $<$ 1 ms (energy deposition rates $>4\times 10^{51}$ erg s*-1*) can produce $>$ 0.2 ${M}_{\odot}$ ${}^{56}\mathrm{Ni}$.

However, extreme magnetars such as these may release a significant fraction of their rotational energy as gravitational waves rather than dipole radiation, so their ${}^{56}\mathrm{Ni}$ yields may still be low. But these events could become much brighter if their ejecta crash into a dense CSM (e.g., Moriya et al., 2013; Whalen et al., 2013c). Furthermore, if a successful CC explosion is followed by the formation of a magnetar with magnetic fields of 1014 - 1015 G and rotation periods $\sim$ 5 ms, it may become very bright at later times (Woosley, 2010; Kasen & Bildsten, 2010; Chen et al., 2016). ${}^{56}\mathrm{Ni}$ blown out by jets at high velocities might also emit observable gamma rays.

In future models the nuclear reaction network could be improved to calculate better yields for magnetar-powered SNe. More realistic initial conditions for the core collapse engine based on dedicated 2D or 3D simulations rather than energy injection by hand could also be implemented. Although we do not include radiation transport our results indicate that multidimensional radiation hydrodynamics will be required to determine how photons are emitted from dense structures created by fluid instabilities (Chen et al., 2014, 2016) and produce more accurate light curves. More realistic prescriptions for the CSM can also be calculated from radiation hydrodynamical models of the ambient H II region and wind cavity of the star. In such profiles, the jet might reach velocities of 0.95$c$ and produce a GRB, mandating special relativistic upgrades to the hydro solver in CASTRO. Another scenario worth consideration is the collision of the asymmetric ejecta with a CSM.

Our models will soon be confronted by more detections of magnetar-powered SNe in SN searches with PTF, Pan-STARRS, LSST and future searches with Euclid and the Wide-Field Infrared Survey Telescope (WFIRST). These exotic explosions may have occurred more frequently in the primeval universe because the Population III initial mass function (IMF) may have been top-heavy (e.g., Bromm et al., 2009; Whalen, 2012; Glover, 2013). Many primordial stars may also have been born with rotation speeds close to the breakup limit (Stacy et al., 2011, 2013) or in binaries (Turk et al., 2009; Stacy et al., 2010; Stacy & Bromm, 2013). Some of these highly energetic transients (e..g. Nakauchi et al., 2012) may also be found in the near infrared by the James Webb Space Telescope and ground-based 30m telescopes and probe the properties of the first stars in the Universe.

We thank the anonymous referee, whose comments and suggestions improved the quality of this paper. The authors also thank Ann Almgren, Weiqun Zhang, and Sergei Blinnikov for technical support with CASTRO and STELLA. K.C. acknowledges the support of an EACOA Fellowship from the East Asian Core Observatories Association and the hospitality of the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293. Work at UCSC was supported by an IAU-Gruber Fellowship, the DOE HEP Program (DE-SC0010676) and the NASA Theory Program (NNX14AH34G). D.J.W. was supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013) via the ERC Advanced Grant ”STARLIGHT: Formation of the First Stars” (project number 339177). YS was supported in part by the Grant-in-Aid for Scientific Research (Nos. 16K17665 and 16H00869). V.B. acknowledges support from NSF grant AST-1413501. CASTRO was developed through the DOE SciDAC program by grants DE-AC02-05CH11231 and DE- FC02-09ER41618. Our numerical simulations were performed at NERSC and the Center for Computational Astrophysics (CfCA) at the National Astronomical Observatory of Japan (NAOJ).