GRAVITATIONAL WAVE SIGNATURES OF HYPERACCRETING COLLAPSAR DISKS

Regarding our precollapse model, we take a 35 M_☉ star (model 35OC) in Woosley & Heger (2006) that is supposed to be one of the most promising GRB progenitor models. To study the effects of rotation systematically, we take the following rotational profiles

$$\begin{equation} \Omega (r,\theta)=\frac{\Omega _0 X_0^4 + \alpha \Omega _\mathrm{lso}(M(X))X^4}{X_0^4+X^4}, \end{equation} \tag{ 1 }$$

where α, Ω₀, and X₀ are model parameters. M(X) is the mass coordinate at the cylindrical radius (X = rsin θ), and Ω_lso is given by Ω_lso = j_lso/X², where j_lso is the specific angular momentum in the last stable orbit of the Schwarzschild BH (e.g., Bardeen et al. 1972; Proga 2005; López-Cámara et al. 2009). By changing α in the range of 0.6 ⩽ α ⩽ 1.2, we compute seven non-magnetized models of α = 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, and 1.2, where they are labeled as models J0.6, J0.7, and so on (see Figure 1). To see the effects on the magnetic fields, we compute two more models: one in which the initial magnetic field (B₀ = 10¹⁰ G or 10¹¹ G) is added to model J0.8 (model J0.8B10 or model J0.8B11) and another in which the initial field is assumed to be purely poloidal and also assumed to be uniform and parallel to the spin axis.

Zoom In Zoom Out Reset image size

Figure 1 shows profiles of the initial angular momentum for some models. As shown, our models are taken to possess much more rapid rotation especially in the range of 1–4 M_☉ in the mass coordinate compared to the original profile of model 35OC (black dashed line). We set such a rotational profile, otherwise the accretion disk cannot survive later than ∼2 s after the onset of gravitational collapse in our simulation. Later on, the accretion disk is swallowed by the central object due to the neutrino cooling which removes the pressure support in the disk. In this case, one cannot account for the duration of long bursts, the activity of which is typically longer than ∼2 s, and can last up to tens of minutes. So, we experimented by choosing to adjust the initial angular momentum as in Equation (1).

Although our models do have higher initial angular momentum, the deviation may not be so serious (compare the light blue line that is for the angular momentum amplified by a factor of two in the original progenitor) considering uncertainties in stellar evolution calculations (such as in the treatment of mass loss, weak interactions, and fluid instabilities (e.g., convection, semiconvection, rotation, and magnetic fields)). To mimic the original progenitor structure, we take X₀ to be the size of the Fe core (≈3000 km). By setting α = 0.8, Ω₀ ≈ 1 rad s⁻¹, the above profile becomes closest to the original profile for mass coordinates larger than 5 M_☉. As a side note, the most prevalent way is to tune the initial angular momentum to satisfy its local centrifugal force to be several percent levels of the local gravitational binding energy (e.g., MacFadyen & Woosley 1999). Our modeling may be more realistic in the sense that the assumed initial angular momentum profiles capture the basic trend obtained in the current progenitor models.

2.1.1. Hydrodynamics

2.2.1. GWs from Matter Motions

To extract the gravitational waveforms from matter motions and magnetic fields, we employ the stress formulae derived in Takiwaki & Kotake (2011) and Kotake et al. (2004). For convenience, we briefly summarize them in the following.

In our axisymmetric case, the only non-vanishing quadrupole term is the plus mode (h₊) in the metric perturbation (e.g., Mönchmeyer et al. 1991), which can be written as

$$\begin{equation} h ({\mbox{\boldmath $X$}},t) = \frac{1}{8}\sqrt{\frac{15}{\pi }}\sin ^2 \theta \frac{A^{E2}_{20}\left(t - \frac{R}{c}\right)}{R} \end{equation} \tag{ 2 }$$

(e.g., Thorne 1980). The quadrupole amplitude of matter GWs A^E2_20matter consists of the following three parts

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

The first term in Equation (3) represents the contribution from non-spherical hydrodynamic motions, which is expressed as

$$\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 \rho _{*}W^2\big({v_r}^2(3 \mu ^2 -1)\nonumber \\ & & + {v_{\theta }}^2(2 - 3 \mu ^2) - {v_{\phi }}^{2} - 6 v_{r} v_{\theta } \,\mu \sqrt{1-\mu ^2}\big), \end{eqnarray} \tag{ 4 }$$

