GRAVITATIONAL WAVE SIGNATURES OF MAGNETOHYDRODYNAMICALLY DRIVEN CORE-COLLAPSE SUPERNOVA EXPLOSIONS

Successful detection of neutrinos from SN1987A paved the way for the neutrino astronomy (Hirata et al. 1987; Bionta et al. 1987) alternative to conventional astronomy by electromagnetic waves. Core-collapse supernovae are now expected to be opening yet another branch of astronomy, gravitational-wave astronomy. Currently, long baseline laser interferometers LIGO (Abbott et al. 2005), VIRGO,³ GEO600,⁴ TAMA300 (Ando & the TAMA Collaboration 2005), and AIGO⁵ with their international networks of observatories are beginning to take data at sensitivities where astrophysical events are predicted (see, e.g., Hough et al. 2005 for a recent review). For these detectors, core-collapse supernovae have been proposed as one of the most plausible sources of gravitational waves (GWs; see, e.g., Kotake et al. 2006; Ott 2009 for recent reviews).

Although the explosion mechanism of core-collapse supernovae has not yet been completely clarified, current multi-dimensional simulations based on refined numerical models show several promising scenarios. Among the candidates is the neutrino heating mechanism aided by convection and standing accretion shock instability (SASI; e.g., Marek & Janka 2009; Bruenn et al. 2009; Scheck et al. 2004, 2008; Suwa et al. 2010), the acoustic mechanism (Burrows et al. 2006, 2007b), and the magnetohydrodynamic (MHD) mechanism (e.g., Ardeljan et al. 2000; Kotake et al. 2004a, 2005; Obergaulinger et al. 2006b; Burrows et al. 2007a; Takiwaki et al. 2009, and references therein). The explosion geometry is expected to be unipolar and bipolar for the former two cases, and bipolar for the MHD mechanism. Since the GW signatures imprint live information of the asphericity at the moment of explosion, they are expected to provide us an important clue to the solution of the supernova mechanism.

So far, most of the theoretical predictions of GWs have focused on the bounce signals in the context of rotational core collapse (e.g., Mönchmeyer et al. 1991; Zwerger & Mueller 1997; Kotake et al. 2003b; Shibata & Sekiguchi 2004; Ott et al. 2004, 2007a, 2007b; Dimmelmeier et al. 2002, 2007, 2008; Scheidegger et al. 2008). For the bounce signals to have a strong and characteristic signature, the iron core must rotate rapidly enough. The waveforms are categorized into the following three types, namely types I, II, and III. Type II and III waveforms have been shown to be less likely to appear than type I waveforms because a combination of general relativity (GR) and electron capture near core bounce suppresses multiple bounces in the type II waveforms (Dimmelmeier et al. 2007; Ott et al. 2007a, 2007b). In general, a realistic nuclear equation of state (EOS) is stiff enough to prohibit type III waveforms. After bounce, asymmetries due to convection (Burrows & Hayes 1996; Müler & Janka 1997; Fryer 2004; Müller et al. 2004), SASI (Marek et al. 2009; Kotake et al. 2007, 2009a, 2009b; Murphy et al. 2009), and g-mode oscillations of protoneutron stars (Ott et al. 2006) are expected to account for sizable GW signals. In general, detection of these GW signals in the postbounce phase (except for the g-mode oscillation) is far more difficult than the bounce signals because they do not possess the same clear signature, but change stochastically with time as a result of chaotically growing convection and SASI in the nonlinear hydrodynamics (Kotake et al. 2009b, 2011; Marek et al. 2009; Murphy et al. 2009).

Rapid rotation, necessary for strong bounce signals, is likely to be obtained in only ∼1% of the massive star population (e.g., Woosley & Bloom 2006). However, this can really be the case for progenitors of rapidly rotating metal-poor stars, which experience so-called chemically homogeneous evolution (Woosley & Heger 2006; Yoon & Langer 2005). The high angular momentum of the core and the strong precollapse magnetic field are preconditioned for the MHD mechanism because the MHD mechanism relies on the extraction of rotational free energy of the collapsing core via magnetic fields. The energetic MHD explosions have received great attention recently as to their possible relevance to magnetars and collapsars (e.g., Harikae et al. 2009, 2010 for collective references), which are presumably linked to the formation of long-duration gamma-ray bursts (GRBs; e.g., Meszaros 2006).

Among the previous studies mentioned above, only a small portion of studies has been spent on determining GW signals with the MHD mechanism (Yamada & Sawai 2004; Kotake et al. 2004b; Obergaulinger et al. 2006b; Cerdá-Durán et al. 2007; Shibata et al. 2006; Scheidegger et al. 2010). This may be because the MHD effects on the dynamics as well as their influence over the GW signals can be visible only for cores with precollapse magnetic fields over B₀ ≳ 10¹² G (Obergaulinger et al. 2006b; Kotake et al. 2004b). Considering that the typical magnetic field strength of GRB progenitors is at most ∼10¹¹–10¹² G (Woosley & Heger 2006), this is already an extreme situation. Interestingly, in a more extreme case of B₀ ∼ 10¹³ G, a secularly growing feature in the waveforms was observed (Obergaulinger et al. 2006a; Shibata et al. 2006; Scheidegger et al. 2010). Moreover, Obergaulinger et al. (2006a) identified as type IV a waveform in which quasi-periodic large-scale oscillations of GWs near bounce are replaced by higher frequency irregular oscillations. Some of these MHD simulations follow adiabatic core collapse in which a polytropic EOS is employed to mimic supernova microphysics. At this level of approximation, the bounce shock generally does not stall and a prompt explosion occurs within a few tens of milliseconds after bounce. Therefore, the main focus in these previous studies has been rather limited to the early postbounce phase (≲ several 10 ms). However, for models with weaker precollapse magnetic fields akin to the current GRB progenitors, the prompt shocks first stall in the core like a conventional supernova model with more sophisticated neutrino treatment (e.g., Burrows et al. 2007a). In such a case, the onset of MHD explosions, depending on the initial rotation rates, can be delayed until ∼100 ms after bounce (Burrows et al. 2007a; Takiwaki et al. 2009). There remains room to study GW signatures in such a case, which we hope to do in this work.

