Gravitational Wave Extraction in Simulations of Rotating Stellar Core Collapse

Gravitational Wave Extraction in Simulations of Rotating Stellar Core Collapse


We perform simulations of general relativistic rotating stellar core collapse and compute the gravitational waves (GWs) emitted in the core bounce phase of three representative models via multiple techniques. The simplest technique, the quadrupole formula (QF), estimates the GW content in the spacetime from the mass quadrupole tensor only. It is strictly valid only in the weak-field and slow-motion approximation. For the first time, we apply GW extraction methods in core collapse that are fully curvature-based and valid for strongly radiating and highly relativistic sources. These techniques are not restricted to weak-field and slow-motion assumptions. We employ three extraction methods computing (i) the Newman-Penrose (NP) scalar , (ii) Regge-Wheeler-Zerilli-Moncrief (RWZM) master functions, and (iii) Cauchy-Characteristic Extraction (CCE) allowing for the extraction of GWs at future null infinity, where the spacetime is asymptotically flat and the GW content is unambiguously defined. The latter technique is the only one not suffering from residual gauge and finite-radius effects. All curvature-based methods suffer from strong non-linear drifts. We employ the fixed-frequency integration technique as a high-pass waveform filter. Using the CCE results as a benchmark, we find that finite-radius NP extraction yields results that agree nearly perfectly in phase, but differ in amplitude by at core bounce, depending on the model. RWZM waveforms, while in general agreeing in phase, contain spurious high-frequency noise of comparable amplitudes to those of the relatively weak GWs emitted in core collapse. We also find remarkably good agreement of the waveforms obtained from the QF with those obtained from CCE. The results from QF agree very well in phase and systematically underpredict peak amplitudes by , which is comparable to the NP results and is certainly within the uncertainties associated with core collapse physics.

04.25.D-, 04.30.Db, 97.60.Bw, 02.70.Bf, 02.70.Hm

I Introduction

Massive stars () end their nuclear burning lives with a core composed primarily of iron-group nuclei embedded in an onion-skin structure of progressively lighter elements. Energy generation has ceased in such a star’s high-density core and relativistically-degenerate electrons provide pressure support against gravity. Silicon shell burning, neutrino cooling, and deleptonization eventually push the core over its effective Chandrasekhar mass. Radial instability sets in, leading to core collapse, accelerated by electron capture and photodisintegration of iron-group nuclei (see, e.g., Bethe (1990); Baron and Cooperstein (1990)).

The collapsing iron core separates into a subsonically collapsing homologous () inner core and supersonically infalling outer core. When the former reaches nuclear density, the nuclear equation of state (EOS) stiffens, dramatically increasing central pressure support and stabilizing the inner core, which, due to its large inertia, overshoots its new equilibrium and then rebounds into the still collapsing outer core, launching the hydrodynamic supernova shock. The acceleration experienced by the inner core in this core bounce is tremendous, leading to the reversal of the collapse velocities of order of its of material on a millisecond timescale.

It was realized early on that the large accelerations encountered in stellar collapse in combination with a source of quadrupole (or higher) order asphericity lead to the emission of a burst of gravitational waves (GWs; see Ott (2009) for a historical overview). Rotation, centrifugally deforming the inner core to oblate shape, is an obvious source of such quadrupole asymmetry and rotating core collapse and bounce is the most extensively studied GW emission process in stellar collapse (see, e.g., Ott et al. (2007a); Dimmelmeier et al. (2008); Scheidegger et al. (2010a, b); Takiwaki and Kotake (2010); Ott et al. (2010) for recent studies and references therein). Alternatively, asymmetries in collapse may arise from perturbations, e.g., due to large convective plumes in the final phase of core nuclear burning, and may lead to GW emission at bounce and/or seed GW-emitting prompt postbounce convection Burrows and Hayes (1996); Fryer et al. (2004); Ott (2009). A multitude of GW emission processes may be active in the postbounce, pre-explosion phase. These include convection/turbulence in the protoneutron star and in the postshock region, nonaxisymmetric rotational instabilities of the protoneutron star, protoneutron star pulsations, instabilities of the standing accretion shock, and asymmetric emission of neutrinos (Ott (2009); Ott et al. (2006); Marek and Janka (2009); Murphy et al. (2009); Kotake et al. (2009); Yakunin et al. (2010) and references therein).

Of the entire ensemble of potential GW emission processes in stellar collapse, rotating core collapse and bounce is arguably the simplest and yields the cleanest signal, depending only on rotation, on the nuclear EOS, and on the mass of the inner core at bounce Dimmelmeier et al. (2008). Moreover, 3D studies have shown that collapsing iron cores with rotation rates in the range of what is physically plausible stay axisymmetric throughout the collapse phase and develop nonaxisymmetric dynamics only after bounce Shibata and Sekiguchi (2005); Ott et al. (2007a); Scheidegger et al. (2010a). Hence, the GW signal of rotating core collapse and bounce is linearly polarized and axisymmetric (2D) simulations are sufficient for its prediction. Unlike postbounce dynamics involving large scale and small scale fluid instabilities of stochastic nature, the GW signal of rotating collapse and bounce can, in principle, be predicted exactly for a given set of initial data. Hence, it has the potential of being used in GW searches using matched-filtering techniques (e.g., Thorne (1987)) or alternative approaches also taking into account detailed signal predictions Röver et al. (2009); Summerscales et al. (2008).