in which ρ_* is the effective density defined as

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

Here, ρ, e, p, c, and G denotes the baryon density, the internal energy, the pressure, the speed of light, and the gravitational constant, respectively, and |b|² = b^μb_μ is related to the energy density of the magnetic fields with b_μ representing the magnetic field in the laboratory frame (e.g., Takiwaki et al. 2009). $W= 1/\sqrt{1-v^kv_k}$ is the Lorentz boost factor with v_k denoting the spatial velocity in spherical coordinates (i = r, θ, ϕ). μ = cos θ is a directional cosine. The second term in Equation (3) represents the contribution from gravity as

$$\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\nonumber \\ & & \times \Bigg[\rho h (W^2+(v_k/c)^2)+ \frac{2}{c^2}\left(p+\frac{\left|b\right|^2}{2}\right)\nonumber\\ && -\frac{1}{c^2}((b^{0})^2+(b_{k})^2)\Bigg]\nonumber \\ & & \times [- r \partial _{r} \Phi (3 \mu ^2 -1) + 3 \partial _{\theta } \Phi \,\mu \sqrt{1-\mu ^2}], \end{eqnarray} \tag{ 6 }$$

where Φ denotes the gravitational potential of the self-gravity. The last term in Equation (3) is the contribution from the magnetic fields as

$$\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 \big[{b_r}^2(3 \mu ^2 -1)\nonumber \\ & & + {b_{\theta }}^2(2 - 3 \mu ^2) - {b_{\phi }}^{2} - 6 b_{r} b_{\theta } \mu \sqrt{1-\mu ^2}\big]. \end{eqnarray} \tag{ 7 }$$

Here, we write the total gravitational amplitude as follows for later convenience,

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

where the quantities of the right hand side of the equation are defined by combining Equations (2) and (3) with Equations (4), (6), and (7). Note that by dropping the O(v/c) terms, the above formulae reduce to the conventional quadrupole formula employed in Newtonian simulations (e.g., Mönchmeyer et al. 1991). In the collapsar disk that we consider in this work, the condition of v_ϕ ≫ v_θ, v_r in Equation (4) is generally satisfied inside the disk. If the disk is in a perfectly stationary state, which means that the centrifugal force (∼ρv²_ϕ) balances the gravitational forces (e.g., Equation (6)), no GWs can be emitted. As will be explained later, the disk attains mass continuously due to mass accretion, the specific angular momentum of which increases outward (e.g., Figure 1). This is the primary reason why non-zero matter GWs from axisymmetrically but dynamically rotating collapsar disks are generated. In the following computations, we assume that the observer is located in the equatorial plane (θ = π/2 in Equation (2)), and also that the distance to the GW source is comparable to nearby GRB-associated CCSNe (R = 100 Mpc) unless stated otherwise.

2.2.2. GWs from Anisotropic Neutrino Emission

To compute the gravitational waveforms from anisotropic neutrino radiation, we follow the formalism pioneeringly proposed by Epstein (1978) and Müler & Janka (1997). In the case of our 2D axisymmetric case, the only non-vanishing component is the plus mode for the equatorial observer,

$$\begin{eqnarray} h_{\nu } &=& \frac{4G}{c^4 R} \int _{0}^{t} dt^{{\prime }} \int _{0}^{\pi } d\theta ^{\prime } \Phi (\theta ^{\prime }) \frac{dl_{\nu }(\theta ^{\prime },t^{\prime })}{d\Omega ^{\prime }}, \end{eqnarray} \tag{ 9 }$$

where Φ(θ') depends on the angle measured from the symmetry axis (θ')

$$\begin{equation} \Phi ({\theta ^{{\prime }}})= \pi \sin \theta ^{{\prime }}(- 1 + 2 | \cos \theta ^{{\prime }}|). \end{equation} \tag{ 10 }$$

As given in Figure 1 of Kotake et al. (2007), this function has positive values in the north polar cap for 0° ⩽ θ' ⩽ 60° and in the south polar cap for 120° ⩽ θ' ⩽ 180°, but becomes negative between 60° < θ' < 120°. To determine the anisotropy in neutrino emission (i.e., dl_ν/dΩ in Equation (9)), we perform a ray-tracing analysis, which was initially proposed to be applicable in Newtonian gravity (Kotake et al. 2009a, 2009b, 2011) and later improved to be applicable to SR and GR (Harikae et al. 2010a, 2010b).