In this study, we choose to take precollapse magnetic fields less than 10¹² G based on a recent GRB-oriented progenitor models. By this choice, it generally takes much longer after bounce compared with the adiabatic MHD models to amplify magnetic fields strong enough to overwhelm the ram pressure of the accreting matter, leading to MHD-driven explosions. Even if the speed of jets in MHD explosions is only mildly relativistic, Newtonian simulations are not numerically stable because the Alfvén velocity ($\propto B/\sqrt{\rho }$) could exceed the speed of light unphysically, especially when the strongly magnetized jets (i.e., large B) propagate to a stellar envelope with decreasing density (ρ). To follow a long-term postbounce evolution numerically stably, we perform special relativistic MHD (SRMHD) simulations (Takiwaki et al. 2009) in which a realistic EOS is employed and the neutrino cooling is taken into account via a multiflavor neutrino leakage scheme. Note in our previous study of GWs in magneto-rotational core collapse (Kotake et al. 2004b) that we were unable to study properties of the GWs in the postbounce phase because the employed Newtonian simulations quite often crashed, especially in the case of strong MHD explosions. To include GR effects in this study, we follow a prescription in Obergaulinger et al. (2006a), which is reported to capture basic features of full GR simulations quite well. By changing precollapse magnetic fields as well as initial angular momentum distributions parametrically, we compute 12 models. By doing so, we hope to study the properties of GWs in MHD explosions systematically and also address their detectability.

Section 2 presents descriptions of the initial models and numerical methods employed in this work. Formalism for calculating the gravitational waveforms in SRMHD is summarized in Section 3. The main results are given in Section 4. We summarize our results and discuss their implications in Section 5.

2.1. Initial Models

We make precollapse models by taking the profiles of density, internal energy, and electron fraction distribution from a rotating presupernova model of E25 in Heger & Langer (2000). This model has a mass of 25 M_☉ at the zero-age main sequence, however, it loses the hydrogen envelope and becomes a Wolf-Rayet (WR) star of 5.45 M_☉ before core collapse. Our computational domain involves the whole iron core of 1.69 M_☉. Note that this model is suggested as a candidate progenitor of long-duration GRBs because type Ib/c core-collapse supernovae that originated from WR stars have an observational association with long-duration GRBs (e.g., Woosley & Bloom 2006).

Since little is known about the spatial distributions of rotation and magnetic fields in evolved massive stars, we add the following profiles in a parametric manner to the non-rotating core mentioned above. For the rotation profile, we assume a cylindrical rotation of

$$\begin{equation} \Omega (X,Z) =\Omega _{0}\frac{X_{0}^2}{X^2+X_{0}^2} \frac{Z_{0}^4}{Z^4+Z_{0}^4}, \end{equation} \tag{ 1 }$$

where Ω is the angular velocity and X and Z denote the distance from the rotational axis and the equatorial plane, respectively. The parameter X₀ represents the degree of differential rotation, which we choose to change in the following three ways: 100 km (strongly differential rotation), 500 km (moderately differential rotation), and 2000 km (uniform rotation). The parameter Z₀ is fixed to 1000 km.

Regarding the precollapse magnetic field, we assume that the magnetic field is nearly uniform and parallel to the rotational axis in the core and dipolar outside. This can be modeled by the following effective vector potential:

$$\begin{equation} A_r=A_\theta =0, \end{equation} \tag{ 2 }$$

$$\begin{equation} A_\phi =\frac{B_0}{2}\frac{r_0^3}{r^3+r_0^3}r\sin \theta, \end{equation} \tag{ 3 }$$

where A_r, A_θ, A_ϕ are the vector potentials in the r, θ, and ϕ directions, respectively, r is radius, r₀ is radius of the core, and B₀ is a model constant (see Takiwaki et al. 2004 for detail). In this study, r₀ is set to 2000 km, which is approximately the size of the precollapse iron core.

By changing the initial angular momentum, degree of differential rotation, and strength of magnetic fields, we compute 12 models. The model parameters are shown in Table 1. The models are named with the first part, B12 or B11, representing the strength of the initial magnetic field, the second part, X1, X5, or X20, indicating the degree of differential rotation (X₀ = 100, 500, 2000 km, respectively), and the third part, β = 0.1 or 1, showing the rotation parameter β. Here, β represents the ratio of the rotational energy to the absolute value of the gravitational energy prior to core collapse. The original progenitor of model E25 in Heger & Langer (2000) has a uniform rotation profile in the iron core and the initial β parameter is ∼0.15%. Thus the initial angular momentum in our models of B11X20β0.1 and B12X20β0.1 is similar to the original one. In a GRB-oriented progenitor of model 35OB (Woosley & Heger 2006), the precollapse magnetic fields reach ∼10¹¹–10¹² G and β ∼ 0.2%, which is not significantly different from the parameters chosen here. Although little is known about the degree of precollapse differential rotation, an extremely strong one—for example, our model series of X1—is unrealistic. These models, albeit rather academic, are examined in order to clearly see the effects of differential rotation.

Table 1. Summary of Initial Models

	β(%)
	0.1%			1%
	X0(km)