Much progress has been made in recent years in the modeling of rotating core collapse and its GW signature. State-of-the-art simulations are general relativistic (GR) Dimmelmeier et al. (2002a); Dimmelmeier et al. (2005); Shibata and Sekiguchi (2004, 2005); Shibata et al. (2006); Ott et al. (2007a, b); Kuroda and Umeda (2010); Dimmelmeier et al. (2007, 2008) and some studies include magnetic fields Shibata et al. (2006); Obergaulinger et al. (2006); Kuroda and Umeda (2010) or finite-temperature EOS, deleptonization, and progenitors from stellar evolutionary calculations Ott et al. (2007a, b); Dimmelmeier et al. (2007, 2008). These improvements in the physics included in core collapse models provide for a more accurate and reliable dynamics underlying the emission of GWs. The calculation of the GW signal itself, however, is still being carried out predominantly in the slow-motion, weak-field quadrupole approximation (e.g., Thorne (1980)) that is of questionable quality, given the extreme densities and velocities involved in core collapse. The quadrupole formula (QF) “extracts” GWs based on matter dynamics alone, is not invariant under general relativistic gauge transformations, treats the emission region as a point source, and suffers from the fact that the definition of the generalized mass quadrupole moment is not unique in GR.

In GR, the GW content of a spacetime can be extracted by means of the perturbative Regge-Wheeler-Zerilli-Moncrief (RWZM) formalism Regge and Wheeler (1957); Zerilli (1970a, b); Moncrief (1974a) which is gauge invariant to first order or via the Newman-Penrose (NP) scalars approach Newman and Penrose (1962); Penrose (1963) which depends on the non-unique choice of the tetrad in which the Newman-Penrose scalars are evaluated. For reliable results, both RWZM and NP require extraction in the wave zone Thorne (1980) at coordinate radii many wavelengths from the source, but even there, coordinate ambiguities exist. The latter are removed only when GWs are extracted at future null infinity (, see Newman and Penrose (1962); Penrose (1963)), where space is asymptotically flat.

Shibata & Sekiguchi Shibata and Sekiguchi (2003) have used simulations of an oscillating polytropic neutron star model to compare QF and finite-radius RWZM results. For the same basic system, Baiotti et al. Baiotti et al. (2009) compared QF, finite-radius RWZM, and finite-radius NP GW extraction with each other and with results from a 1D perturbation analysis. Both studies found that in the context of neutron star oscillations, the phase of the waveforms obtained with the quadrupole approximation agrees exceptionally well with that of the RWZM and NP extraction methods. Shibata & Sekiguchi, using their particular choice of the generalized quadrupole moment, found a systematic underprediction of the GW amplitudes by the QF. Baiotti et al. Baiotti et al. (2009), who studied multiple incarnations of the QF, found either underprediction or overprediction of the amplitude, both by up to , depending on the particular choice of QF. Nagar et al. Nagar et al. (2005) studied the performance of RWZM and QF-based GW extraction from oscillating polytropic tori and found qualitatively similar results, and quantitative differences in amplitudes and integrated emitted energies between and , again depending on the choice of quadrupole moment.

RWZM and NP GW extraction and comparisons with the QF approximation for GWs emitted in core collapse spacetimes have proven difficult. On the one hand, the emitted GWs are weak: Typical strain amplitudes are , where D is the distance to the source, and typical emitted energies are of order Ott (2009), many orders of magnitude lower than what is expected, for example, from double neutron star coalescence Oechslin et al. (2007) or binary black hole mergers Reisswig et al. (2009a); Lousto et al. (2010). On the other hand, the GWs have typical frequencies of and corresponding wavelengths of , hence require extraction at large coordinate radii where the grid resolution of core collapse simulations is typically too low to allow extraction of the relatively low-amplitude GWs emitted in core collapse (see, e.g., the discussion in Ott (2006)). Shibata & Sekiguchi, in Shibata and Sekiguchi (2005), were able to extract GWs with the RWZM formalism from an extreme core collapse model that developed a rotationally-induced large-scale nonaxisymmetric deformation after bounce, emitting GWs with . For this model, they found that the QF accurately predicts the GW phase, but underestimates the strain amplitude by . Due to the aforementioned difficulties, these authors were unable to compare RWZM with QF for more moderate, axisymmetric models. Cerdá-Durán et al. Cerdá-Durán et al. (2005) performed core collapse simulations using a second-order post-Newtonian (2PN) extension of the conformal-flatness approximation to GR. Exploiting an approximate relationship of the non-conformal 2PN part of the metric to its GW part Cerdá-Durán et al. (2005), they were able to extract GWs from their 2PN metric in standard axisymmetric rotating core collapse models. They found very close agreement (to a few percent in strain amplitude) between QF and 2PN GW signals for almost all considered collapse models. Siebel et al. Siebel et al. (2003) performed nonrotating axisymmetric core collapse simulations by employing evolutions based on a fully general relativistic null cone formalism. They added nonspherical perturbations to the star, leading to the emission of GWs which they were able to extract with the Bondi news function at . Comparisons to the QF suggested a significant discrepancy in amplitude and frequency from the more reliable Bondi news result.

The results of Shibata & Sekiguchi Shibata and Sekiguchi (2005) and of Cerdá-Durán et al. Cerdá-Durán et al. (2005) provide some handle on the performance of the QF approximation in core collapse spacetimes. The former study, while being performed in full GR, considered only a single extreme model. In addition, the authors were forced to extract GWs with RWZM at too small radii for completely reliable results. The latter study, while considering a broader ensemble of models, was restricted to 2PN without considering full GR, leaving room for doubts about the quality of their GW extraction technique. Finally, the results of Siebel et al. Siebel et al. (2003) were limited to axisymmetry without rotation and are unreliable in the presence of strong shocks Siebel et al. (2003).

In this study, we readdress GW extraction from rotating core collapse spacetimes. We perform GR hydrodynamics simulations of rotating core collapse, for the first time in the core collapse context extracting GWs with RWZM, NP, and multiple QFs and comparing the results of these methods. In addition, and also for the first time in the present context, we utilize the Cauchy-Characteristic Extraction (CCE) approach Winicour (2005); Bishop et al. (1997); Reisswig et al. (2009b); Reisswig et al. (2010); Babiuc et al. (2010) that propagates the GW information to for completely gauge independent and unambiguous GW extraction.

