A&A 490, 231-241 (2008)
DOI: 10.1051/0004-6361:20078577

Gravitational waves from 3D MHD core collapse simulations

S. Scheidegger - T. Fischer - S. C. Whitehouse - M. Liebendörfer

Department of Physics, University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland

Received 30 August 2007 / Accepted 2 July 2008

Abstract
We present the gravitational wave analyses from rotating (model s15g) and nearly non-rotating (model s15h) 3D MHD core collapse supernova simulations at bounce and during the first couple of ten milliseconds afterwards. The simulations are launched from 15 $M_{\odot}$ progenitor models stemming from stellar-evolution calculations. Gravity is implemented by a spherically symmetric effective general relativistic potential. The input physics uses the Lattimer-Swesty equation of state for hot, dense matter and a neutrino parametrisation scheme that is accurate until the first few ms after bounce. The 3D simulations allow us to study features already known from 2D simulations, as well as nonaxisymmetric effects. In agreement with recent results, we find only type I gravitational wave signals at core bounce. In the later stage of the simulations, one of our models (s15g) shows nonaxisymmetric gravitational wave emission caused by a low T/|W| dynamical instability, while the other model radiates gravitational waves due to a convective instability in the protoneutron star. The total energy released in gravitational waves within the considered time intervals is $1.52\times10^{-7}~M_{\odot}$ (s15g) and $4.72\times10^{-10}~M_{\odot}$ (s15h). Both core collapse simulations indicate that corresponding events in our Galaxy would be detectable either by the LIGO or Advanced LIGO detector.

Key words: gravitational waves - supernovae: general - hydrodynamics - neutrinos - stars: rotation - stars: neutron

1 Introduction

Gravitational wave astronomy could well make its first observations soon as running earth-based facilities such as LIGO[*], VIRGO[*], GEO600[*], TAMA[*], and AIGO[*] are able to reach the required high sensitivities. For a review of gravitational wave detection see e.g. (Hough & Rowan 2007). The successful observation of gravitational waves (GW) would be a major breakthrough by opening a new window on the universe, allowing us to observe electromagnetically hidden regions (Thorne 1995). One of the most promising gravitational wave sources is the stellar core collapse leading to a supernova explosion. The observation of both the GW and the neutrino signal from a galactic supernova would reveal hitherto hidden details about the explosion scenario and impose constraints on the nuclear and weak interaction physics under conditions that cannot be obtained in terrestrial experiments.

As a core of a star of $ M \geq 8~M_{\odot} $ reaches the end of its stellar evolution, it becomes gravitationally unstable as soon as thermonuclear burning produces a significant amount of iron-group nuclei. Next to photo-disintegration, electron captures on free protons and nuclei reduce the mostly electron-supported pressure, and the core eventually starts to collapse. The collapse continues until nuclear densities of $\sim$ $2\times10^{14}$ g/cm3 are reached, depending on the equation of state (EoS). As soon as we enter this density regime, the core overshoots its equilibrium position and bounces back. A sound wave immediately forms and steepens into a shock wave that propagates outwards. However, the shock stalls within $\sim$5 ms after core bounce due to the large energy loss caused by the dissociation of nuclei into free nucleons at a cost of $\sim$8.8 MeV per nucleon and by neutrino emission connected to copious electron captures on the emerging free protons. It continues to propagate outwards to radii around 100-200 km as standing accretion shock. Thereafter some ms after the prompt explosion mechanism has failed, a delayed explosion mechanism by neutrino heating is thought to occur, e.g. Bethe (1990).

The idea of reviving the stalled shock again via neutrino reactions behind and ahead of the shock has long been investigated as a possible explosion mechanism (Janka et al. 2001; Bethe & Wilson 1985). In the dissociated matter behind the shock, electron-flavour neutrino captures are the dominant heating reactions, while in the accreting matter ahead of the shock, electron-neutrino absorption and neutrino-nucleon scattering are the dominant interactions that preheat the infalling unshocked material (Bruenn & Haxton 1991). In addition, neutrino-electron scattering and pair annihilation may also contribute to the heating. However, the effects of preheating and neutrino-electron scattering on the shock revival are insignificantly weak compared to the effect from the capture of neutrinos on free nucleons in the heating region below the shock. Equally as important as the neutrino heating is the neutrino cooling of matter that settles on the protoneutron star (PNS) (Janka 2001).

As a core collapse supernova is also likely to show aspherical features close to its centre (Leonard et al. 2006), where the matter assumes a high density, a tiny fraction of the released binding energy can be emitted via gravitational radiation. The following features have been suggested as possible causes of the asymmetries in the energy-matter distribution needed to emit GW: the rotational stellar core collapse, convection in the high-density protoneutron star, nonaxisymmetric rotational instabilities, fluid instabilities in the lower-density hot mantle surrounding it (possibly triggering neutron star oscillation modes), and an anisotropic neutrino emission (for recent reviews see e.g. Fryer & New 2006; Kotake et al. 2006).