B0	100 km	500 km	2000 km	100 km	500 km	2000 km
1011 G	B11X1β0.1	B11X5β0.1	B11X20β0.1	B11X1β1	B11X5β1	B11X20β1
1012 G	B12X1β0.1	B12X5β0.1	B12X20β0.1	B12X1β1	B12X5β1	B12X20β1

Notes. Model names are labeled by the precollapse magnetic fields and rotation. β represents the ratio of initial rotational energy to the absolute value of the initial gravitational energy. From left to right in the table, Ω₀ in unit of rad s⁻¹ (Equation (1)) is 24, 2.8, 0.95, 76, 8.9, and 3.0, respectively. Note that X₀ and B₀ are defined in Equations (1) and (3), respectively.

Download table as:  ASCIITypeset image

2.2. Special Relativistic Magnetohydrodynamics

Numerical results in this work are calculated by the SRMHD code developed in Takiwaki et al. (2009). In the following, we first mention the importance of SR and then briefly summarize the numerical schemes.

The Alfvén velocity of MHD jets propagating into the outer layer of the iron core can be estimated as v_A ∼ 10¹⁰ cm s⁻¹(B/10¹³ G)(ρ/10⁵ g cm⁻³)^−1/2, where ρ and B are the typical density and the magnetic field along the rotational axis. It can be readily inferred that the Alfvén velocity could exceed the speed of light unphysically in Newtonian simulations. SR corrections are also helpful in correctly capturing the dynamics of infalling material in the vicinity of the protoneutron star because their free-fall velocities and rotational velocities become close to the speed of light. Such conditions are quite ubiquitous in MHD explosions. We have learned that even if the propagation speeds of the jets are only mildly relativistic (at the least), SR treatments are quite important for keeping the stable numerical calculations in good accuracy over a long-term postbounce evolution (e.g., Harikae et al. 2009).

The MHD part of our code is based on the formalism of De Villiers et al. (2003). The state of the relativistic fluid element at each point in the space time is described by its density, ρ, specific energy, e, velocity, v_i, and pressure, p, and the magnetic field in the laboratory frame is described by the four-vector $\sqrt{4\pi }b^{\mu }={^*F}^{\mu \nu }U_{\nu }$, where *F^μν is the dual of the electro-magnetic field strength tensor and U_ν is the four-velocity. The basic equations of the SRMHD code are written as

$$\begin{equation} \frac{\partial {D}}{\partial {t}} +\frac{1}{\sqrt{\gamma }}\partial _{i}{\sqrt{\gamma }Dv^i} = 0 \end{equation} \tag{ 4 }$$

$$\begin{eqnarray} &&\frac{\partial {E}}{\partial {t}} +\frac{1}{\sqrt{\gamma }}\partial _{i}{\sqrt{\gamma }Ev^i} = -p\frac{\partial {W}}{\partial {t}} -\frac{p}{\sqrt{\gamma }}\partial _{i}{\sqrt{\gamma }W v^i} -{\cal L}_\nu \quad \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \frac{\partial {(S_i-b^tb_i)}}{\partial {t}} &+&\frac{1}{\sqrt{\gamma }} \partial _{j}{\sqrt{\gamma }(S_i v^j-b_ib^j)} = -\frac{1}{2} (\rho h (Wv_k)^2 \nonumber\\ && - (b_k)^2)\partial _{i}{\gamma ^{kk}} -\big(\rho h W^2- {b^t}^2\big) \partial _{i}{\Phi _\mathrm{tot}}\nonumber \\ & & -\partial _{i}{\left(p+\frac{|b|^2}{2}\right)} \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \frac{\partial {(Wb^i - Wb^tv^i)}}{\partial {t}} +\partial _{j}{(Wv^jb^i-Wv^ib^j)} &=&0 \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \partial ^{k}{{\partial _{k}{\Phi }}}&=&4\pi \biggl [ \rho h (W^2+(v_k)^2) + 2\left(p+\frac{\left|b\right|^2}{2}\right) \nonumber \\ && -((b^{0})^2+(b_{k})^2) \biggl ], \end{eqnarray} \tag{ 8 }$$

where $W={1}/{\sqrt{1-v^kv_k}}$, D = ρW, E = eW, and S_i = ρhW²v_i are the Lorentz boost factor, auxiliary variables corresponding to density, energy, and momentum, respectively. All of them are defined in the laboratory frame. Equations (4)–(6) represent the mass, energy, and momentum conservation, respectively. In Equation (6), note that the relativistic enthalpy, h = (1 + e/ρ + p/ρ + |b|²/ρ), includes the magnetic energy. Equation (7) is the induction equation in SR. In solving the equation, the method of characteristics is implemented to propagate accurately all modes of MHD waves (see Takiwaki et al. 2009 for more detail). Equation (8) is the Poisson equation for the gravitational potential of Φ, which is solved by the modified incomplete Cholesky conjugate gradient method. For an approximate treatment of GR gravity, Φ_tot in Equation (6) includes a GR correction to Φ as in Buras et al. (2006). We employ a realistic EOS based on the relativistic mean field theory (Shen et al. 1998).

We approximate the neutrino cooling by a neutrino multiflavor leakage scheme (Epstein & Pethick 1981; Rosswog & Liebendörfer 2003), in which three neutrino types, electron neutrino (ν_e), electron antineutrino ($\bar{\nu }_{e}$), and heavy-lepton neutrinos (ν_μ, $\bar{\nu }_{\mu }$, ν_τ, $\bar{\nu }_{\tau }$, collectively referred to as ν_X), are taken into account. The implemented neutrino reactions are electron capture on proton and free nuclei, positron capture on neutron, and photo-, pair, plasma processes (Fuller et al. 1985; Takahashi et al. 1978; Itoh et al. 1989, 1990). A transport equation for the lepton fractions (namely $Y_e-Y_{e^{+}},Y_{\nu _e}, Y_{\bar{\nu }_e}$, and $Y_{{\nu _{X}}})$ is solved in an operator-splitting manner (see Equation (7) in Takiwaki et al. 2009 for more details). $\cal {L}_{\nu }$ in Equation (5) represents the neutrino cooling rate summed over all the reactions, which can also be estimated by the leakage scheme (see Epstein & Pethick 1981; Rosswog & Liebendörfer 2003; Kotake et al. 2003a for more detail).