In choosing our models set, we are guided by Cerdá-Durán et al. Cerdá-Durán et al. (2005), and draw precollapse configurations from the set of Dimmelmeier et al. (2002a). These models are GR -polytropic iron cores in rotational equilibrium and we evolve them with an analytic hybrid polytropic/-law EOS used in many previous studies of rotating core collapse Zwerger and Müller (1997); Dimmelmeier et al. (2002a); Obergaulinger et al. (2006); Shibata and Sekiguchi (2004, 2005); Shibata et al. (2006); Cerdá-Durán et al. (2005). For physically accurate GW signal predictions to be used in GW data analysis, a microphysically more complete treatment is warranted. Fortunately, recent results of studies employing such modeling technology (e.g., Ott et al. (2007a, b); Dimmelmeier et al. (2007, 2008); Abdikamalov et al. (2010)) show that, with a proper choice of EOS parameters, hybrid-EOS models are able to qualitatively and to some extent quantitatively reproduce the GW signals obtained with the much more complex and computationally intensive microphysical studies. Hence, for the purpose of this study, we resort to the simpler hybrid-EOS models.

Our simulations employ the open-source Zelmani GR core collapse simulation package Ott et al. (2009) that is based on the Cactus Computational Toolkit Goodale et al. (2003); cac () and the Einstein Toolkit ein (). While using the full GR formalism, we limit our simulations to an octant of the 3D cube, using periodic boundary conditions on two of the inner faces of the octant and reflective boundary conditions on the third face. This limits 3D structure to even and that are multiples of , which is not a limitation for the current study, since rotating core collapse and the very early postbounce evolution are likely to proceed nearly axisymmetrically Ott et al. (2007a); Scheidegger et al. (2008); Scheidegger et al. (2010a). We note that, even though the GW signal in rotating core collapse is dominated by the ’+’ polarization mode, there is no reason to expect different behavior for other GW multipoles or polarizations and our results should translate to the non-axisymmetric case.

The results of our simulations indicate that NP extraction yields results that agree well with those obtained from the most sophisticated CCE method. We observe differences in amplitude of , depending on the model, while the agreement in phase is nearly perfect. We also find that the RWZM formalism yields unphysical high-frequency signal components that make this method less suitable for core collapse simulations where the signal is very weak. Finally, we note that the quadrupole approximation yields surprisingly close results to those obtained from CCE. While the phases nearly perfectly agree, the amplitude shows differences of .

This paper is structured as follows. In Sec. II, we discuss our methodology, initial data, and EOS details. Section III discusses the various GW extraction methods that we employ. In Sec. IV, we present our results and discuss them in detail. Finally, in Sec. V, we summarize and review our findings.

Ii Methods

We adopt the ADM foliation of spacetime York (1979a). All equations assume unless noted otherwise. In the following, Latin indices run from 1 to 3 while Greek ones run from 0 to 3. We adhere to abstract index notation. is the 4-metric, is the 3-metric and the extrinsic curvature.

ii.1 Infrastructure and Mesh Refinement: Cactus and Carpet

Our code uses the Cactus Computational Toolkit Goodale et al. (2003); cac () to manage the complexity inherent in large software projects. Cactus is an open source high-performance computing environment designed for scientists and engineers; its modular structure enables parallel computation across different architectures, and facilitates collaborative code development between different groups. Indeed, our code uses a set of components of the public EinsteinToolkit Schnetter (2008); ein (), a community project developing and supporting open software for relativistic astrophysics, such as e.g. the curvature and hydrodynamics evolution methods described below. Many improvements made in the course of the research for this paper were contributed back to the community.

In particular, Cactus allows us to clearly separate between physics components and computational components in our code. Distributed memory parallelism in Cactus is provided by a driver component which implements the data structures discretising the manifold on which the computational state vector lives. In our case, this is the Carpet driver Schnetter et al. (2004, 2006); Carpet () providing block-structured adaptive mesh refinement (AMR) and multi-block discretization. Carpet parallelizes using a hybrid approach combining MPI and OpenMP, where inter-node communication is handled via MPI and intra-node communication via OpenMP or via MPI, depending on the particular system and on details of the simulation setup.

Carpet implements Berger-Oliger style AMR Berger and Oliger (1984), where the fine grids are aligned with coarse grids, refined by factors of two. Carpet also implements subcycling in time, where finer grids take two time steps for every coarse grid step. The latter greatly improves efficiency, but also introduces significant complexity into the time evolution method. The refined regions can be chosen and modified arbitrarily, which we use here to add additional, finer levels during evolution as successively higher resolutions are required to capture the collapse. This is described in more detail in Ott (2006).

We use fifth-order accurate spatial interpolation for spacetime variables and third order essentially non-oscillatory interpolation for hydrodynamics variables. Time interpolation, which is necessary to provide boundary conditions to fine levels at times where there is no coarse level, is second-order accurate. We apply no time refinement between levels 3 and 4, which corresponds to reducing the Courant-Friedrichs-Lewy factor on levels 3 and coarser by a factor of 2. This increases the accuracy on level 3 where we extract gravitational waves. In total, we use 9 refinement levels (including the base grid), an outer boundary radius of () and a finest zone size of in our baseline resolution.

ii.2 Curvature Evolution: McLachlan

Evolution System

We evolve the full Einstein equations in a split (a Cauchy initial boundary value problem) York (1979b), using the BSSN formulation Alcubierre et al. (2000), a slicing Alcubierre et al. (2003), and -driver shift condition Alcubierre et al. (2003). This leads to the following set of evolved variables:


Our exact evolution equations are as described by Eqs. (3) to (10) of Brown et al. (2009), which we list here for completeness:


with the momentum constraint damping constant set to . The stress energy tensor is incorporated via the projections