The understanding of GW emission from differentially rotating core collapse has evolved with time due to the improving input physics used in simulations. In 2D axisymmetric computations by Müller (1982), iron cores, a Newtonian hydrodynamics code, and a tabulated finite temperature EoS were used. The results allowed the link between rotation and the efficiency of GW emission to be recognised. Having taken more micro-physics into account in Newtonian gravity, such as electron capture on protons and a simplified neutrino transport scheme (Moenchmeyer et al. 1991), it was possible to distinguish two different wave-form characteristics, later known as types I and II. Zwerger & Müller (1997) performed an extensive parameter study of a wide variety of models, using progenitors in rotational equilibrium, Newtonian gravity, and a polytropic EoS, but neglecting electron capture and neutrino physics. The results led to the introduction of types I, II, and III wave forms and their quantitative distinction in dependence on the stiffness of the EoS and rotation. A type I waveform is characterised by a large amplitude peak at core bounce and subsequent damping ring-down oscillations. It appears as if the influence of the initial angular momentum is small. In this case the core bounce occurs by the stiffening of the equation of state and is pressure-dominated. A type II signal occurs if the core bounce around nuclear density is dominated by strong centrifugal forces, and it has several distinct peaks caused by multiple `centrifugal' core bounces with following coherent re-expansion phases of the inner core. The type III signature is defined by a high positive peak at bounce followed by some smaller oscillations with very short periods. It appears when there is fast, pressure-dominated core collapse due to very efficient electron capture. Continued improvement in the input physics (i) by using realistic equations of state or neutrino physics; (ii) by adding magnetic fields or incorporating GR effects; and (iii) by extending the dimensionality from 2D to 3D, permitted us to deepen and underline the qualitative and quantitative understanding of GW forms from rotational core collapse around bounce time (Müller et al. 2004; Dimmelmeier et al. 2007; Janka & Müller 1996; Kotake et al. 2003; Dimmelmeier et al. 2005; Kotake et al. 2004; Fryer et al. 2002; Rampp et al. 1998; Ott et al. 2007b; Dimmelmeier et al. 2002). The most recent simulations by Ott et al. (2007a) and Dimmelmeier et al. (2007) point out that realistic input physics (GR, a micro-physical EoS, approximative description of deleptonisation) might lead solely to type I wave forms.

The investigation of GW emission resulting from convectively driven small-scale aspherities inside the PNS, the ``hot-bubble'' and anisotropic neutrino emission are still hampered by the requirement of computationally expensive neutrino transport in the postbounce evolution. Nevertheless, simulations that include accurate neutrino transport in axisymmetric models were explored with respect to GW in Müller et al. (2004) based on state-of-the-art progenitor models, GR-corrections, and sophisticated equations of state. The postbounce phase is believed to provide an interestingly large signal for hundreds of milliseconds, or even longer than 1 s, with a frequency distribution between $\approx $ 300-1200 Hz in case of convection, while the dominant frequencies lasting from anisotropic neutrino emission lie at about $\approx $10 Hz. Nonaxisymmetric dynamics has recently been investigated with neutrino physics approximations (Ott et al. 2007a). In these computations it was possible to show without the addition of any seed perturbation that differentially rotating protoneutron stars are likely to become dynamically unstable at low $\beta=T/\vert W\vert$-values (ratio of rotational to gravitational energy), leading to strong narrow-band GW emission (Saijo & Yoshida 2006; Ou & Tohline 2006; Ott et al. 2005; Saijo et al. 2003; Watts et al. 2005).

Another mechanism for gravitational wave emission in core collapse supernovae was proposed in Ott et al. (2006). It was pointed out that PNS core g-mode oscillations might be a major source of gravitational radiation. The release of GW in their simulations is highly correlated with the fundamental core g-mode, leading to a spectrum that peaks at twice the frequency of the aforementioned g-mode.

In this article we present the gravitational wave analysis from two 3D magneto-hydrodynamics (MHD) core collapse simulations. The simulations are based on a nuclear equation of state. The deleptonisation by electron capture and neutrino emission are accurately parametrised up to about 5 ms after bounce. Important general relativistic corrections to the gravitational potential are taken into account as well. The simulations are launched from progenitor models at the end of stellar evolution calculations (Woosley & Weaver 1995). One of the initial models (s15g) rotates moderately fast at the onset of the collapse with $2\pi$ rad/s for the central angular frequency, while the second model (s15h) rotates more slowly as estimated in Heger et al. (2005). We performed contiguous simulations from collapse into the postbounce phase and study three-dimensional effects on the emission of gravitational waves. The GW-producing asymmetries in three-dimensional simulations could be significantly different from the asymmetries present in axisymmetric models. As only two other groups have been performing three-dimensional supernova models with a nuclear equation of state and neutrino physics approximations for the prediction of gravitational wave signals (Dimmelmeier et al. 2007; Ott et al. 2007b; Fryer et al. 2004), our simulations provide an independent confirmation of these results. Rather than augmenting the database with numerous GW pattern templates, we aim to underline the importance of the neutrino physics in supernova GW signal predictions and to encourage the search for GW signals from Galactic core collapse supernova events.

This paper is organised as follows. In Sect. 2 we define the numerical methods and input physics used in our supernova models. We also describe the gravitational wave extraction formalism applied for the analysis of our three-dimensional data sets. Section 3 is devoted to the resulting gravitational waves, their properties, and the possible detectability of signals. In Sect. 4 we discuss the uncertainties in the numerical models and try to separate the robust features from details that are subject to change in future improved models. A conclusion is given in Sect 5.

2 Method

2.1 Description of the hydrodynamical models s15g and s15h

Only a few recent multi-dimensional collapse simulations have made the effort to include neutrino physics. These schemes are either very computationally expensive (Dessart et al. 2006; Buras et al. 2003) or rely on simplifications of the neutrino transport and its micro-physics (Fryer & Warren 2004a; Kotake et al. 2003). The complete Boltzmann neutrino transport equation can only be solved in spherical symmetry (Mezzacappa 2005, and references therein). Hence, three-dimensional hydrodynamical models must rely on neutrino transport approximations.

A simple and computationally efficient parametrisation of the deleptonisation in the collapse phase can be derived from a tabulation or fit of the electron fraction $ Y_{\rm e} $ as a function of density. Spherically symmetric models with Boltzmann neutrino transport (Liebendörfer et al. 2005) show that the $ Y_{\rm e} $ in different layers in the homologously collapsing core follow a similar deleptonisation trajectory, i.e. reach similar $ Y_{\rm e} $ values as a function of density. As changes in $ Y_{\rm e} $ can only be due to electron captures, it is possible to deduce corresponding changes in entropy and to calculate the neutrino stress from the emitted neutrinos (Liebendörfer 2005).

  \begin{figure} \par\includegraphics[width=8.8cm]{8577fig1.eps} \end{figure} Figure 1: Comparison of the almost non-rotating 3D model s15h with the spherically symmetric model G15, which is based on general relativistic three-flavour Boltzmann neutrino transport (Liebendörfer et al. 2005). From the upper left to the lower right we compare as a function of enclosed mass: the density-, the velocity, the $ Y_{\rm e} $, and the entropy profiles. The solid lines show the results of model s15h and the dotted lines the results of model G15. The thin lines represent a time instance at 5 ms before bounce and the thick lines represent a time instance at 5 ms after bounce. Excellent agreement is found in all four quantities - with one exception. The parametrised neutrino leakage scheme cannot model the neutrino burst. The neutrino burst causes a prominent $ Y_{\rm e} $-dip and additional cooling in the G15 data.
Open with DEXTER

Figure 1 shows a comparison of the 3D parametrised run s15h, launched from an almost non-rotational progenitor model, with a spherically symmetric general relativistic model G15 that is based on general relativistic Boltzmann neutrino transport (Liebendörfer et al. 2005). The 1D results of core collapse are accurately reproduced by the 3D run. The parametrised neutrino physics presents a significant improvement with respect to adiabatic simulations and may even rival neutrino transport schemes that neglect neutrino-electron scattering. However, the accuracy breaks down with the launch of the neutrino burst at a few milliseconds after bounce. With the currently implemented scheme, accretion flows in the postbounce phase only deleptonise down to $ Y_{\rm e}\sim 0.3 $instead of $ Y_{\rm e}\sim 0.15 $. Neutrino heating is neglected altogether. This turns the quantitative model of the collapse phase into a qualitative model of the postbounce phase.

The hydrodynamics of the 3D simulations is based on a simple and fast cosmological MHD code (Pen et al. 2003) that has been parallelised, improved, and adapted to the supernova context. A realistic equation of state (Lattimer & Swesty 1991) is used and the monopole term of the gravitational potential is implemented by a spherically symmetric mass integration that includes general relativistic corrections (Marek et al. 2006). The 3D computational domain consists of a central cube of 6003km3 volume, treated in equidistant Cartesian coordinates with a resolution of 1 km. The 3D code has buffer zones that require a prescription of the conditions in the immediate neighborhood of the 3D computational domain at each time step. To obtain these conditions, we embed the 3D computational domain in a larger spherically symmetrical computational domain that is treated by a spherically symmetric hydrodynamics code ``Agile'' (Liebendörfer et al. 2002). After each time step, the buffer zones at the boundary of the 3D hydrodynamics code are filled with the current conditions of the spherically symmetric solution. As long as the infall velocity at the boundary is subsonic, we spherically average the conditions in the 3D code and feed them back to the mass shells of Agile that are enclosed by the 3D domain. Once the infall velocities become supersonic, this step can be omitted because no hydrodynamic signal can leave the 3D domain anymore before the shock passes the boundary after the onset of an explosion.

For both of our models we use a 15  $M_{\odot}$ progenitor star of Woosley & Weaver (1995) as initial model. As the rotation rates of inner stellar cores are not very well-known, we assign angular momentum according to a simple parametrisation with a shellular quadratic cutoff at 500 km radius. The angular momentum is assumed to be conserved until the infalling layers enter the 3D computational domain. For one model, which we name s15h, we set an initial angular velocity of $ \Omega =0.3 $ rad/s as obtained in a stellar evolution model that takes account of the effects of rotation and magnetic fields (Heger et al. 2005). Because this value represents a rather slow rotation, we performed an alternative model, s15g, with an initial angular velocity of $ \Omega =2\pi $ rad/s. This is still a slow rotation rate compared to values assumed in many parametrised studies of the GW signal, but it allows for a clear distinction from model s15h. The initial $ \beta _{i}=T/\vert W\vert $ shows $ 0.26\% $ (s15g) and $ 0.00059\% $ (s15h), respectively. Both models collapse in a similar manner to spherically symmetric simulations. After bounce, the declining strength of the bounce-shock results in a negative entropy gradient and the protoneutron star becomes convectively unstable on a short time scale. The simulation s15g was carried out until about 70 ms postbounce, whereas s15h run $\sim$100 ms postbounce.

We have set the magnetic fields to the values suggested in Heger et al. (2005). The advection of the magnetic field into the boundary of the 3D computational domain is technically very delicate because the field lines assigned to the buffer zones have to ensure a vanishing divergence at the transition to the slightly evolved magnetic field just inside the boundary. We found a quite pragmatic solution where the field is first assigned according to above analytical setup of a divergence-free initial field. In a second step, the divergence at the boundary is analysed and corrected such that the transition is divergence-free. Due to its weakness, the magnetic field does not noticeably influence the hydrodynamical evolution during the early postbounce evolution.

2.2 Extracting gravitational radiation in three dimensions

Throughout this article we work in cgs-units and use the following values for the speed of light, $c=2.997\times 10^{10}$ cm s-1, the gravitational constant, $G=6.672\times10^{-8}$ cm3 g-1 s-2, and the parsec, 1 pc $=3.086\times 10^{18}$ cm. We do not assume any symmetry. In this case, the gravitational wave field hijTT can be decomposed into two orthogonal polarisations with amplitudes A+ and $A{\times}$, see e.g. (Misner et al. 1973; Finn & Evans 1990):

\begin{displaymath}h_{ij}^{TT}(\textbf{X},t)=\frac{1}{R}(A_{+}e_{+}+A_{\times}e_{\times}), \end{displaymath} (1)

where R is the distance to the source. In spherical coordinates, the unit polarisation tensors are given by

\begin{displaymath}e_{+}=e_{\theta}\otimes e_{\theta}-e_{\phi}\otimes e_{\phi} \end{displaymath} (2)

and

\begin{displaymath}e_{\times}=e_{\theta}\otimes e_{\phi}+e_{\phi}\otimes e_{\theta}. \end{displaymath} (3)

In the slow-motion limit (Misner et al. 1973; Finn & Evans 1990) the amplitudes A+ and $A_{\times }$ are given by linear combinations of the second time derivative of the transverse traceless mass quadrupole tensor:
$\displaystyle A_{+} = \ddot{t}_{\theta\theta}-\ddot{t}_{\phi\phi}$     (4)
$\displaystyle A_{\times} = 2\ddot{t}_ {\theta\phi}.$     (5)

In the Cartesian orthonormal basis, the quadrupole tensor is given by

 \begin{displaymath}t_{ij}^{TT}=\frac{G}{c^4}\int {\rm d}V \rho\left(x_{i}x_{j}-\frac{1}{3}\delta_{ij}r^2\right). \end{displaymath} (6)

Note that its time derivatives are evaluated at a retarded time. However, shortcomings exist that are related to numerical high-frequency noise and the r2 momentum-arm, which make the performance of a direct evaluation of Eq. (6) poor, as discussed in Finn & Evans (1990). Therefore, we use an alternative post-Newtonian expression from Blanchet et al. (1990), where the second-order time derivatives of the quadrupole moment are transformed into hydrodynamical variables that are known from the core collapse simulation:

 \begin{displaymath}\ddot{t}_{ij}^{TT}=\frac{G}{c^4}\int {\rm d}V \rho\left(2v_{i}v_{j}-x_{i}\partial_{j}\Phi-x_{j}\partial_{i}\Phi\right), \end{displaymath} (7)

where $\Phi$ is the gravitational potential. This expression allows the evaluation of the wave field from data at only one time instance; moreover, the integral has a compact support. This formula is not however gauge invariant and only valid in the Newtonian slow-motion limit. Nevertheless, it was shown that this approximation is accurate enough when compared to other approaches (Shibata & Sekiguchi 2003).

The polarisation modes can be obtained explicitly from a coordinate transformation, for example for the $\theta\theta$-component:

\begin{displaymath}t_{\theta\theta}=t_{ij}^{TT}\frac{\partial x^{i}}{\partial\theta}\frac{\partial x^{j}}{\partial\theta}\cdot \end{displaymath} (8)

This leads to the following non-vanishing components (Oohara et al. 1997):

\begin{eqnarray*}t_{\theta\theta}& = & (t_{xx}^{TT}\cos^{2}\phi+t_{yy}^{TT}\sin^... ...hi)+t_{xz}^{TT}\sin\theta\sin\phi-t_{yz}^{TT}\sin\theta\cos\phi. \end{eqnarray*}


For simplicity, we evaluate below only the gravitational wave amplitudes for $\theta=\phi=0$ (denoted below with the subscript I)
$\displaystyle A_{+\rm I} = \ddot{t_{xx}}-\ddot{t_{yy}}$     (9)
$\displaystyle A_{\times \rm I} = 2\ddot{t_{xy}}$     (10)

and for $\theta=\frac{\pi}{2}$, $\phi=0$ (denoted as II)
$\displaystyle A_{+\rm II} = \ddot{t_{zz}}-\ddot{t_{yy}}$     (11)
$\displaystyle A_{\times\rm II} = -2\ddot{t_{yz}}.$     (12)

We notice that $A_{\rm +II}$ corresponds to the non-vanishing quadrupole amplitude $A_{20}^{\rm E2}$ of axisymmetric models in the multipole expansion of the radiation field (Thorne 1980).

The energy carried away by gravitational radiation can be calculated by the following expression:

 
                  $\displaystyle E_{\rm GW}$ = $\displaystyle \frac{c^3}{5G}\int \left[\frac{{\rm d}}{{\rm d}t}(I_{ij}-\frac{1}{3}\delta_{ij}I_{ll})\right]^2 {\rm d}t$ (13)
  = $\displaystyle \frac{2c^3}{15G}\int {\rm d}t[\dot{I_{xx}}^2+\dot{I_{yy}}^2 + \dot{I_{zz}}^2$  
    $\displaystyle - \dot{I}_{xx}\dot{I}_{yy}- \dot{I}_{xx}\dot{I}_{zz}-\dot{I}_{yy}\dot{I}_{zz}$  
    $\displaystyle +3(\dot{I}_{xy}^2+\dot{I}_{xz}^2+\dot{I}_{yz}^2)],$  

where $I_{ij}=\ddot{t}_{ij}$.

Whenever the question is raised whether the emitted signal of a source would possibly be detectable by earth-based facilities, such as e.g. LIGO, one is particularly interested in the detector dependent characteristical frequencies $f_{\rm c}$ and amplitudes $h_{\rm c}$ (calculated according to Eq. (31) in Hawking & Israel 1989) and the ``signal-to-noise ratios'' (SNR). We estimate the SNR for optimal filtering searches as (Müller et al. 2004; Flanagan & Hughes 1998)

  
      $\displaystyle \left(\frac{S}{N}\right)^{2}_{\rm optimal~ filter}$ = $\displaystyle 4\int \frac{\vert\hat{h}(\nu)\vert^{2}}{S_{h}(\nu)}{\rm d}\nu$ (14)
  = $\displaystyle 4\int \left[\frac{\vert\hat{h}(\nu)\vert\nu^{1/2}}{h_{\rm rms}(\nu)}\right]^2{\rm d}\log\nu,$ (15)

where $\hat{h}(\nu)$ is the Fourier transform of the gravitational wave amplitude, $ h(\nu) = h_{+}(\nu) + h_{\times}(\nu) $, and $S_{h}(\nu) [1/{\rm Hz}]$ is the power spectral density of strain noise in the detector and $h_{\rm rms} [{\rm Hz}^{-1/2}]\equiv\sqrt{S_{h}(\nu)}$ stands for the rms noise level of the detector at a given frequency $\nu$. The Fourier spectra were normalised by using Parseval's theorem.

Below, we determined the SNRs for an optimal orientation of detector and source. As detector sensitivity we used the present performance of one single LIGO instrument and the improved, possible future performance of an Advanced LIGO (AdvLIGO) detector (Shoemaker 2007). Note that Advanced LIGO possesses several adjustable frequency responses. While the ``burst'' selection provides a broad range of nearly maximum sensitivity (optimal for model s15g), an ``nsns''-tuned instrument is likely to be used for sources that radiate at lower frequencies (model s15h).

3 Results


  \begin{figure} \par\includegraphics[width=8.8cm]{8577fig2.eps} \end{figure} Figure 2: Time evolution of the central densities of the models s15g and s15h. At core bounce the following values are assumed: $\rho _{\rm c}=3.8 \times 10^{14}$ g/cm3 for model s15g and $\rho _{\rm c}=4.5 \times 10^{14}$ g/cm3 for model s15h. The density decrease in the postbounce phase is most likely due to numerical dissipation of steep density gradients in the neutron star.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm]{8577fig3.eps} \end{figure} Figure 3: Time evolution of the quadrupole amplitudes $A_{\rm +I}$, $A_{\rm\times I }$, $A_{\rm +II}$, and $A_{\rm\times II}$ for model s15g.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm]{8577fig4.eps} \end{figure} Figure 4: Time evolution of the quadrupole amplitudes $A_{\rm +I}$, $A_{\rm\times I }$, $A_{\rm +II}$, and $A_{\rm\times II}$ for model s15h.
Open with DEXTER

The computed quadrupole wave amplitudes of models s15g and s15h are displayed in Figs. 3, 4, and 6. A quantitative overview of key properties and results is given in Table 1 and Fig. 5. Detector dependent quantities are summarised in Table 2 and Fig. 8.

Table 1: Summary of GW-related quantities for the models s15g and s15h.

3.1 Model s15g

Model s15g undergoes a rotational core collapse as we assume an initial central angular velocity of $ \Omega =2\pi $ rad/s. During the stage of contraction, the core spins up massively while becoming oblate. The collapse gets abruptly halted due to the stiffening of the EoS above nuclear densities. Afterwards the core rebounds and drives a hydrodynamical shock wave outwards. These conditions give rise to strong time-dependent variations in the physical quantities that affect the quadrupole tensor (see Eq. (7)). Its behaviour is then directly reflected by the model's gravitational wave signature. Since the core collapse proceeds nearly axisymmetrically, the only GW amplitude considerably driven by the rotationally induced large-scale asymmetries is $A_{\rm +II}$. It exceeds the other s15g-wave trains $A_{\rm +I}$, $A_{\rm\times I }$, and $A_{\rm\times II}$ by 1-2 orders of magnitude, as one can see in Fig. 3 for times around bounce. The initial GW signal is emitted just before core bounce, when the rapid infall of matter and the spin-up of the core is dominant. It shows a prebounce rise. Then, the instantaneous slowdown of matter at core bounce leads to a prominent negative peak, which is followed by a ring-down behaviour that lasts for the first few ms postbounce. This generic type I wave characteristics is displayed in the lower left panel of Fig. 3. Performing a Fourier transform of the GW signal around bounce (-5 ms <t<5 ms), we find a spectrum with a very narrow bandwidth peaking around 893 Hz. This is in good agreement with the recent findings of Müller et al. (2004), but more than 150 Hz higher than in Dimmelmeier et al. (2007) and Ott et al. (2007a). It has been suggested that the difference stems from full relativistic calculations shifting the bounce spectrum to lower frequencies in comparison to the ones using an effective gravitational potential (Dimmelmeier 2007).


  \begin{figure} \par\includegraphics[width=8.8cm]{8577fig5.eps} \end{figure} Figure 5: Time evolution of the GW luminosity in model s15g. In this 80 ms window, about $90\%$ of the energy emitted in GW stems from the $A_{\rm +II}$ contribution around bounce.
Open with DEXTER

After this first and predominantly axisymmetric stage of GW-emission, the occurrence of a low T/|W|-instability revives the gravitational wave signal again around $\approx $20 ms post-bounce (Saijo & Yoshida 2006; Ott et al. 2007a; Watts et al. 2005; Saijo et al. 2003; Ott et al. 2005; Ou & Tohline 2006). Low T/|W| dynamical instabilities are triggered in differentially rotating systems such as neutron stars in situations where the patten speed $\sigma_{\rm p}=\sigma/m$ of an unstable mode m matches the local angular velocity at a point in the star (see Fig. 7), commonly called corotation point. (The modes are, as in Watts et al. (2005), assumed to behave harmonically as $\exp[-i(\sigma t -m\phi)]$, where $\sigma$ is the mode's eigenfrequency.) It permits the azimuthal fluid modes to amplify. This non-axisymmetric process yields a quasi-periodic GW signal with a rather constant time variation, leading to a narrow-band emission at 905 Hz that lasts until the end of our simulation, as one can see, particularly in the upper panels of Fig. 3 for times t>20 ms. The analysis method we use to observe the growth of nonaxisymmetric structures decomposes the density at a fixed radius R and constant z-component into its azimuthal Fourier components as done before e.g. in Ou & Tohline (2006) and Ott et al. (2007a):

\begin{displaymath}\rho (R,z,\phi) =\sum_{m = -\infty}^{\infty}C_{m}(R,z){\rm e}^{{\rm i}m\phi}, \end{displaymath} (16)

where the complex Fourier amplitudes are defined by

\begin{displaymath}C_{m} = \frac{1}{2\pi} \int_{0}^{2\pi}\rho(R,z,\phi){\rm e}^{-{\rm i}m\phi}{\rm d}\phi. \end{displaymath} (17)

In Fig. 6 the normalized mode amplitudes Am=|Cm|/C0 are monitored to measure the growth of unstable modes. In our model s15g, we find $m=\{1,2,3\}$-modes being triggered, with the so-called m=2 bar-mode growing fastest. Furthermore we state that the $m=\{1,2,3\}$ modes all possess the same pattern speed. The close relation between the m=2 bar-mode instability and the emission of gravitational waves can be seen in the following features. First, the sudden onset of GW emission along the pole, which must be due to completely nonaxisymmetric dynamics, coincides with the amplitude of the m=2 mode reaching approximately the absolute amplitude of the m=4 mode caused by the grid. Second, the dominant frequency of emission corresponds perfectly to the eigenfrequency of the m=2 mode. Finally, the two GW-polarisations + and $\times$ are phase-shifted by $\pi /2$, as one would expect of a perfect, monochromatic GW-source such as a rotating bar. These findings in the context of the low T/|W| instability and supernova dynamics stand in remarkable agreement with the recent ones of Ott et al. (2007a). For low $\beta$-unstable models similar to s15g, they find narrow-band GW emission at $\approx $920-930 Hz. The main difference to our calculations is the point that the dominant mode that was found in those computations was the m=1 mode. As a closing remark to model s15g, we state that the time evolution of the energy emitted by gravitational radiation (see Fig. 5) fits the behaviour of the waves. It demonstrates a high peak around bounce at $1.3\times10^{51}$ erg/s followed by a ringdown and an oscillating renaissance at about 1047-1048 erg/s for times t>20 ms.

Table 2: Summary of detector-dependent quantities.


  \begin{figure} \par\includegraphics[width=15.98cm]{8577fig6.eps} \end{figure} Figure 6: The upper panel displays the GW emission of the A+ and the $A_{\times }$-amplitudes along the pole. The lower panel shows the normalized mode amplitudes Am for $m=\{1,2,3,4\}$ extracted at a radius of 20 km. The sudden onset of the nonaxisymmetric GW-signal along the pole as soon as the m=2 mode amplitude almost reaches the absolute size of the m=4 grid mode background. Both A+ and $A_{\times }$ oscillate at the same frequency of 905 Hz, phase shifted by $\pi /2$, as one would expect from a rotating bar.
Open with DEXTER

3.2 Model s15h

The slow-rotating model s15h undergoes a quasi-spherically symmetric core collapse, consequently showing fairly weak type I amplitudes around bounce (see Fig. 4). On the other hand, the marginally present centrifugal forces allow a core bounce at higher central densities ( $4.5 \times10^{14}$ g/cm3) than in the previously discussed model ( $3.8\times10^{14}$ g/cm3) as one can see in Fig. 2. The general dependence of $A_{\rm +II}$ on rotationally induced asymmetries gets obvious in the following feature: while the $A_{\rm +II}$ amplitude from s15h at bounce is more than one order of magnitude smaller as in s15g, all other amplitudes, namely $A_{\rm +I}$, $A_{\rm\times I }$, and $A_{\rm\times II}$ are of the same order of magnitude as their counterparts from model s15g, emphasising a general feeble subordination to rotationally caused aspherities around bounce.

Prompt convection, based on the negative entropy gradient left behind the stalling shock starts just a bit before 10 ms after bounce. It leads to a convectively driven rise in the GW signal. Its characteristic is dominated by the stochastic process of convective mass motion: neither a clear signal type nor any correlation between the two polarisations as can be found here, as one expects from this kind of matter motion. However, there is a peak in the wave train at $\approx $7 ms, most visible in $A_{\rm +II}$, which we could not attribute to a physical feature of core collapse. It is most likely the result of a grid-alignment effect at a radius of 30-50 km. When the fluid first breaks its spherical symmetry in this region, this leads to a numerical quadrupole moment in the fluid aligned with the main coordinate axes. The overall spectra arising from the s15h-GW amplitudes is qualitatively and quantitatively different to the one of s15g, because it ranges from some Hz to about 1000 Hz with major contributions below 500 Hz, clearly dominated by post-bounce convective motions (see Fig. 8). The signal lasting from later times (t>20 ms) mainly takes place at a frequency range from roughly 250 Hz down to about about 10 Hz. This broad band emission of GW mirrors the wide spread and inhomogeneous velocity distribution of the convective motion. The total energy emitted in GW is nearly three orders of magnitude lower than in the previously discussed model (see Table 1), and the maximum luminosity reaches only $\approx $ ${3\times10^{47}}$ erg/s at core bounce. The net size of our convectively driven amplitudes (some centimetres), as well as the frequency band of emission and the total amount of energy released in gravitational radiation, fit the results of Müller et al. (2004) qualitatively well, where they performed 2D core collapse simulations of slowly rotating progenitor models with similar input physics.


  \begin{figure} \par\includegraphics[width=8.8cm]{8577fig7.eps} \end{figure} Figure 7: Angular velocity profiles of model s15g ( upper profile) and s15h ( lower profile) along the positive x-axis in the equatorial plane at $t-t_{\rm b}\approx 6$ ms (S15g) and $t=t_{\rm b}$ (s15h), respectively. The m=2 mode pattern speed of s15g is given by the horizontal line. First, both models' innermost $\approx $10 km are nearly in perfect solid body rotation, and second, model s15g is in corotation at a radius of about 10 km, spinning with some 2850 [rad/sec].
Open with DEXTER

3.3 Detectability

As in Müller et al. (2004), in Fig. 8 we compare the quantity $\vert\hat{h}(\nu)\vert\nu^{1/2}$ in Eq. (15), evaluated at a Galactic distance of 10 kpc for a single LIGO detector's $h_{\rm rms}$. Model s15g's signal enhancement in the equatorial plane around 900 Hz finds its physical origin in the core bounce, while the polar signal around the same frequency is dominated by the m=2 dynamical instability. The computed SNRs, the characteristic frequencies $f_{\rm c}$, and amplitudes $h_{\rm c}$ are displayed in Table 2. The results indicate a fair chance of detecting model s15g within our Galaxy at the present performance of LIGO. The prospect of observing a narrow-band and long-lasting GW signal from the nonaxisymmetric dynamical instability (s15g) is enhanced, increasing the detection limit for a facility running at current sensitivity up to about 100 kpc. Furthermore, Fig. 8 reveals that the SNR of the nearly non-rotating supernova model s15h, where the bounce signal is only marginally present, is most probably to small to be traced by LIGO, although the peak sensitivity of current detectors lies in the frequency range that our results show is dominated by convection. The increased detector sensitivities of Advanced LIGO should allow the catching of events such as the rotational core collapse model s15g at distances around 1 Mpc and the nearly spherically symmetric s15h at about 100 kpc.

4 Uncertainties in the underlying models

Many branches of physics contribute to the models that enable the prediction of GW signals from core-collapse supernovae. Therefore, it is a difficult task to objectively judge the relative importance of uncertainties in GW signal templates. Perhaps we can start by distinguishing features that appear in generic models from features that appear in deterministic models. A feature of generic models would be a feature that appears for all representative models out of a given class, while a feature of a deterministic model would only appear in a specific simulation. We believe that most current supernova models should be seen as generic models that aim to represent a class of stars, rotation rates, etc., instead of a specific event. Hence, we should concentrate on the generic features of the predicted GW signals and defer the discussion of deterministic features to future, improved supernova models. Among the generic features we can further distinguish qualitative statements and quantitative statements. For example, the finding that the A $ _{\rm +II} $ quadrupole from a rotating collapse model shows a pre-bounce rise before a high negative peak as shown in Fig. 3 is a generic, qualitative statement (Dimmelmeier et al. 2007), while the information that the Fourier transform peaks at 893Hz is a quantitative statement.


  \begin{figure} \par\includegraphics[width=16.58cm]{8577fig8.eps} \end{figure} Figure 8: Spectral energy distribution of the GW signals at a distance of 10 kpc in comparison with the present LIGO strain sensitivity and the possible performance of Advanced LIGO (broadband tuning). The strain sensitivity curves have kindly been provided by Shoemaker (2007). The left panel shows the result for the rotational model s15g and the right panel shows the result for the non-rotational model s15h. Optimal orientation between source and detector is assumed and the polarizations are combined according to Eq. (15). The dashed lines show the GW signal for an observer on the rotational axis of the source, the solid lines show the GW signal for an observer in the equatorial plane.
Open with DEXTER

It is almost impossible to draw the line between generic and deterministic features based on a single simulation. But even if a full series of models is available, it is difficult to exclude some generic features of this series not being influenced by a generic limitation of the underlying numerical algorithms. Hence, it is very important that different groups do not just develop the one and best-suited code for a given problem, but rather a variety of different numerical approaches so that the results can be compared so as to reveal the common generic features of the GW signal from core-collapse simulations.

The last stage of the evolution of a massive star proceeds through very different phases. First, there is the stellar core collapse and bounce at nuclear density. Then there is a possibly extended accretion phase that eventually leads to the supernova explosion. These two phases involve different conditions of matter and pose different challenges to the numerical modelling. The quality and reliability of the models depends very strongly on the investigated phase. In the following, we discuss the uncertainties of the models for each phase separately.

4.1 Core-collapse and bounce

The models of stellar core collapse and bounce provide a link between progenitor star properties and the first strong emergence of a GW signal. Many previous studies analysed the GW signal based on idealised input physics in order to study the qualitative dependence of the core-bounce signal on progenitor star properties (e.g. Zwerger & Müller 1997). For example, polytropic equations of state are used and the influence of neutrino interactions on the collapse dynamics are ignored. This approach produced a variety of GW bounce signals as a function of parameters like the angular momentum profile of the progenitor star, the adiabatic index of the equation of state, or the presence of strong magnetic fields. As there is no GW emission in spherical symmetry, these studies have always been conducted in a multi-dimensional setting.

However, the input physics relevant during the collapse of the stellar core is very rich on the microscopic level (Martínez-Pinedo et al. 2006), and this dynamical phase is very sensitive to small perturbations. The evolution of the collapsing core depends significantly on the adiabatic index of the equation of state, general relativistic effects, the electron capture rates on free protons and nuclei, and the coherent scattering opacities of neutrinos off the different nuclei. Input physics improvements have been explored in the context of a long-term modelling effort to clarify the supernova explosion mechanism in spherically symmetric models. Through the ever-increasing power of computers, these improvements can now be carried over to the multi-dimensional predictions of GW signals (Müller et al. 2004). Few groups have implemented equations of state that contain nuclear input physics and added our simple and efficient neutrino physics parametrisation scheme to get a more accurate mapping of the progenitor properties onto the expected GW signal from the bounce at nuclear matter densities (Dimmelmeier et al. 2007; Cerda-Duran et al. 2007). It was possible to show by two- and three-dimensional general relativistic simulations that only type I signals are expected to occur (Dimmelmeier et al. 2007; Ott et al. 2007b), which we independently confirm by the models presented in this paper.

On the other hand, the uncertainties in the progenitor properties and the equation of state above nuclear density are rather large. Hence, even with the most accurate numerical scheme it is not possible to calculate a definitive GW signal from core bounce. The improved models of the latest generation help to accurately map progenitor properties and equation of state properties to a potentially observable GW signal. This may lead to constraints from future GW observations. For example, the amplitude and timing of the GW signal with respect to the neutrino signal strongly depends on the rotation rate of the inner stellar core. It could also depend on the size of the collapsing core (weak interactions) and on asymmetric perturbations induced by convection in the stellar envelope at the time of collapse. The peak frequency of the Fourier-transformed signal contains information on the compressibility of nuclear matter at bounce. Hence, it should be possible to compare GW wave signals for different physical equations of state. One should also investigate a potential impact of strong magnetic fields on the asymmetries of the collapse that might become visible in the GW signal (Obergaulinger et al. 2006; Kotake et al. 2004).

4.2 Postbounce phase

The homologously collapsed stellar core forms a neutron star before it probably enters an extended accretion phase. During this postbounce phase layers, from the outside of the original iron core fall into the standing accretion shock that results from core bounce. They are shock-heated and dissociated, and settle on the protoneutron star. For several reasons, this phase is much more difficult to accurately capture in a numerical model than the collapse phase. One problem is that the neutrinos not only stream away from the surface of the neutron star, but their fractional absorption behind the standing accretion shock also leads to an essential feedback to the accretion dynamics. The hot layers around the protoneutron star show several three-dimensional fluid instabilities that make an accurate treatment of the radiative coupling between the deleptonisation of the protoneutron star and the accreting layers very delicate. Additionally, one has to keep in mind that the density contrast between the centre of the neutron star, and these hot layers behind the standing accretion shock may easily exceed five orders of magnitude, which imposes severe time-step constraints on codes that treat the whole domain consistently for the extended evolution time until the supernova explosion is thought to be launched. As the time step of the hydrodynamics part is limited by the ratio of the zone width to the signal speed, larger time steps can be taken if the resolution is low in the regions of highest sound speed. As described above, we can currently perform simulations with an equidistant resolution of 1 km. This resolution is conveniently high in the outer layers behind the standing accretion shock, but rather low at the surface and in the interior of the protoneutron star. We avoid numerical artifacts from the steep density gradient at the surface of the protoneutron star by analytically considering the hydrostatic density gradient in the TVD advection scheme. However, the low resolution might still suppress an efficient coupling of low protoneutron star modes to higher dynamical modes (Weinberg & Quataert 2008). A more flexible grid for our code is under development in order to adapt the resolution to the local conditions.

For these reasons, we consider the technical uncertainties in the postbounce phase to be on a similar level to the uncertainties of the input physics, which of course continue to be present after core bounce. Only very few studies have predicted gravitational waves based on a postbounce model that includes sophisticated neutrino physics in axisymmetry (e.g. Müller et al. 2004). Neutrino physics approximations have been used in earlier three-dimensional models e.g. Fryer et al. (2004). In our current models, we continue the hydrodynamical simulation to the postbounce phase, but are aware that the neutrino physics parametrisation cannot handle the neutronisation burst that sets in shortly after bounce, and it does not feature the deep electron fraction trough that develops behind the shock due to the copious electron captures. Furthermore, the neutrino parametrisation scheme does not feature any neutrino heating. Hence, it is important to distinguish the results from the collapse phase and bounce, where we believe that the neutrino physics parametrisation is a viable and reasonably accurate approach, from the dynamics of the postbounce phase, where the neutrino physics parametrisation is not sufficiently accurate. The data from the postbounce evolution in our models has to be understood as an idealised exploration of potential GW features after the bounce signal. This is currently the state-of-the-art in three-dimensional simulations, and it is handled the same way in the models of Ott et al. (2007b).

To improve the models also in the postbounce phase, we have developed the isotropic diffusion source approximation (IDSA) for the neutrino transport (Liebendörfer et al. 2007) for future models. The scheme has been implemented and tested in spherical symmetry, but is not yet fully functional in three dimensions. From basic considerations and comparison to the results in axisymmetry with accurate neutrino transport (Müller et al. 2004), we expect that the additional neutrino emission will lead to a more compact neutron star that is accreting matter from an envelope with stronger fluid asymmetries. Compared to our current simulations, one could expect a stronger GW signal in the late postbounce phase with a shorter periodicity in the signal for comparable initial rotation rates of the progenitor. However, this will be explored in more detail with the next generation of postbounce supernova models.

5 Conclusion

In this paper we presented the GW analysis of two 3D MHD core collapse supernova simulations, which differ only in the amount of initial rotation. Model s15h started with an angular velocity of $ \Omega =0.3 $rad/s, model s15g with $ \Omega =2\pi $ rad/s. We incorporated progenitor stars from stellar evolution calculations, spherically symmetric GR effects, the Lattimer-Swesty EoS, magnetic fields and a neutrino parametrisation scheme that is accurate until the first few ms after bounce. The particular choice of the initial angular momentum allowed clear distinctions between gravitational wave features that stem from rotation or nonaxisymmetric motions.

The amplitudes for the three directions and polarisations, $ A_{\rm +I},A_{\rm\times I},A_{\rm\times II} $, are not very sensitive to rotation and show a similar size of several centimetres. Apparently, they initially couple only weakly to rotationally induced large-scale asymmetries in the mass-energy distribution. The only GW amplitude that is strongly correlated to axisymmetric rotation at the time of core bounce turns out to be $A_{\rm +II}$ in the $ \theta =\pi /2 $, $\phi=0$-direction. In model s15g, it exceeds all other amplitudes and in particular the corresponding one from model s15h by more than one order of magnitude. It shows a clear type I characteristics and implies that a rotational core collapse stays axisymmetric in the early postbounce phase, as lately discussed in Ott et al. (2007a). At the bounce stage of simulation s15g, the dominant band of emission peaks around 893 Hz (s15g), which is in good agreement with the recent findings from Müller et al. (2004), but roughly 150 Hz higher than in Dimmelmeier et al. (2007) and Ott et al. (2007a). Within the investigated time window of about 80 ms duration, the channel $A_{\rm +II}$ in model 15g accounts during the bounce and ring-down phases for $\approx $$90\%$ of the total energy release in GW emission.

The models of the later postbounce phase are still affected by significant technical and physical uncertainties. They should be understood as a qualitative outlook for future models that will have to include better neutrino transport and a higher resolution of the protoneutron star. In the postbounce phase, nonaxisymmetric dynamics starts to play an important role. The wave trains are revived either through convective or through nonaxisymmetric instabilities in the protoneutron star. The GW amplitudes from both models show a sustained signal of approximately constant size, but with a different physical origin. In model s15g, the occurrence of a so-called low T/|W| instability of dominant m=2 character around 20 ms post-bounce leads to the prolonged narrow band GW emission at a frequency of 905 Hz, which is in good agreement with Ott et al. (2007a) in all points, except their dominant mode is m=1. On the other hand, the onset of nonaxisymmetric dynamics in model s15h is triggered by convective motions due to a negative entropy gradient. In our current models, the characteristic frequency of emission linked to convective features takes place below 250 Hz. This number may require adaption after the inclusion of more accurate neutrino transport. Furthermore, no sign of a decay, wether due to a dynamical instability or convection, is present in the long-term evolution of the wave trains, as already observed in (Müller et al. 2004; Ott et al. 2007a).

Finally we conclude that gravitational waves from the discussed dynamical features, namely the rotational core bounce and its subsequent ringdown and low T/|W| instability related to a core collapse supernova like model s15g could possibly be detected by a current LIGO instrument within a distance of 10 kpc. For model s15h, it seems more likely that LIGO would only catch signals from later stages of the supernova evolution, meaning frequency distributions caused by convective motion. The future improvement of detector sensitivities in combination with several adjustable frequency responses of Advanced LIGO will give us the opportunity to tune the instruments to a particular event/source and therefore allow looking much deeper into space. This will enhance the chance of observing core-collapse supernovae similar to our models up to distances of 1 Mpc.

Acknowledgements
We would like to thank R. Käppeli and F.-K. Thielemann for supporting this work in numerous ways. Further acknowledgements go to H. Dimmelmeier and C. D. Ott for enlightening discussions and detailed comments and to D. Shoemaker for his very helpful notes on LIGO. The simulations have been carried out on the McKenzie cluster at CITA (Dubinski et al. 2003), the Athena cluster at the University of Basel, and the Swiss Supercomputing Centre CSCS. We acknowledge support by the Swiss National Science Foundation under grant No. 200020-105328/1 and PP002-106627/1.

References

 
Copyright ESO 2008