In our two-dimensional simulations, the spherical coordinates are employed with 300(r) × 60(θ) grid points to cover the computational domain assuming axial and equatorial symmetry. The radial grid is nonuniform, extending from 0 to 4000 km with a finer grid near the center. The finest grid for the radial direction is set to 1 km. The polar grid uniformly covers from θ = 0 to π/2. The finest grid for the polar direction is 25 m. The numerical tests and convergence with this choice of the numerical grid points are given in Section 5 of Takiwaki et al. (2009).

To measure the strength of explosion, we define the explosion energy as:

$$\begin{equation} E_{\mathrm{exp}} = \int _{\mathrm{D}} dV e_{\mathrm{local}} = \int _{\mathrm{D}} dV (e_{\mathrm{kin}}+e_{\mathrm{int}}+e_{\mathrm{mag}}+e_{\mathrm{grav}}). \end{equation} \tag{ 9 }$$

Here e_local is the sum of e_kin, e_int, e_mag, and e_grav which represent kinetic, internal, magnetic, and gravitational energy, respectively, and are defined as

$$\begin{eqnarray} e_{\mathrm{kin}}&=&\rho W \left(W-1\right), \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} e_{\mathrm{int}}&=& e W^2+p(W^2-1), \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} e_{\mathrm{mag}}&=& |b|^2\left(1-\frac{1}{2W^2}\right) -\frac{{b^0}^2}{2W^2}, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} e_{\mathrm{grav}}&=&-\rho h W^2 \Phi, \end{eqnarray} \tag{ 13 }$$

and D in Equation (9) represents the domain where the local energy (e_local) is positive, indicating that the matter is gravitationally unbound. The explosion energy is evaluated when the MHD jets pass through 1000 km along the polar direction.

To extract the gravitational waveforms in MHD explosions, we present the stress formula in SRMHD for later convenience. As shown below, this can be done straightforwardly by extending the Newtonian MHD formulation presented in Kotake et al. (2004b).

From the Einstein equation, one obtains the following formula as a primary expression for the leading part of the gravitational quadrupole field emitted by the motion of a fluid in the post-Newtonian approximation (e.g., Mönchmeyer et al. 1991; Finn & Evans 1990; Blanchet et al. 1990),

$$\begin{eqnarray} h_{ij}^{{\rm TT}}({\mbox{\boldmath $X$}},t) &=& \frac{4G}{c^4 R}P_{ijkl} ({\mbox{\boldmath $N$}}) \int d^3 x T_{kl}, \end{eqnarray} \tag{ 14 }$$

where G and c are the gravitational constant and the velocity of light, respectively, T_kl is the energy momentum tensor of the source, and $R \equiv (\delta _{ij} X^i X^j)^{1/2} = |{\mbox{\boldmath $X$}}|$ is the distance between the observer and the source. P_ijkl, with ${\mbox{\boldmath $N$}} = {\mbox{\boldmath $X$}}/R$, denotes the transverse-traceless (TT) projection operator on the plane orthogonal to the outgoing wave direction ${\mbox{\boldmath $N$}}$ (e.g., Mönchmeyer et al. 1991).

T_ij consists of the three parts, namely perfect fluid, electromagnetic field, and gravitational potential as follows:

$$\begin{eqnarray} T_{ij}&=& {T_{ij}}_{\rm (hyd)} + {T_{ij}}_{\rm (mag)} + {T_{ij}}_{\rm (grav)}. \end{eqnarray} \tag{ 15 }$$

The first term, which we refer as the hydrodynamic part, is explicitly written as

$$\begin{eqnarray} {T_{ij}}_{\rm (hyd)} &=& \rho _{*} W^2 v_i v_j + p \delta _{ij}, \end{eqnarray} \tag{ 16 }$$

where ρ_* is effective density defined as

$$\begin{equation} \rho _{*}= \rho + \frac{e + p + |b|^2}{c^2}. \end{equation} \tag{ 17 }$$

The second term in Equation (15) represents the contribution from the magnetic fields as

$$\begin{eqnarray} {T_{ij}}_{\rm (mag)} &=& -b_i b_j. \end{eqnarray} \tag{ 18 }$$

The last term in Equation (15) is the contribution from the gravitational potential,

$$\begin{eqnarray} {T_{ij}}_{\rm (grav)} &=& \frac{1}{4\pi G}\Biggl (\Phi _{,i}\Phi _{,j} - \frac{1}{2}\delta _{ij}\Phi _{,m}\Phi ^{,m}\Biggr), \end{eqnarray} \tag{ 19 }$$

where Φ corresponds to the self-gravity in Equation (8).

In our axisymmetric case, only one non-vanishing quadrupole term remains in the metric perturbation, namely ℓ = 2, m = 0 in terms of the pure-spin tensor harmonics, as

$$\begin{equation} h_{ij}^{{\rm TT}}({\mbox{\boldmath $X$}},t) \stackrel{\ell = 2, m = 0}{=} \frac{1}{R}A^{{\rm E2}}_{20}\left(t - \frac{R}{c}\right)T^{{\rm E2},20}_{ij}(\theta,\phi), \end{equation} \tag{ 20 }$$

where T^{E2, 20}_ij(θ, ϕ) is