Applying the formalism in Lindquist (1966), the Boltzmann equation for the neutrino occupation probability f_ν(_ν, Ω) for a given neutrino energy (_ν) along a specified direction of Ω can be expressed as

$$\begin{equation} \frac{df_{\nu }(\epsilon _{\nu },\Omega)}{d\lambda } = n [Q_e(1-f_{\nu }) - \kappa f_{\nu }] = n[Q_e - \kappa ^{*} f_{\nu }], \end{equation} \tag{ 11 }$$

where $n(\mbox{\boldmath $x$})$ is the proper number density of the external medium with which neutrinos interact and thus measured in its own local rest frame, _ν is the neutrino energy measured in the local proper frame, λ denotes an affine parameter along the geodesics, Q_e and κ represents neutrino emissivity and absorptivity, and (1 − f) represents the Pauli blocking term (e.g., see Harikae et al. 2010b for more detailed information to derive the equation). Here, we consider only ν_e and $\bar{\nu }_e$ for simplicity, and the energy and momentum transfer via neutrino scattering are neglected, which is not only difficult to treat with the ray-tracing technique but is also a major undertaking in the radiative transport problem in general. In the final expression of Equation (11), κ* is defined to represent the effective absorptivity. As for the opacity sources of neutrinos, electron capture on proton and nuclei, positron capture on neutron, and neutrino scattering with nucleon and nuclei are included (Fuller et al. 1985; Takahashi et al. 1978; Bruenn et al. 2010). Here, κ* is estimated as κ* = Σ[n_target · σ(_ν)] with n_target and σ(_ν) being the target number density of each reaction and the corresponding cross-sections, respectively.

According to Zink (2008), the formal solution of Equation (11) can be given as

$$\begin{eqnarray} f_{\nu }(\epsilon, \Omega) &=& \int ^{\lambda _{\rm out}}_{\lambda _{0}} n(\lambda ^{\prime \prime }) Q_e(\lambda,^{\prime \prime } \epsilon _{\nu })\nonumber\\ &&\times \exp \Biggl [{{-}\int ^{\lambda _{\rm out}}_{\lambda ^{\prime \prime }} n(\lambda ^{\prime }) \kappa ^{*}(\lambda ^{\prime }) d\lambda ^{\prime } }\Biggr ] d\lambda ^{\prime \prime }, \end{eqnarray} \tag{ 12 }$$

which is referred to as the rendering equation of the radiation transport problem. We perform a line integral along the geodesics from every point on the surface of neutrinospheres (λ₀) from which neutrinos can escape freely, up to the outermost boundary of the computational domain (λ_out). On the neutrinospheres, the neutrino distribution function is assumed to take a Fermi–Dirac type distribution f^FD[ = 1/(_ν/k_BT + 1)] with a vanishing chemical potential.¹¹ For the neutrino energy bins (_ν), we use 16 logarithmically spaced energy bins reaching from 3 to 300 MeV.

The path integration in Equation (12) is done explicitly along the geodesics. In doing so, we determine each integration step by restricting the maximum change of neutrino opacity for all the neutrino energy bins to be less than 10%. By also using an adaptive mesh refinement approach, our ray-trace code was proven to safely pass several test problems (see Harikae et al. 2010b for more details). For the ray-tracing calculation, 300(r) × 40(θ) × 32(ϕ) mesh points are cast. The radiation field typically changes more slowly than the hydrodynamics. We perform the ray-tracing calculation every 50 hydrodynamic steps to reduce the computational cost. Note that the obtained waveforms did not change significantly when the interval is varied every 100 or 200 steps.

With f(_ν, Ω) at the outermost boundary, which is obtained using the above procedure, the emergent neutrino energy fluxes along the specified direction of Ω can be estimated as

$$\begin{equation} \frac{dl_{\nu }(\Omega,\epsilon _{\nu })}{d\Omega dS} = \int f(\epsilon _{\nu },\Omega)\cdot (c\epsilon _{\nu }) \cdot \frac{\epsilon _{\nu }^2 d\epsilon _{\nu }}{(2 \pi \hbar c)^3}. \end{equation} \tag{ 13 }$$