We have introduced the notation . All quantities with a tilde refer to the conformal 3-metric , which is used to raise and lower indices. In particular, and refer to the covariant derivative and the Christoffel symbols with respect to . The expression denotes the trace-free part of the expression inside the parentheses, and we define the Ricci tensor contributions


This is a so-called -variant of BSSN. The evolved gauge variables are lapse , shift , and a quantity related to the time derivative of the shift. The gauge parameters , , , and are determined by our choice of a slicing:


and -driver shift condition:


The expression attenuates the -driver depending on the radius as described below.

The -driver shift condition is symmetry-seeking, driving the shift to a state that renders the conformal connection functions stationary. Of course, such a stationary state cannot be achieved while the metric is evolving, but in a stationary spacetime the time evolution of the shift and thus that of the spatial coordinates will be exponentially damped. This damping time scale is set by the gauge parameter (see Eq. 23) which has dimension (inverse time). As described, e.g., in Müller and Brügmann (2010); Schnetter (2010), this time scale may need to be adapted in different regions of the domain to avoid spurious high-frequency behavior in regions that otherwise evolve only very slowly, e.g., far away from the source.

Here we use the simple damping mechanism described in Eq. (12) of Schnetter (2010), which is defined as


with a constant defining the transition radius between the interior, where , and the exterior, where falls off as . Eq. 23 describes how appears in the gauge parameters. In this paper we use ().

We implement the above BSSN equations and gauge conditions in the McLachlan code Brown et al. (2009); ES- () which is freely available as part of the EinsteinToolkit. McLachlan is auto-generated from the definition of the variables and equations in the Mathematica format by the Kranc code generator Lechner et al. (2004); Husa et al. (2006); kra (). Kranc is a suite of Mathematica packages comprising a computer algebra toolbox for numerical relativists. Kranc can be used as a “rapid prototyping” system for physicists or mathematicians handling complex systems of partial differential equations, and through integration into the Cactus framework one can also produce efficient production codes.

We use fourth-order accurate finite differencing for the spacetime variables and add a fifth-order Kreiss-Oliger dissipation term to remove high frequency noise. We use a fourth-order Runge-Kutta time integrator for all evolved variables.

Initial Conditions

We set up our initial condition from the ADM variables , , lapse , and shift , as provided by the initial data discussed in Sec. II.4. From these we calculate the BSSN quantities via their definition, setting , and using cubic extrapolation for at the outer boundary. This extrapolation is necessary since the are calculated from derivatives of the metric, and one cannot use centered finite differencing stencils near the outer boundary. We assume that one could instead also use one-sided derivatives to calculate on the boundary.

The extrapolation stencils distinguish between points on the faces, edges, and corners of the grid. Points on the faces are extrapolated via stencils perpendicular to that face, while points on the edges and corners are extrapolated with stencils aligned with the average of the normals of the adjoining faces. For example, points on the edge are extrapolated in the direction, while points in the corner are extrapolated in the direction. Since several layers of boundary points have to be filled for higher order schemes (e.g., three layers for a fourth order scheme), we proceed outwards starting from the innermost layer. Each subsequent layer is then defined via the points in the interior and the previously calculated layers.

Boundary Conditions

During time evolution, we apply a Sommerfeld-type radiative boundary condition to all components of the evolved BSSN variables as described in Alcubierre et al. (2000). The main feature of this boundary condition is that it assumes approximate spherical symmetry of the solution, while applying the actual boundary condition on the boundary of a cubic grid where the face normals are not aligned with the radial direction. This boundary condition defines the right hand side of the BSSN state vector on the outer boundary, which is then integrated in time as well, so that the boundary and interior are calculated with the same order of accuracy.

The main part of the boundary condition assumes that we have an outgoing radial wave with some speed :


where is any of the tensor components of evolved variables, the value at infinity, and a spherically symmetric perturbation. Both and depend on the particular variable and have to be specified. This implies the following differential equation:


where . The spatial derivatives are evaluated using centered finite differencing where possible, and one-sided finite differencing elsewhere. We use second order stencils in our implementation.

In addition to this main part, we also account for those parts of the solution that do not behave as a pure wave, e.g., Coulomb type terms caused by infall of the coordinate lines. We assume that these parts decay with a certain power of the radius. We implement this by considering the radial derivative of the source term above, and extrapolating according to this power-law decay.

Given a source term , we define the corrected source term via


where is the normal vector of the corresponding boundary face. The spatial derivatives are evaluated by comparing neighbouring grid points, corresponding to a second-order stencil evaluated in the middle between the two neighbouring grid points. We assume a second-order decay, i.e., we choose .

As with the initial conditions above, this boundary condition is evaluated on several layers of grid points, starting from the innermost layer. Both the extrapolation and radiative boundary condition algorithms are implemented in the publicly available NewRad component of the Einstein Toolkit.

This boundary condition is only a coarse approximation of the actual decay behavior of the BSSN state vector, and it does not capture the correct behavior of the evolved variables. However, we observe that this boundary condition leads to stable evolutions if applied sufficiently far from the source. Errors introduced at the boundary (both errors in the geometry and constraint violations) propagate inwards with the speed of light Brown et al. (2009). Gauge changes introduced by the boundary condition, which are physically not observable, propagate faster, with a speed up to for our gauge conditions.

ii.3 General-Relativistic Hydrodynamics: GRHydro

We employ the open-source GR hydrodynamics code GRHydro that is part of the EinsteinToolkit ein () and is an updated version of the code Whisky described in Baiotti et al. (2005).

The equations of ideal GR hydrodynamics evolved by GRHydro are derived from the local GR conservation laws of mass and energy-momentum,


where denotes the covariant derivative with respect to the 4-metric. is the mass current with the 4-velocity and the rest-mass density . is the stress-energy tensor. The quantity is the specific enthalpy, is the fluid pressure and is the specific internal energy.

We choose a definition of the 3-velocity that corresponds to the velocity seen by an Eulerian observer at rest in the current spatial 3-hypersurface York (1983),