$$\begin{equation} T^{{\rm E2},20}_{ij}(\theta,\phi) = \frac{1}{8}\sqrt{\frac{15}{\pi }}\sin ^2 \theta \end{equation} \tag{ 21 }$$

(e.g., Thorne 1980). The projection operator in Equation (14) acts on T_ij as

$$\begin{equation} P_{ijkl} T^{kl} = 2 T^{zz} - T^{xx} - T^{yy}. \end{equation} \tag{ 22 }$$

Transforming Equation (14) to the spherical coordinates, and expressing b_i and v_i in terms of unit vectors in the r, θ, ϕ directions, we obtain for A^E2₂₀ the expression

$$\begin{equation} {A_{20}^{\rm {E} 2}} = {A_{20}^{\rm {E} 2}}_{\rm (hyd)} + {A_{20}^{\rm {E} 2}}_{\rm (mag)} + {A_{20}^{\rm {E} 2}}_{\rm (grav)}, \end{equation} \tag{ 23 }$$

where

$$\begin{eqnarray} {A_{20}^{\rm {E} 2}}_{\rm (hyd)} &=& \frac{G}{c^4} \frac{32 \pi ^{3/2}}{\sqrt{15 }} \int _{0}^{1}d\mu \int _{0}^{\infty } r^2 \,dr {f_{20}^{\rm {E} 2}}_{\rm (hyd)}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} {f_{20}^{\rm {E} 2}}_{\rm (hyd)} &=& \rho _{*}W^2({v_r}^2(3 \mu ^2 -1) + {v_{\theta }}^2(2 - 3 \mu ^2) - {v_{\phi }}^{2}\nonumber\\ && - 6 v_{r} v_{\theta } \,\mu \sqrt{1-\mu ^2}); \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} {A_{20}^{\rm {E} 2}}_{\rm (grav)} &=& \frac{G}{c^4} \frac{32 \pi ^{3/2}}{\sqrt{15 }} \int _{0}^{1}d\mu \int _{0}^{\infty } r^2 \,dr {f_{20}^{\rm {E} 2}}_{\rm (grav)}, \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} {f_{20}^{\rm {E} 2}}_{\rm (grav)} &=& \Big[\rho h (W^2+(v_k/c)^2) + \frac{2}{c^2}\left(p+\frac{\left|b\right|^2}{2}\right) -\frac{1}{c^2}\big((b^{0})^2 \nonumber\\ && +(b_{k})^2 \big)\Big] \left[- r \partial _{r} \Phi (3 \mu ^2 -1) + 3 \partial _{\theta } \Phi \,\mu \sqrt{1-\mu ^2}\right];\nonumber\\ \end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray} {A_{20}^{\rm {E} 2}}_{\rm (mag)} &=& - \frac{G}{c^4} \frac{32 \pi ^{3/2}}{\sqrt{15 }} \int _{0}^{1}d\mu \int _{0}^{\infty } r^2 \,dr {f_{20}^{\rm {E} 2}}_{\rm (mag)}, \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} {f_{20}^{\rm {E} 2}}_{\rm (mag)} &=& \big[{b_r}^2(3 \mu ^2 -1) + {b_{\theta }}^2(2 - 3 \mu ^2) - {b_{\phi }}^{2}\nonumber\\ && - 6 b_{r} b_{\theta } \mu \sqrt{1-\mu ^2}\big], \end{eqnarray} \tag{ 29 }$$

where μ = cos θ. For later convenience, we write the total GW amplitude as

$$\begin{equation} h^{\rm TT} = h^{\rm TT}_{\rm (hyd)} + h^{\rm TT}_{\rm (mag)} + h^{{\rm TT}}_{\rm (grav)}, \end{equation} \tag{ 30 }$$

where the quantities of the right hand are defined by combining Equations (20) and (23) with Equations (25), (27), and (29). By dropping O(v/c) terms, the above formulae reduce to the conventional Newtonian stress formula (e.g., Mönchmeyer et al. 1991). In the following computations, the observer is assumed to be located in the equatorial plane (θ = π/2 in Equation (21)), and the source is assumed to be located near our Galactic center (R = 10 kpc).

The gravitational waveforms obtained in this work can be categorized into two types, which we call the increasing type and the cancellation type for convenience. Note that the latter type does not mean a new waveform, as will be explained later in this section. The former type was presented in previous literature (Obergaulinger et al. 2006a; Shibata et al. 2006; Scheidegger et al. 2010), however, its properties have not yet been clearly understood. In Section 4.1, we first overview their characters, which are peculiar in the case of MHD explosions. In Section 4.2, we analyze their properties by carefully comparing each contribution in Equation (30) to the total GW amplitudes. Then in Section 4.3, we perform the spectra analysis and discuss their detectability.

4.1. Properties of Waveforms in the MHD Exploding Models

Figure 1 shows examples of the two categories, which we refer to as the increasing type (left panels) and the cancellation type (right panels), respectively. In the increasing type, the total wave amplitudes (red line) have a monotonically increasing trend after bounce (t − t_b = 0 in the figures) while in the cancellation type (right panels), the total amplitudes after bounce stay almost zero. This is because the contribution from the magnetic fields (blue line, Equation (29)) cancels with the contribution from the sum of the hydrodynamic and gravitational parts (green line, Equations (25) and (27)). Regardless of the difference in the two types, it is common that the magnetic contribution (blue line) increases almost monotonically with time. Not surprisingly, the bounce GW signals (t − t_b ≲ 20 ms) are categorized into the so-called type I or II waveforms. Note here that the MHD simulations are terminated at around 100 ms after bounce for all the computed models. This is simply because the GW amplitudes in a later phase decrease because the MHD shock comes out of the computational domain and the enclosed mass in the domain becomes smaller.