By summing up the energy fluxes with the weight of the area in the plane perpendicular to the rays (:dS), we can find dl_ν/dΩ along the specified direction Ω,

$$\begin{equation} \frac{dl_{\nu }({\Omega })}{d\Omega } = \int \frac{dl_{\nu }({\Omega })}{d\Omega dS} dS. \end{equation} \tag{ 14 }$$

Repeating the above procedures, dl_ν(Ω)/dΩ can be estimated for all directions and by which we can find the amplitudes of the GWs from neutrinos through Equation (9).

To capture hydrodynamic features in our models, Figure 4 shows the evolution of neutrino luminosities for some selected models. In t ∼ 1.0 s after we start our simulations (t ≡ 0 s), all the models experience rapid infall and the subsequent shock formation in the center, which leads to a drastic increase and the subsequent decrease in the neutrino luminosities.¹²

Zoom In Zoom Out Reset image size

Here, we take model J0.8 as a reference, because the precollapse angular momentum is adjusted to be closest to the original progenitor. The top left panel in Figure 5 shows a snapshot at t = 0.42 s near the shock formation. Note that the accretion mass in this epoch is typically greater than 2–3M_☉ in the center. The maximum mass of the NS of the Shen EOS is in the same mass range (Shen et al. 1998, see also O'Connor & Ott 2011; Kiuchi & Kotake 2008). Recent full GR simulations show that the mass of the central object just before collapsing to a BH is typically ≲ 2.3 M_☉ (e.g., Ott et al. 2011) with typical radius of several km (as inferred from their Figure 3). These evidences might support our very crude assumption of the prompt BH formation that is modeled by setting the BH initially in the center, although such assumption can only be tested by full GR simulations using the same progenitor model. Later on, the density configuration (compare the left half in each panel) deforms to become more oblate with time due to accretion of material with higher angular momentum outside (e.g., Figure 1). As will be explained in the later section, the luminous accretion disk is the primary source of GWs from anisotropic neutrino emission.

Zoom In Zoom Out Reset image size

The top right panel in Figure 5 shows a snapshot at t = 1.89 s. The bluish regions near the spin axis of the accretion disk (left half, density) correspond to the so-called polar funnel regions. The entropy becomes highest near the surface of the accretion disk due to the shock heating when the accreting material hits the wall of the disk (right panel). Comparing Figure 3 (left panel) to Figure 4, the GWs from neutrinos typically deviate from zero later than t ≳ 2–3 s, when the (total) neutrino luminosities become as high as ∼10⁵² erg s⁻¹. Later on, the neutrino luminosities show a gradual increase with time, reflecting the increase in the mass accretion to the newly formed accretion disk.

In Figure 4, the neutrino luminosity for models J0.6 (red line) and J0.7 (green line) is shown to steeply decrease at ∼3.3 s and 7.5 s, respectively. This is because the accretion disk is swallowed to the center (bottom panel in Figure 6) mainly because the pressure support in the accretion disk is reduced by the neutrino cooling. The disappearance of the accretion disks is also the reason for the sudden decrease in the GW signals observed both in the neutrino and matter sectors (e.g., Figures 2 and 3). Before the disappearance, the accretion disk is observed to show rapid expansion and contraction (as indicated in the top panels of Figure 6). At the same time, the mass flux on the Lagrange point around the disk changes violently, which may be related to the so-called runaway instability (e.g., Abramowicz et al. 1983; Font & Daigne 2002). Except for these slowing rotating models, the neutrino luminosities gradually settle to become nearly constant typically later than t ∼ 6 s (Figure 4). The saturation of the neutrino luminosity is because the disk is already deleptonized and the neutrino emission there is suppressed (e.g., Harikae et al. 2009 for more details). The neutrino luminosities at this epoch become as high as 10⁵²–10⁵³ erg s⁻¹, which touches the level of ∼10% of the accretion luminosities (see also Chen & Beloborodov 2007; Sekiguchi & Shibata 2011). Apparently the accretion disk is neutrino-cooling dominated.

Zoom In Zoom Out Reset image size

From Figure 4, the neutrino luminosity of model J0.8 is shown to be the highest (e.g., blue line). Models with higher initial angular momentum have more extended disks due to larger centrifugal forces (compare the bottom panels in Figures 5 and 7), leading to lower density and temperature in the disks. This is the reason why the neutrino luminosities for models J1.0 (light blue line) and J1.2 (yellow line in Figure 4) are lower in this order. Reflecting this, the GWs from neutrinos become highest for model J0.8 (left panel of Figure 3; see also |h_ν| in Table 1). Regarding the matter GWs, it should be noted that their maximum absolute amplitudes are obtained not for the most rapidly rotating model (model J1.2), but for model J0.8 (right panel of Figure 3). This situation may be akin to the GW signals emitted near core-bounce in CCSNe (see Kotake et al. 2006, 2011; Ott 2009 for recent reviews). Too much initial angular momentum works to suppress the matter compression leading to the smaller mass-quadrupole moment. As a result, the matter GWs become maximum for models with moderately rotating models (e.g., Mönchmeyer et al. 1991; Yamada & Sato 1995; Zwerger & Müller 1997; Kotake et al. 2003b; Shibata & Sekiguchi 2004; Ott et al. 2004, 2007; Dimmelmeier et al. 2002, 2007, 2008; Scheidegger et al. 2010).

Zoom In Zoom Out Reset image size

In this section, we move on to look more into the detail of the properties of gravitational waveforms mentioned in the previous section. By taking model J0.8 as a reference, we first focus on the neutrino GWs.

Figure 8 shows the local neutrino energy fluxes (dl_ν/(dΩdS), Equation (13)) at t = 9.0 s (e.g., bottom panel of Figure 5) seen from the polar (left) and the equatorial direction (right). Note that the polar direction is taken to be parallel to the spin axis of the accretion disk. The top panel is for the Minkowski geometry, which corresponds to a purely Newtonian case. Reflecting the shape of the accretion disk in axisymmetry, the contours of the neutrino flux are deformed oblately when seen from the equator (top right), while the neutrino flux seen from the polar direction (top left) looks like a superposition of circles in a concentric fashion, the center of which corresponds to the spin axis. Note that the small central cavity in the top left panel corresponds to the inner boundary of our SR simulation.

Zoom In Zoom Out Reset image size

The middle panels are calculated for the Minkowski geometry including SR corrections (see Harikae et al. 2010b for more detail). Comparing the middle right to the top right panel, a clear difference is the left–right asymmetry in the middle right panel. This is due to the SR beaming effects. Since the rotational velocity of the accretion disk is perpendicular to the polar direction, the SR beaming effect suppresses the neutrino emission toward the polar region. As a result, the neutrino fluxes become dark (right-hand side in the middle right panel) because the direction of the rotating material is opposite the direction to the observer. To study the GR effects on the neutrino luminosity, we place 4 M_☉ in the central region, which mimics the event horizon of the BH. The effects of the BH spin are examined by setting the Kerr parameter by hand to be a = 0 and a = 0.999 for the Schwarzschild and the extreme Kerr geometry, respectively. The bottom panels of Figure 8 are for the extreme Kerr geometry, which, however, looks very similar to the middle panels.

For a more detailed comparison, the left panel of Figure 9 shows the neutrino luminosity per solid angle dL_ν/dΩ for the Minkowski (indicated by "SR"), Schwarzschild ("GR(a = 0))," and extreme Kerr geometry ("GR(a = 0.99))." The most important message in this panel is that in every case, the neutrino luminosity seen from the direction near the parallel to the spin axis (θ = 0) becomes higher than the one seen from the equatorial direction (θ = π/2). This is because the cross-section of the pancake-like accretion disk seen from the spin axis becomes larger compared to the one seen from the horizontal direction (Figure 8). Remembering again that Φ(θ') in Equation (10) is positive near the north and south polar caps, the dominance of the polar neutrino luminosities makes the neutrino GWs positive in the polar cap regions (right panel in Figure 9). This is the reason why the neutrino GWs increase monotonically with time (e.g., left panel of Figure 3). Here, it is worth mentioning that from neutrino luminosity (L_ν in Figure 4), their typical duration (Δt_ν), anisotropy in the neutrino radiation (α_ν, read from the left panel of Figure 9), the GW amplitudes from neutrinos can be estimated from Equation (9) to be

$$\begin{equation} h_{\nu } \approx 10^{-22} \Bigl (\frac{\alpha _{\nu }}{0.2}\Bigr)\Bigl (\frac{L_{\nu }}{10^{53} {\rm \,erg} {\rm \,s}^{-1}}\Bigr) \Bigl (\frac{\Delta t_{\nu }}{10 {\rm \,s}}\Bigr) \Bigl (\frac{D}{100 {\rm \,Mpc}}\Bigr)^{-1}, \end{equation} \tag{ 15 }$$

which is in good agreement with the numerical results (Figure 2).

Zoom In Zoom Out Reset image size

Comparing the SR case with the GR case in the left panel of Figure 9, the neutrino luminosity becomes smaller by ∼20%–30% in the GR case due to the GR redshift and also due to the bending effects. Comparing the green line with the blue line, the frame-dragging effect due to the BH spin barely affects the emergent neutrino luminosity. This is probably because the GR effect on the neutrino luminosity becomes important only in the central regions very close to the BH (e.g., Figure 15, Harikae et al. 2010b). Our results indicate that a ray-tracing calculation in the Schwarzschild geometry is at least needed to accurately estimate the neutrino GWs in the collapsar's environment.

Now we briefly analyze the matter GWs in model J0.8. As summarized in Section 2.2, the matter GWs are estimated by the sum of the hydrodynamic and gravity parts in Equations (4) and (6) for non-magnetized models. Among the kinetic energy terms (∝ρv_iv_j) in Equation (4), the largest contribution comes from the rotational energy (i.e., −ρ_*W²v²_ϕ), which is shown in Figure 10 (green line). In contrast to the negative contribution by this term, the gravity part (blue line labeled "Gravity" in Figure 10) is shown to make a positive contribution. If a star rotates perfectly stationary, the centrifugal forces balance out the gravitational force, leading to no GWs. In the collapsar disk studied in this work, the disk attains mass continuously due to mass accretion with its specific angular momentum increasing outward. This is the reason why the disk is not perfectly stationary, leading to a non-zero GW emission from matter motions (see pink line labeled by "Matter" in Figure 10). However, their GW amplitudes are much smaller compared to those from anisotropic neutrino emission (compare the pink line with the light blue line (neutrino GWs) in Figure 10).

Zoom In Zoom Out Reset image size

The left panel of Figure 11 shows the normalized contribution of each term in A^E2₂₀ at t = 9 s for model J0.8, which is estimated by the volume integral of A^E2₂₀ within a given sphere enclosed by a certain radius. The region within a radius between 50 km and 160 km is shown to contribute to the production of GWs, which corresponds to a high-density region in the accretion disk (see the bottom panel in Figure 5). It can be also shown that the radial gradient of the gravitational potential (green line indicated "Gravity (r)" in the plot) closely cancels with the contribution from the centrifugal force (red line) in the disk. As shown in the right panel of Figure 11, the enclosed mass there becomes as big as ∼4 M_☉ as a result of the hyperaccreting activity lasting ∼ 9 s until that, and the typical rotational period there is time on the order of milliseconds (red line).

Zoom In Zoom Out Reset image size

Putting these numbers to the standard GW stress formula (e.g., Shapiro & Teukolsky 1983), an upper bound of the matter GW amplitudes may be estimated to be

$$\begin{eqnarray} h_{\rm matter} &=& \frac{2G}{c^4 D} \ddot{I_{ij}} \sim \frac{2G}{c^4 D} \frac{{\it MR}^2}{T^2}\lesssim 10^{-23} \Bigl (\frac{\epsilon }{0.1}\Bigr) \Bigl (\frac{100 {\rm \,Mpc}}{D}\Bigr)\nonumber\\ &&\times\Bigl (\frac{M}{4 \,M_{\odot }}\Bigr) \Bigl (\frac{R}{100 {\rm \,km}}\Bigr)^2 \Bigl (\frac{T}{4 {\rm \,ms}}\Bigr)^{-2}, \end{eqnarray} \tag{ 16 }$$

where D is the distance to the source, $\ddot{I_{ij}}$ is the second time derivative of the quadrupole moment of I_ij, and M, R, and T represents the typical mass and radius of the accretion disk and the timescale, which we may take to be the rotational period, respectively, and is the degree of the non-sphericity, which we take optimistically to be 10% in the case of the accretion disk. This estimate, roughly consistent with the numerical results (e.g., right panel of Figure 3), also shows that anisotropic neutrino emission, producing the GW amplitude on the order of ≈10⁻²² for a source of 100 Mpc, is the primary source in the long-term evolution of collapsars (∼10 s; see Equation (15)).

We now discuss the effects of magnetic fields on the waveforms by taking model J0.8B11, which has the strongest initial magnetic field in our models. In Figure 12, the pink line represents the GWs from magnetic fields (Equation (7)). But, first of all, let us briefly comment on the burst-like feature of the green line (contributed from matter motions) at t ∼ 0.24 s in the figure. This may look similar to the burst of GWs emitted at the moment of BH formation possibly followed by a damped sinusoidal oscillation (Seidel 1990, 1991), but this can only be captured in full GR simulations (Baiotti et al. 2007; Ott et al. 2011; Sekiguchi & Shibata 2011; Kuroda et al. 2012). As already seen in the luminosity plot of Figure 4, the above burst simply corresponds to the shock formation at the center in our SR models. In our main point, we analyze the increasing trend of the GWs from magnetic fields in the following.

Zoom In Zoom Out Reset image size

The left panel of Figure 13 depicts a snapshot of density (left half), entropy (right top), and plasma β (right bottom, the ratio of magnetic to the matter pressure) at the final simulation time (t = 541 ms) for model J0.8B11. The high entropy regions in the slightly off-axis region (seen as reddish in the top right panel) correspond to the magnetohydrodynamically driven outflows that are pushed outward by the twisted toroidal magnetic fields (also seen as reddish in the bottom right panel, indicated "Toroidal").

Zoom In Zoom Out Reset image size

The right panel of Figure 13 shows contributions to the total GW amplitudes (Equation (9), in which the left-hand-side panels are for the sum of the hydrodynamic and gravitational parts (indicated "Matter"), namely, log (± [A^E2_20(hyd) + A^E2_20(grav)]) (left top(+)/bottom(−) (Equations (4) and (6)), and the right-hand-side panels are for the magnetic part, namely, log  [(± A^E2_20(mag))] (right top(+)/bottom(−)) (e.g., Equation (7)). By comparing the top two panels, it can be seen that the positive contribution overlaps with 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 (4)). The magnetic part also contributes to the positive trend (see top right half in the right panel (labeled by mag(+))). This comes from the toroidal magnetic fields (e.g., +b_ϕ² in Equation (7)), which dominantly contribute to drive MHD explosions.

Unfortunately, our MHD code becomes numerically unstable when the strong MHD jets propagate to a stellar mantle with decreasing density, which prevents us from studying the resulting GWs for a longer period. The neutrino GWs are much smaller than the other GW sources (e.g., Figure 12)) simply due to the shorter simulation time. We expect that the increasing trends of the magnetic fields remain as the MHD shocks propagate farther out (e.g., Takiwaki & Kotake 2011). To confirm it, we need to implement a numerical technique especially developed to solve the force-free fields, which is a major undertaking (e.g., McKinney 2006).

Finally, Figure 14 depicts the GW spectra for models J0.8 (left) and J1.0 (right). The GW spectra in the frequency domain between 1 and 100 Hz becomes slightly larger for model J0.8 than for model J1.0. This reflects a more efficient release of the gravitational binding energy for the moderately rotating case (model J0.8) as already mentioned in Sections 3.2 and 3.4. For the cosmological distance scale of GRBs (∼100 Mpc), these low-frequency GW signals are unfortunately very hard to detect even by advanced detectors (like the advanced LIGO or KAGRA/LCGT) whose sensitivity is severely limited by seismic noises (Ando & the TAMA Collaboration 2005; Abbott et al. 2005; Weinstein 2002; Kuroda & LCGT Collaboration 2010). The good news is that these signals could be detectable by a recently proposed future space-based interferometer-like Fabry–Perot-type DECIGO (Kawamura et al. 2006; black line in Figure 14). Two low-luminosity LGRBs were already observed at distances of ∼40 Mpc (980425, Galama et al. 1998) and ∼130 Mpc (060218, Ferrero et al. 2006). Their local rate, being much higher than that of normal bursts, is expected to be as large ≈0.1D₁₀₀ yr⁻¹ with D₁₀₀ representing the distance normalized by 100 Mpc (see discussions in Corsi & Mészáros 2009). Our results suggest that the GW astronomy of collapsars could become a reality with the DECIGO-class GW detectors, hopefully in the near future.

Zoom In Zoom Out Reset image size