where is the Lorentz factor. In terms of the 3-velocity, the contravariant 4-velocity is then given by


and the covariant 4-velocity is


The GRHydro scheme is written in a first-order hyperbolic flux-conservative evolution system for the conserved variables , , and in terms of the primitive variables ,


where is the determinant of . The evolution system then becomes




Here, and are the 4-Christoffel symbols. The above equations are solved in semi-discrete fashion. The spatial discretization is performed by means of a high-resolution shock-capturing (HRSC) scheme employing a second-order accurate finite-volume discretization. We make use of the Marquina flux formula for the local Riemann problems and piecewise-parabolic cell interface reconstruction (PPM). For a review of such methods in the GR context, see Font (2008). The time integration and coupling with curvature are carried out with the Method of Lines Hyman (1976).

ii.4 Equation of State and Initial Stellar Models

For the purpose of this study, we employ the simple analytic hybrid EOS Janka et al. (1993); Dimmelmeier et al. (2002a) that combines a 2-piece piecewise polytropic pressure with a thermal component , i.e., . To model the stiffening of the EOS at nuclear density , we assume that the polytropic index jumps from below nuclear density to above. As detailed in Dimmelmeier et al. (2002b), it is possible to construct an EOS that is continuous at ,


In this, denotes the total specific internal energy which consists of a polytropic and a thermal contribution. [cgs] is the polytropic constant for a polytrope of relativistic degenerate electrons at . The thermal index corresponds to a mixture of relativistic () and non-relativistic () gas. This EOS mimics the effects of the stiffening of the physical EOS at and can handle the significant thermal pressure contribution introduced by shock heating in the postbounce phase. Provided appropriate choices of EOS parameters (e.g., Dimmelmeier et al. (2007)), the hybrid EOS leads to qualitatively correct collapse and bounce dynamics. Consequently, this leads to GW signals that are similar in morphology, characteristic frequencies and amplitudes to those computed from more microphysically complete simulations Ott et al. (2007a); Dimmelmeier et al. (2008); Abdikamalov et al. (2010).

We employ () polytropes in rotational equilibrium generated via Hachisu’s self-consistent field method Komatsu et al. (1989a, b) that not only provides fluid, but also spacetime curvature initial data. The polytropes are set up with the rotation law discussed in Zwerger and Müller (1997); Dimmelmeier et al. (2002a) and are parametrized via the differential rotation parameter and the initial ratio of rotational kinetic energy to gravitational binding energy . While being set up as marginally stable polytropes with , during evolution, the initial sub-nuclear polytropic index is reduced to to accelerate collapse. Following previous studies Zwerger and Müller (1997); Dimmelmeier et al. (2002a); Ott et al. (2007b), we use in the super-nuclear regime.

From the initial stellar configurations of Zwerger and Müller (1997); Dimmelmeier et al. (2002a) we draw a subset of three models that cover the range of astrophysically expected GW signals from rotating iron core collapse Dimmelmeier et al. (2008) and accretion-induced collapse Abdikamalov et al. (2010). Our choices have been used previously in a comparison study of full GR and conformally-flat simulations Ott et al. (2007b):

  • Model A1B3G3 is in near uniform rotation with , has , and, once mapped to the evolution grid, uses a sub-nuclear adiabatic index . Its GW signal is of the standard “Type-I” morphology Zwerger and Müller (1997); Dimmelmeier et al. (2002a); Dimmelmeier et al. (2007) and of moderate strength (see Dimmelmeier et al. (2002a); Ott et al. (2007b) for details).

  • Model A3B3G3 also uses . It is strongly differentially rotating, with its initial central angular velocity dropping by a factor of two over . This, in combination with , leads to rapid rotation in the inner core, resulting in a very strong GW signal at core bounce and dynamics that are significantly affected by centrifugal effects. It produces a “Type-I” GW signal with a centrifugally-widened broad peak at core bounce.

  • Model A1B3G5 has the same rotational setup as model A1B3G3, but its sub-nuclear adiabatic index is reduced to . This leads to rapid collapse, to a very small inner core at core bounce, and to a weak “Type-III” GW signal Zwerger and Müller (1997); Dimmelmeier et al. (2002a) akin to that potentially emitted by an accretion-induced collapse event Abdikamalov et al. (2010).

For convenience, key model parameters are summarized in Table 1.

Model Type
A1B3G3 I (weak)
A3B3G3 I (strong)
Table 1: Initial parameters of differentially rotating stellar cores used for the core collapse simulations. The models are described by three quantities: the degree of differential rotation , the ratio of rotational to potential energy, and the sub-nuclear adiabatic index during the collapse. For convenience we also report the wave-signature type of the three models and the mass present on the computational grid.

Iii Gravitational Wave Extraction

iii.1 The Quadrupole Approximation

The quadrupole approximation is the only means of extracting GWs in Newtonian or conformally-flat GR simulations, but has found wide application also in GR simulations of stellar collapse Shibata and Sekiguchi (2004, 2005); Ott et al. (2007b, a).

The coordinate-dependent quadrupole formula estimates the GW strain seen by an asymptotic observer by considering exclusively the quadrupole stress-energy source. It neglects any non-linear GR effects. This approximation is valid strictly only in the weak-field and slow-motion limit Thorne (1980) where spacetime is essentially flat.

The quadrupole formula is given in the transverse-traceless (TT) gauge by




denotes the reduced mass-quadrupole tensor. Since we are working in the weak-field, slow-motion approximation, the placement of tensor indices is arbitrary. is not uniquely defined in GR and the choice of the density variable is not obvious. Following Dimmelmeier et al. (2002a); Shibata and Sekiguchi (2004); Dimmelmeier et al. (2008); Ott et al. (2007b, a), we set , because, (i), this is the conserved density variable in our code, and (ii), is the natural volume element. See Baiotti et al. (2009) for other potential choices and their relative performance for GWs from oscillating polytropes.