Zoom In Zoom Out Reset image size

Figure 2 depicts a classification of the computed models, in which "C" and "I" indicate the cancellation and increasing type, respectively. I* indicates the mixture of the two types, which we call the intermediate type. The figure shows that the bifurcation of the two types is predominantly determined by the precollapse magnetic fields so that the models with stronger magnetic fields (B₀10¹² G) are basically classified as the increasing type. For models shown in bold, which are slow rotators with uniform rotation in our models, the field amplification works less efficiently than for models with stronger differential rotation (such as X₀ = 100, 500 km; see also Cerdá-Durán et al. 2008). This suppresses the increase in the wave amplitudes due to magnetic fields, which gives rise to the intermediate state between the two types.

Zoom In Zoom Out Reset image size

Figure 3 is a summary showing intervals measured from the stall of the bounce shock (top) and from core bounce (bottom) until the MHD-driven revival of the stalled bounce shock. As already mentioned in Takiwaki et al. (2009), generation of the postbounce MHD jets proceeds in the following two ways. One is launched relatively promptly after the stall of the bounce shock, typically earlier than ∼30 ms (see models in red) and another is launched rather later after bounce (≳ 30 ms; models in green). In this sense, our computed models could be roughly categorized into two types, namely the prompt MHD explosion (those in red) or the delayed MHD explosion (green). This simply reflects that it takes longer time for the weakly magnetized models to amplify the magnetic pressure behind the stalled shock enough strong to overwhelm the ram pressure of accreting matter. By comparing Figure 2 to Figure 3, the two characters in the waveforms have a rough correlation with the difference of the explosion dynamics.

Zoom In Zoom Out Reset image size

Figure 4 shows a summary regarding the explosion energy (e.g., Equation (9)). Comparing Figure 2 to Figure 4, the explosion energy for the increasing type (models in yellow in Figures 2 and 4) is higher compared to the cancellation type (models in light blue). Given the rotation rate of β = 0.1%, the explosion energy is smallest for the intermediate type (models in orange) because of the insufficient field amplification as mentioned above. It can also be shown that the explosion energies for the models with the increasing trend generally exceed 10⁵¹ erg (models in yellow).

Zoom In Zoom Out Reset image size

Having summarized the waveform classification together with the explosion dynamics, we move on to look in more detail at what causes the difference between the two types in the next section.

4.2. Analysis of Waveforms

4.2.1. Increasing Type

By taking model B12X1β0.1 as a reference, we first focus on the increasing type waveform. The left panel of Figure 5 depicts distributions of entropy (left half) and plasma β (right half, the ratio of matter to magnetic pressure) at 100 ms after bounce. It can be seen that the outgoing jets indicated by the velocity fields (arrows in the left half) are driven by the magnetic pressure behind the shock (see bluish region, i.e., low plasma β in the right half of the panel).

Zoom In Zoom Out Reset image size

The right panel of Figure 5 shows contributions to the total GW amplitudes (Equation (27)), in which the left-hand panels are for the sum of the hydrodynamic and gravitational part, namely log (± [f^E2_20(hyd) + f^E2_20(grav)]) (left top(+)/bottom(−)) (Equations (25) and (27)), and the right-hand panels are for the magnetic part, namely log (± f^E2_20(mag)) (right top(+)/bottom(−)) (e.g., Equation (29)). By comparing the top two panels in Figure 5, it can be seen that the positive contribution overlaps the regions where the MHD outflows exist. The major positive contribution is from the kinetic term of the MHD outflows with large radial velocities (e.g., +ρ_*W²v_r² in Equation (25)). The magnetic part also contributes to the positive trend (see the top right half in the right panel; labeled by mag(+)). This comes from the toroidal magnetic field (e.g., +b_ϕ² in Equation (29)), which dominantly contributes to driving MHD explosions. The magnetic contribution was already mentioned in Kotake et al. (2004b). This study furthermore adds that the kinetic energy of MHD outflows more importantly contributes to the positive trend.

Figure 6 shows a normalized cumulative contribution of each term in A^E2₂₀, which is estimated by the volume integral of A^E2₂₀ within a given sphere enclosed by a certain radius. It can be seen that the respective contributions of the hydrodynamic and gravity parts (indicated by "hyd and grav") are prominent for a radius outside ∼1000 km, which stems from exploding regions with large kinetic energy as mentioned above. Note that the secular drift observed in the increasing waveform may come from ambiguities in estimating the gravitational potential in the stress formula as pointed out by Dimmelmeier et al. (2002). To exclude such potential issues, we plot the waveform calculated by the first moment of momentum–density formalism (e.g., Finn & Evans 1990). As is shown in Figure 7, the increasing trend is also seen in the waveform estimated by the first-moment formalism (green line). As a side note, a smoother curve is obtained for the stress formula (red line), about which Mönchmeyer et al. (1991) pioneeringly mentioned that the numerical evaluation of the time derivative sometimes makes the waveform noisy.⁶

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.2. Cancellation Type

Now we proceed to focus on the cancellation-type waveform by taking model B11X1β1 as an example. The top panel in Figure 8 shows density distributions (left half) with the plasma β (right half). Note that the jet head of MHD outflows is located at ∼400 km along the pole which sticks out of the plots in Figure 8 (compare the difference in scales of Figure 5). In fact, the bluish regions (low β) around the pole have a lower density because material there has already been blown up due to the passage of MHD-driven shocks. The middle panel of Figure 8 shows the term-by-term contribution to A^E2₂₀ as in the right panel of Figure 5. From the right half of the panel, it is shown that the magnetic contribution is dominantly positive (labeled by mag(+)), which is also the case with the increasing type as mentioned in the previous section. Note that the negative contribution from magnetic fields (bluish region in the bottom right half of the middle panel) comes from the regions that are off-axis from the propagation of the MHD jets. In these regions, the poloidal components of magnetic fields are stronger than the toroidal ones, which makes the negative contribution mainly through −[b_r²(3μ² − 1) + b_θ²(2 − 3μ²)] in Equation (29).