The reduced mass-quadrupole tensor can be computed directly from the computed distribution . In order to eliminate the effects of numerical noise when differentiating Eq. (37) twice in time, we make use of the continuity equation to obtain the first time derivative of Eq. (37) without numerical differentiation Finn and Evans (1990); Blanchet et al. (1990); Dimmelmeier et al. (2005),


where we set as defined by Eq. 29. Note that we have switched to contravariant variables in the integrand as these are the ones present in the code. Finally, the remaining time derivative needed for evaluating the quadrupole GW strain (Eq. 36) is performed numerically.

In order to assess the sensitivity of the predicted waves on the particular choice of the velocity variable in Eq. (38), we implement two modified versions. In variant VS, we use Shibata et al.’s definition of the 3-velocity (e.g., Shibata and Sekiguchi (2004)) that differs from ours by a gauge term. In variant PV, we follow Dimmelmeier et al. (2002b) and employ physical velocity components (individually bound to ) that, in Cartesian coordinates, are given by , assuming that the 3-metric is nearly diagonal (which is the case in our gauge; see Sec. II.2).

iii.2 The Regge-Wheeler-Zerilli-Moncrief Formalism

A particular Ansatz for analyzing gravitational radiation in terms of odd and even multipoles in the far-field of the source was originally developed by Regge, Wheeler Regge and Wheeler (1957) and Zerilli Zerilli (1970a, b), respectively. Moncrief subsequently provided a gauge-invariant reformulation Moncrief (1974b) (see Nagar and Rezzolla (2005) for a review). Assuming that, at large distances from the source, the GW content of the spacetime can be viewed as a linear perturbation to a fixed background, we can write


where is the fixed background metric and its linear perturbation. The background metric is usually assumed to be of Minkowski or Schwarzschild form, which we can write as


By splitting the spacetime into timelike and radial and angular parts, it is possible to decompose the metric perturbation into odd and even multipoles, i.e., we can write


The even and odd multipole components are defined according to their behavior under a parity transformation . Odd multipoles transform as while even multipoles transform as . Both multipole components can be expanded in terms of vector and tensor spherical harmonics (e.g., Thorne (1980)).

Given the Hamiltonian of the perturbed Einstein equations in ADM form Arnowitt et al. (1962), it is then possible to derive variational principles for the odd and even-parity perturbations Moncrief (1974b) to give equations of motions that are similar to wave equations with a scattering potential.

The solutions to the odd and even-parity wave equations are given by the Regge-Wheeler-Moncrief and the Zerilli-Moncrief master functions, respectively. The odd-parity Regge-Wheeler-Moncrief function reads


and the even-parity Zerilli-Moncrief function reads


where , and where




These master functions depend entirely on the spherical part of the metric given by the coefficients and , and the perturbation coefficients for the individual metric perturbation components , , , , , , , , and which can be obtained from any numerical spacetime by projecting out the Schwarzschild or Minkowski background Camarda and Seidel (1999). For example, the coefficient can be obtained via


where is the radial component of the numerical metric represented in the spherical-polar coordinate basis, are spherical harmonics, and is the surface line element of the extraction sphere. The coefficient represents the spherical part of the background metric and can be obtained by projection of the numerical metric component on over the extraction sphere


Similar expressions hold for the remaining perturbation coefficients.

The odd- and even-parity master functions Eq. (42) and Eq. (43) can be straight-forwardly related to the gravitational-wave strain and are given by


where are the spin-weight spherical harmonics. We note that this relation is strictly true only at an infinite distance from the source. Since our numerical domain is finite in size, we choose some, ideally large, but finite radius. In Sec. IV.3, we check how well the GWs extracted with the RWZM formalism asymptote with increasing extraction radius.

In the present work, our models exclusively trigger the even-parity master function , and is zero. In this case, we can simplify Eq. (49) and obtain


relating the strain directly to .

iii.3 Newman-Penrose Scalars

Another method for calculating the gravitational waveforms is based on the conformal structure of asymptotically flat spacetimes as established by Bondi, Sachs and Penrose Bondi et al. (1962); Sachs (1962); Penrose (1963). This method is conveniently represented in terms of spin-weighted scalars as introduced by Newman and Penrose Newman and Penrose (1962). In the following we refer to it as NP extraction. According to the peeling theorem Bondi et al. (1962); Sachs (1962), a certain component of the conformal Weyl tensor obeys the slowest fall-off from the source, and hence is identified as outgoing gravitational radiation:


Here, the slowest fall-off is obeyed by the NP scalar , which is defined as1


where is the conformal Weyl tensor associated with the 4-metric and , are part of a null-tetrad Penrose (1963); Newman and Penrose (1962) which satisfies while all other inner products vanish. In addition, this tetrad is related to the 4-metric via . At future null infinity , the topology of the spacetime is a time succession of spheres, . Hence the simplest choice for the null tetrad at is such that it resembles the unit sphere metric. Moreover, the simplest choice for a coordinate system at is given by the Bondi gauge Bondi et al. (1962); Sachs (1962), which makes use of an areal radius coordinate.

In most current numerical relativity simulations, the radiation is computed at a finite radius where the Bondi coordinates are usually not imposed. Rather, we use the gauge as evolved by the slicing and -driver conditions discussed in Sec. II.2. In practice, we impose a simple polar-spherical coordinate system with constant coordinate radius , which does not take into account the background geometry, and hence does not make use of an areal radius. Thus, the gravitational radiation as computed on these coordinate spheres is not measured in the correct gauge, and leads to a systematic error that needs to be assessed. Note that it is principally possible to transform to the correct gauge Lehner and Moreschi (2007).

In our construction of an approximate tetrad, we follow common practice (e.g. Pollney et al. (2009a); Bruegmann et al. (2008); Baker et al. (2002)) and use a triad of spatial vectors obtained via a Gram-Schmidt orthonormalization starting from


where are Cartesian coordinates of the computational grid and is the three-dimensional Levi-Civita symbol. The tetrad is given in terms of this triad and the timelike normal vector by


A straightforward calculation shows that we are thus able to express exclusively in terms of the “3+1” variables according to


where the electric and magnetic parts of the Weyl tensor are defined as Friedrich (1996)


Here the denotes the Hodge dual and is the projection operator. The Gauss-Codazzi equations (see e.g. Alcubierre (2008)) enable us to calculate the electric and magnetic parts from the “3+1” variables according to


In a given numerical simulation, we calculate from Eq. (59) on a set of coordinate spheres defined by . On each of these spheres, we use spin-weighted spherical harmonics of spin weight in order to decompose the resulting wave signal into multipoles


In all our simulations, the wave signal is dominated by the , mode whose angular dependence is given by


The NP scalar is related to the gravitational wave strain via


It is convenient to decompose the two GW polarizations into multipoles in analogy to Eq. (64)


These multipoles are related to those of the NP scalar by


Note that the final result is not fully gauge-invariant and contains an unknown amount of systematic error. The reasons are two-fold: First, we did not choose a proper Bondi null tetrad on our extraction spheres, and, second, the extraction spheres have finite radius, thus are neglecting non-linear back-scattering effects of the gravitational field in the wave zone. However, since our coordinate frame asymptotically approaches the Minkowski spacetime, both errors can be minimized by performing extrapolation to future null infinity , using a set of extraction spheres at finite radii. Unfortunately, even if the extrapolation is accurate, an uncertain amount of residual error may remain. In Sec. IV.3, we check how well the extracted waves approximate their asymptotic shape and magnitude.

iii.4 Cauchy-Characteristic Extraction

To circumvent the problem of finite-radius extraction and to eliminate this systematic error, we apply the technique of CCE Bishop et al. (1999a); Babiuc et al. (2005); Winicour (2005); Reisswig et al. (2009b); Reisswig et al. (2010); Babiuc et al. (2010) to obtain the NP scalar as discussed in the previous section2, in this case directly evaluated at future null infinity . The CCE technique couples an exterior characteristic evolution of the full Einstein equations to the interior strong-field 3+1 Cauchy evolution of the spacetime.

Characteristic evolutions are based on null-hypersurface foliations of spacetime and have the advantage of allowing for a compactification of the radial coordinate component, thus allowing to include future null infinity on the computational grid Winicour (2005). Unfortunately, the characteristic formulation gives rise to the formation of caustics, i.e., the null rays on which the coordinate system is based can intersect in strong-field regions, leading to coordinate singularities. The scheme is therefore not well suited for the evolution of the actual GW source. Characteristic evolutions, on the other hand, are well adapted to the far-field region of spacetime and can efficiently evolve the metric fields out to where it is possible to obtain (and, hence, ) in a mathematically unambiguous and gauge-invariant way Bishop et al. (1997); Babiuc et al. (2005, 2009); Bishop and Deshingkar (2003).

We therefore proceed as follows: we evolve the interior region containing the collapsing matter with the standard Cauchy formulation as described in Sec. II.2 and II.3. During the Cauchy evolution, we store the 3-metric components including lapse and shift on coordinate spheres with fixed radius defining the world-tube .

This world-tube forms the inner boundary for the subsequent characteristic evolution of the Einstein equations. The full 4-metric can be reconstructed from the stored 3-metric components, the lapse and the shift at the inner boundary. Upon construction of proper initial data on an initial null hypersurface, which here we simply assume to be conformally flat, we have then fully specified any data necessary to evolve the fields out to . More details on the exact mathematical procedure can be found in Bishop et al. (1999a); Winicour (2005). The characteristic field equations are solved numerically using the PITT null evolution code Bishop et al. (1997). The numerical implementation of CCE including results from binary black hole mergers is discussed in Reisswig et al. (2009b); Reisswig et al. (2010); Babiuc et al. (2010). For the characteristic computational grid, we use points in the radial direction and points in each angular direction for the two stereographic patches covering the sphere. The characteristic timestep size equals that of the Cauchy evolution.

After each characteristic timestep, the NP scalar is evaluated directly from the metric at and transformed to the desired Bondi gauge Babiuc et al. (2009). Thus, the CCE method is free of gauge and near-zone effects and represents the most rigorous extraction technique. However, there is still some remaining systematic error that is due to the presence of matter at the world-tube locations. Since the current set of characteristic equations does not take into account any form of matter contribution, a non-zero stress-energy tensor introduces an unknown error. We therefore have to perform checks of the dependence of the waveforms on the world-tube locations. In principle, it is possible to also incorporate matter on the characteristic side Bishop et al. (1999b), which we leave to future work.

We note that CCE does not remove the artifical outer grid boundary from the Cauchy evolution. Thus, inconsistencies arising from this boundary can, in principle, still influence the interior domain. It is possible to circumvent this problem by enlarging the computational domain so that the outer boundary is causally disconnected from the world-tube locations (see, e.g., Pollney et al. (2009a) in the context of binary black holes). In simulations of core collapse, however, this is currently not computationally feasible, but experiments with varied outer boundary locations have shown that boundary effects are negligible for our current choice of boundary location.

Finally, we point out that inconsistencies in the characteristic and Cauchy initial data may lead to a loss of some non-linear effects. Even though we expect these problems to be very small (see Reisswig et al. (2010)). These and the outer boundary issues highlighted in the above can be fully accounted for only by employing Cauchy-characteristic matching (CCM) (e.g. Winicour (2005)). This technique uses the characteristic evolution as a generator for Cauchy boundary data at the world-tube, i.e. the world-tube becomes a two-way boundary between Cauchy and characteristic evolution. In practice, CCM has not yet been successfully implemented.

iii.5 Remarks on Integration and Physical Units