Zoom In Zoom Out Reset image size

Looking at the sum of the hydrodynamic and gravitational part (left in the middle panel), a large negative contribution comes from regions near the rotational axis (in red, bottom left). The bottom panel of Figure 8 further shows a contribution from the hydrodynamic (left) and gravitational parts (right), separately. Regarding the hydrodynamic part, the negative contribution is highest in the vicinity of the equatorial plane which closely coincides with the oblately deformed protoneutron star (red in the bottom left in the bottom panel). This is because the negative contribution comes from the centrifugal forces (e.g., the term related to the rotational energy, −ρv_ϕ² in Equation (25)). For the gravitational part (right side in the bottom panel), a large negative contribution comes from regions in the vicinity of the rotational axis (bottom right). This comes from the term of −r∂_rΦ(3μ² − 1) in Equation (25), which is the radial gradient of the gravitational potential. Remembering that ∂_rΦ > 0 is generally satisfied in self-gravitating objects, the sign of this term is determined by μ = cos θ, which is a directional cosine measured from the rotational axis. As a result, the gravitational part makes a negatively contribution in the vicinity of the rotational axis (μ ∼ 1). These two factors coming both from the hydrodynamic and gravity part make the negative contribution to the GW amplitudes. In fact, the negative contribution is highest near the surface of the protoneutron star (∼10 km in Figure 9). Outside it, the magnetic part becomes almost comparable to the sum of the hydrodynamic and gravitational parts, which creates the cancellation type shown in the right panels of Figure 1.

Zoom In Zoom Out Reset image size

4.3. Spectrum Analysis

Now we move on to performing a spectral analysis (Figure 10). For the two types (left panel, increasing type; right panel, cancellation type), the peak amplitudes in the spectra are around 1 kHz, which comes from the GWs near at bounce. It can be seen that the spectra for lower frequency domains (below ∼100 Hz) are much larger for the increasing type (left panels) compared with the cancellation type (right panels). This reflects a slower temporal variation of the secular drift inherent to the increase-type waveforms (e.g., Figure 1).

Zoom In Zoom Out Reset image size

As a measure for characterizing the dominance in the lower frequency domains, we define $\tilde{h}_{\rm low}$, which represents average amplitudes below 100 Hz (see Figure 11). Although the peak amplitudes, $\tilde{h}_{\rm peak}$, in the spectra have no clear correlation with the two types, we point out that the final GW amplitudes (the first column) and the $\tilde{h}_{\rm low}$ (the third column) are much larger for the increasing type (in yellow) compared with the cancellation type (in light blue). In Figure 10, the peak amplitudes near 1 kHz are, irrespective of the two types, marginally within the detection limits of the currently running detector of the first LIGO and the detection seems more feasible for the next-generation detectors such as the LCGT and the advanced LIGO for a Galactic supernova. It is true that the GWs in the low frequency domains mentioned above are relatively difficult to detect because of seismic noises, but a recently proposed future space interferometer like the Fabry–Perot type DECIGO is designed to be sensitive in those frequency regimes (Kawamura et al. 2006; Kudoh et al. 2006). The sensitivity curve of the detector is plotted with the black line in Figure 10. Our results suggest that these low-frequency signals, if observed, could be one important messenger of the increasing type waveforms that are likely to be associated with MHD explosions exceeding 10⁵¹ erg.

Zoom In Zoom Out Reset image size

By performing a series of two-dimensional SRMHD simulations, we studied signatures of GWs in the MHD-driven core-collapse supernovae. In order to extract the gravitational waveforms, we presented a stress formula including contributions both from magnetic fields and special relativistic corrections. By changing the precollapse magnetic fields and initial angular momentum distributions parametrically, we computed 12 models. As for the microphysics, a realistic EOS was employed and the neutrino cooling was taken into account via a multitype neutrino leakage scheme. With these computations, we found that the total GW amplitudes show a monotonic increase after bounce for models with a strong precollapse magnetic field (10¹² G) that also have a rapid rotation imposed. We pointed out that this trend stems both from the kinetic contribution of MHD outflows with large radial velocities and from the magnetic contribution dominated by the toroidal magnetic fields that predominantly trigger MHD explosions. For models with weaker initial magnetic fields, the total GW amplitudes after bounce stay almost zero because the contribution from the magnetic fields cancels with the contribution from the hydrodynamic counterpart. These features can be clearly understood with a careful analysis of the explosion dynamics. It was pointed out that the GW signals with an increasing trend, possibly visible to the next-generation detectors for a Galactic supernova, would be associated with MHD explosions exceeding 10⁵¹ erg.