The NP scalar must be integrated twice in time to yield the strain , which introduces an artificial “memory” Berti et al. (2007), i.e., a non-linear drift of the signal so that the wavetrain deviates from an oscillation about zero. This behavior cannot be explained by the two unknown integration constants resulting, at most, in a linear drift.

As suggested in Reisswig and Pollney (2010), this non-linear drift is a consequence of random-walk-like behavior induced by numerical noise. In the present work, we make use of methods that are strictly valid only in pure vacuum and at our extraction spheres the average matter densities are non-zero (see Table 3). This systematic error can lead to an additional artificial low-frequency drift. In order to eliminate this effect, we use fixed-frequency integration (FFI) as proposed in Reisswig and Pollney (2010). The NP variable is Fourier transformed, the resulting spectrum is divided by for frequencies and divided by otherwise. An inverse Fourier transform then yields the strain essentially free of spurious drifts and oscillations, given a proper choice of .

Finally, we need to address the question of units. The gravitational wave strain is by construction dimensionless. For comparison of waveforms at different extraction radii, it is convenient to compensate for the fall-off of the strain, where is the distance from the observer to the source, and to work with and instead. In most of the following, we convert from code units, which are in , to cgs units when stating and plotting numerical results. The conversion factors we use are for length, and for time. For simplicity, we state the radii of GW extraction spheres and world-tube radii in code units. These and their corresponding cgs values are listed in Table 3.

Iv Results

In this section, we compare the most reliable extraction method that contains the least amount of systematic errors, CCE, with the various other curvature-based extraction methods, i.e., RWZM and NP extraction (both at finite radii). We also perform a comparison with the quadrupole approximation which has been employed in virtually all core collapse simulations to date.

This section is structured as follows. First, we review briefly the morphology of the gravitational waveforms expected from rotating core collapse and bounce.

Second, we elaborate in more detail on the method with which we obtain the gravitational strain from the quantities measured during the simulation. This is important, since the derived strain typically contains severe non-linear drifts making a proper analysis largely impossible without significant preprocessing.

Third, we assess the accuracy of each individual method, i.e., we analyze the radial dependence of the NP and RWZM extraction methods, since they are strictly valid only in an asymptotic frame at an infinite distance from the source where any contributions from the stress-energy tensor vanish. Since the matter densities are non-zero at the CCE world-tube locations, we also analyze the radial dependence of the waveforms extracted via CCE.

Fourth, we compare the results obtained via NP and RWZM extraction, and the approximate QF with results obtained via CCE.

Finally, we perform a convergence check on the computed waveforms by using a set of three different resolutions.

iv.1 Morphology of Rotating Core Collapse Waveforms

The core collapse models considered in this work remain nearly axisymmetric during collapse and emit GWs predominantly via the even-parity spherical harmonic mode. This mode has a maximum on the equator and, hence, we plot all waveforms as seen by an observer in the equatorial plane. We write where is the spin-weighted spherical harmonic with spin . Note that the mode is axisymmetric and thus, the equatorial strain has no -dependence.

We convert to cgs units by using the transformation as discussed in Sec. III.5, and we align the maxima of the waveforms such that they occur at , corresponding roughly to the time of core bounce in each model.

The waveforms of the three models are shown in Figs. 1 (model A1B3G3), 2 (model A1B3G5), and 3 (model A3B3G3). All models exhibit a very similar behavior. Prior to core bounce , the matter undergoes an aspherical accelerated collapse. Due to this aspherical acceleration, the GW signal is monotonically rising until it peaks when the contracting inner core is drastically decelerated. This deceleration is caused by the sudden stiffening of the EOS as a result of nuclear repulsive forces which emerge when nuclear densities are reached. During this deceleration, the GW signal becomes rapidly negative, reaching its second peak (the “bounce peak”) roughly when the core rebounds. Subsequently, the inner core undergoes a relaxation phase (ring-down) in which it loses its remaining pulsation energy by launching secondary shocks. This results in an oscillatory GW signal that decreases in amplitude as the core approaches its final equilibrium.

While the overall morphology of the GWs emitted by the three models is the same, there are subtle differences that are worth commenting on. Models A1B3G3 and A3B3G3 produce so-called type-I signals Zwerger and Müller (1997) with a single pronounced major peak at core bounce. Since model A3B3G3 is more rapidly spinning, its inner core is more deformed, and hence produces a stronger GW signal at core bounce than model A1B3G3. Model A1B3G5 has a very small inner core at bounce and produces a type-III signal that is characterized by a much less pronounced bounce peak and generally low-amplitude GW emission. Note that type-II signals, characterized by multiple wide and pronounced bounce peaks seen in early work Zwerger and Müller (1997); Mönchmeyer et al. (1991); Kotake et al. (2003); Yamada et al. (1999); Ott et al. (2004) have been demonstrated to disappear in simulations using general relativity and a proper electron-capture treatment Ott et al. (2007a); Dimmelmeier et al. (2008).

iv.2 Computing the Strain

We first consider the computation of the strain from the RWZM formalism. Since our models emit GWs predominantly in the even mode , the computation of the strain from the even- and odd-parity RWZM master functions reduces to Eq. (50) so that no time integral is necessary to obtain . However, we still notice an unphysical drift in the waves. Since the RWZM master functions are computed at a finite distance from the source, we have the following systematic errors (summarized as the “finite-radius error”): (i) a non-vanishing matter density at the extraction spheres, (ii) near-zone effects and (iii) gauge ambiguities. The latter error arises as a result of deviations from the Bondi gauge (see Lehner and Moreschi (2007) for an improvement). The artificial drift is part of the finite-radius error, since it is becoming less pronounced with increasing extraction radius.

In order to reduce the contribution of these artificial low-frequency components, we first transform to the Fourier domain, multiply by in order to take the first time derivative, and then apply fixed FFI Reisswig and Pollney (2010) to obtain