Although the presented simulations have utilized the leakage scheme to approximate the deleptonization, it would be more accurate (especially before bounce) to employ a formula developed by Liebendörfer (2005), which was designed to fit one-dimensional Boltzmann results. Figure 12 shows snapshots at around 33 ms after bounce for model B12X20β0.1 in which the leakage scheme (left panel) or the Y_e prescription (right panel) is employed, respectively. Note that "G15" is taken in our simulation until bounce among the parameter sets in (Liebendörfer 2005)⁷ and is switched to the leakage scheme at the postbounce phase. As is shown, shock revival also occurs for the model with the Y_e prescription (right panel). Figure 13 depicts the magnetic pressure (red line) versus ram pressure (blue line) along the polar axis (left panel) or the equatorial plane (middle) near the rebirth of the stalled shock for the model with the Y_e prescription (top panels) or the leakage scheme (bottom panels). For the equator, the magnetic pressure is much less than the ram pressure (middle panel), while the magnetic pressure amplified by the field wrapping along the pole becomes as high as the ram pressure of the infalling material at the shock front, leading to the MHD-driven shock formation (see right panels). Regardless of the two different deleptonization schemes, these important features associated with the MHD explosions are shown to be quite similar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now we mention a comparison between the obtained results and relevant MHD simulations. Model R4E1CF in Scheidegger et al. (2010), whose precollapse rotational parameter is β = 0.5% with a uniform rotation imposed and whose initial poloidal magnetic field is set to 10¹² G, is close to our model B12X20β0.1. From their Figure 23, the jet propagates to ∼300 km along the rotational axis at around ∼18 ms after bounce. In our counterpart model, the MHD-driven shock revives after around 30 ms after bounce, and it reaches 300 km at around ∼10 ms, which is equivalent to ∼40 ms after bounce. Considering that our model (β = 0.1%) is a slower rotator than model R4E1CF (β = 0.5%), the delay of the shock revival for our model seems reasonable. Model s20A1B5-D3M12 in Cerdá-Durán et al. (2008), whose precollapse angular velocity is 4 rad s⁻¹ (the rotational parameter should be close to β = 0.1%) with uniform rotation and whose initial poloidal magnetic field is 10¹² G, is close to our model B12X20β0.1. From their Figure 10, the MHD jet propagates to ∼240 km at 51 ms after bounce. The employed EOS is the same as ours (the Shen EOS), while the deleptonization scheme taken in their study was the Y_e formula (Liebendörfer 2005). As mentioned above, the dynamics are rather close to our corresponding model. Among the computed models in Burrows et al. (2007a), model M15B12DP2A1H, which has a precollapse angular velocity of π rad s⁻¹ (the rotational parameter should be close to β = 0.06%) with an initial dipolar magnetic field of 10¹² G, is close to our model B12X20β0.1. The interval before the launch of the MHD shock for their model is 80 ms after bounce (e.g., their Table 1), which is much later than it is for our model (28 ms after bounce). This may be due to the larger initial angular momentum (β = 0.1%) assumed in our study. Model A3B3G5-D3M13 in Obergaulinger et al. (2006b), which has a rotational parameter of β = 0.9% with a differential rotation imposed (the radial cutoff is 500 km) and the initial poloidal magnetic field is 10¹³ G, is closer to our model B12X1β0.1. The MHD jet reaches to 500 km at around 7 ms after bounce, which is also the case of our counterpart model. As discussed above, our results are compatible to the ones obtained in the relevant foregoing results.

The major limitation of this study is the assumption of axisymmetry. Recently, it was reported in three-dimensional (3D) MHD core-collapse simulations (Scheidegger et al. 2008, 2010) that the fast growth of the spiral SASI hinders the efficient amplification of the toroidal fields, which could suppress the formation of jets that are rather easily realized in two-dimensional simulations. As a sequel of this study, we plan to investigate the 3D effects in SRMHD. Regarding a resolution dependence of our results, Figure 14 indicates that our standard resolution is adequate to follow the evolution of the computed models. However, it is not at all sufficient to capture the magneto-rotational instability (MRI; e.g., Balbus & Hawley 1998). At least 10–100 times finer mesh points are required for resolving the MRI (Obergaulinger et al. 2009), which may require some adaptive mesh refinement treatment, a very important component that remains to be improved. If the MRI could be resolved, the increasing type waveform could also emerge for models with weak initial magnetic fields because a more efficient field amplification could be captured. Although the GR effects were treated only by rough approximation, we think that the GR should not drastically change our results qualitatively because the central protoneutron stars will not collapse to a black hole during our simulation time as inferred from a simple argument of the compactness of the inner core. As for the microphysics, the neutrino heating is not included in this study. However, the inclusion of the neutrino heating may play a minor role in the waveforms since the timescales before the neutrino-driven explosions set in (e.g., Marek & Janka 2009; Bruenn et al. 2009; Suwa et al. 2010, and references therein) are much longer than the MHD explosions.

Zoom In Zoom Out Reset image size

As one possible extension of this study, we think it interesting to study GWs from anisotropic neutrino radiation in the MHD models. Extrapolating the outcome obtained in previous studies (e.g., Müller et al. 2004; Kotake et al. 2009a), we anticipate that the total amplitudes become larger when we include the neutrino GWs. This is because the neutrino GWs may make a positive contribution since the neutrino emissions from the oblately deformed protoneutron stars become stronger toward the polar direction (Janka & Moenchmeyer 1989; Kotake et al. 2003a; Walder et al. 2005; Ott et al. 2008). If this is really the case, the total GW amplitudes especially for the increasing type should be much larger, possibly making its detection more promising. Furthermore, the neutrino signals from MHD explosions may have a sharp directional dependence through neutrino oscillations, reflecting the aspherical propagation of the shock to the stellar envelope (Kawagoe et al. 2009). Taking a correlation analysis between the neutrino and GW signals could help to reveal the hidden nature of the central engines. In fact, several observational proposals have been made in this direction recently (van Elewyck et al. 2009; Aso et al. 2008; Leonor et al. 2010). The MHD-driven core-collapse supernovae, albeit rather minor among typical type II supernova explosions, seem to still contain a number of rich research subjects.

We are grateful to K. Sato and S. Yamada for continuing encouragements. We also thank our anonymous referee for a number of suggestions that helped us to improve our manuscript. Numerical computations were in part carried out on XT4 and a general common use computer system at the center for Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan. This study was supported in part by Grants-in-Aid for Scientific Research from the Ministry of Education, Science, and Culture of Japan (nos. 19104006, 19540309, and 20740150).