The influence of model parameters on the prediction of gravitational wave signals from stellar core collapse
Key Words.:Gravitational waves – (Stars:) supernovae: general – Hydrodynamics – Neutrinos – Stars: rotation – Stars: neutron
We present a gravitational wave (GW) analysis of an extensive series of three-dimensional magnetohydrodynamical core-collapse simulations. Our 25 models are based on a 15 progenitor stemming from (i) stellar evolution calculations, (ii) a spherically symmetric effective general relativistic potential, either the Lattimer-Swesty (with three possible compressibilities) or the Shen equation of state for hot, dense matter, and (iii) a neutrino parametrisation scheme that is accurate until about 5ms postbounce. For three representative models, we also included long-term neutrino physics by means of a leakage scheme, which is based on partial implementation of the isotropic diffusion source approximation (IDSA). We systematically investigated the effects of the equation of state, the initial rotation rate, and both the toroidal and the poloidal magnetic fields on the GW signature. We stress the importance of including of postbounce neutrino physics, since it quantitatively alters the GW signature. Slowly rotating models, or those that do not rotate at all, show GW emission caused by prompt and proto-neutron star (PNS) convection. Moreover, the signal stemming from prompt convection allows for the distinction between the two different nuclear equations of state indirectly by different properties of the fluid instabilities. For simulations with moderate or even fast rotation rates, we only find the axisymmetric type I wave signature at core bounce. In line with recent results, we could confirm that the maximum GW amplitude scales roughly linearly with the ratio of rotational to gravitational energy at core bounce below a threshold value of about . We point out that models set up with an initial central angular velocity of rad or faster show nonaxisymmetric narrow-band GW radiation during the postbounce phase. This emission process is caused by a low dynamical instability. Apart from these two points, we show that it is generally very difficult to discern the effects of the individual features of the input physics in a GW signal from a rotating core-collapse supernova that can be attributed unambiguously to a specific model. Weak magnetic fields do not notably influence the dynamical evolution of the core and thus the GW emission. However, for strong initial poloidal magnetic fields (G), the combined action of flux-freezing and field winding leads to conditions where the ratio of magnetic field pressure to matter pressure reaches about unity which leads to the onset of a jet-like supernova explosion. The collimated bipolar out-stream of matter is then reflected in the emission of a type IV GW signal. In contradiction to axisymmetric simulations, we find evidence that nonaxisymmetric fluid modes can counteract or even suppress jet formation for models with strong initial toroidal magnetic fields. The results of models with continued neutrino emission show that including of the deleptonisation during the postbounce phase is an indispensable issue for the quantitative prediction of GWs from core-collapse supernovae, because it can alter the GW amplitude up to a factor of 10 compared to a pure hydrodynamical treatment. Our collapse simulations indicate that corresponding events in our Galaxy would be detectable either by LIGO, if the source is rotating, or at least by the advanced LIGO detector, if it is not or only slowly rotating.
Gravitational wave (GW) astronomy may soon become a reality and will allow humankind to address questions about many different astrophysical objects that are hidden from the electromagnetic detection. Within the past few years, the first generation of the ground-based GW detectors LIGO (USA), VIRGO (Italy), GEO600 (Germany), and TAMA (Japan) have got very close to or even reached design sensitivity and collected partially coincident data, as discussed by Abbott et al. (2009). Lately, the two 4 km LIGO detectors were upgraded to sensitivities increased by a factor of 2-3 (enhanced LIGO, see Adhikari et al. (2006)) and resumed observations in 2009. Upgrades of the three LIGO interferometers and VIRGO are expected to be completed by 2014 and will increase the observable volume by a factor of 1000. Operating at such a high level of precision, GW detectors are sensitive to many different sources, such as compact binary coalescence from black-holes (BH) and neutron stars up to distances of several 100 megaparsecs, but might also very likely produce the first detections of black-hole and neutron star mergers, core-collapse supernovae, and neutron star normal mode oscillations, as recently reviewed in Sathyaprakash & Schutz (2009).
The broad majority of GW sources can be subdivided into two classes: (i) mathematically well-posed events such as BH-BH coalescence and (ii) scenarios that involve matter. The first category of sources can be modelled accurately, the waveform can be calculated with high precision, and thus a matched-filter analysis can be applied. On the other hand, the latter kind of GW sources can only be modelled imperfectly, as matter effects carry large physics uncertainties and open up a huge parameter space for initial conditions. Moreover, even if there were ‘perfect’ models, there is still a turbulent, stochastic element in some GW emission mechanisms that makes it strictly impossible to compute templates. As a result, more general, un-modelled burst analysis techniques are used, as summarised in Abbott et al. (2009) and references therein. One of the scenarios where burst data analysis must be applied is the stellar core collapse, where in addition to the matter effects even the fundamental explosion mechanism is not fully settled. It is still unclear which processes may convert some of the released gravitational binding energy of order erg into the typically observed kinetic and internal energy of the ejecta of erg = 1 Bethe [B] needed to blow off the stellar envelope. While most state-of-the art simulations investigate a mechanism based on neutrino heating in combination with hydrodynamical instabilities, such as Marek & Janka (2009), others explore the alternative magneto-rotational mechanism (Leblanc & Wilson 1970; Kotake et al. 2004a; Burrows et al. 2007; Mikami et al. 2008; Takiwaki et al. 2009), or the acoustic mechanism which is driven by strong core g-modes (Burrows et al. 2006, 2007). For a recent review see Janka et al. (2007). Takahara & Sato (1988), Gentile et al. (1993), and lately also Nakazato et al. (2008) and Sagert et al. (2009) reported that a QCD phase transition may power a secondary shock wave which triggers a successful hydrodynamical explosion.
Beside neutrinos, which have already been observed in the context of stellar core collapse of SN1987A (Hirata et al. 1988), GWs could provide access to the electromagnetically hidden compact inner core of some such cataclysmic events. Since a core-collapse supernova is expected to show aspherical features (Leonard et al. 2006), there is reasonable hope that a tiny amount of the released binding energy will be emitted as GWs which could then provide us with valuable information about the angular momentum distribution (Dimmelmeier et al. 2008) and the baryonic equation of state (EoS) (Marek et al. 2009), both of which are uncertain. Furthermore, they might help to constrain theoretically predicted SN mechanisms (Ott 2009b). Given the actual sensitivities of today’s operating ground-based GW observatories (Whitcomb 2008) combined with knowledge of waveforms from state-of-the art modelling, GWs from a Galactic core-collapse supernova may be considered to be within the detector limits. The collapsing iron core and the subsequently newly formed PNS can be subject to a whole variety of asymmetric hydrodynamical and nonaxisymmetric instabilities which give rise to GW emission.
In this context probably most attention during the past three decades has been paid to rotational core collapse and bounce dynamics. In this particular phase, the core spins up due to angular momentum conservation, resulting in an oblate and time-dependent deformation that leads to strong GW emission (Müller 1982; Mönchmeyer et al. 1991; Janka & Müller 1996; Zwerger & Müller 1997; Rampp et al. 1998; Kotake et al. 2003, 2004b; Obergaulinger et al. 2006; Ott et al. 2004, 2007a, 2007b; Dimmelmeier et al. 2002, 2007, 2008; Scheidegger et al. 2008). The steady improvement of the models (e.g. the inclusion of GR, a micro-physical EoS, treatment of neutrino physics) recently led to a theoretically well understood single and generic, so-called type I wave form which is characterised by a large negative peak at core bounce, followed by ring-down oscillations that damp quickly (Ott et al. (2007a, b); Dimmelmeier et al. (2007, 2008); with detailed references therein). Despite the reduction to a single wave form, the combined information of the GW amplitude and the location of the narrow peak of the GW spectral energy density in frequency space contains information that makes it still possible to constrain progenitor- and postbounce rotation, but can barely distinguish between different finite-temperature EoS (see e.g. Dimmelmeier et al. (2008), who used the EoS of Lattimer & Swesty (1991) and Shen et al. (1998b)).
Gravitational waves from magneto-rotational collapse were considered in detail by Kotake et al. (2004b); Shibata et al. (2006) and Obergaulinger et al. (2006). As main differences compared to simulations which do not include magnetic fields it was found by these groups that only in the very special case of precollapse fields as strong as G the overall dynamics can be influenced. As Heger et al. (2005) have argued, such strong fields are unlikely to occur in standard core-collapse supernova progenitors. The GW amplitude is then affected by (i) time-dependent magnetic fields, which contribute considerably to the overall energy density, and (ii) by bipolar magnetohydrodynamic (MHD) jet outflows which give rise to a so-called type IV signal with memory (Obergaulinger et al. 2006). Physically, such a memory effect in the GW signal arises from the temporal history of asymmetric matter outflow, leaving behind a constant offset in the amplitude (Thorne 1989).
Recently it has been argued through numerical simulations of equilibrium neutron star models or full core-collapse simulations that PNSs with a high degree of differential rotation can be subject to nonaxiymmetric rotational instabilities at low values (= , ratio of rotational to gravitational energy), leading to strong narrow-band GW emission (Saijo et al. 2003; Watts et al. 2005; Ott et al. 2005; Saijo & Yoshida 2006; Ou & Tohline 2006; Ott et al. 2007b; Cerdá-Durán et al. 2007; Ott et al. 2007a; Scheidegger et al. 2008). However little is known about the true nature of the instability at present. Previous work has failed to establish an analytical instability criterion, and the dependence of the instability on PNS rotation rate and degree of differential rotation is still unclear, as it was pointed out by Ott (2009b).
Current state-of-the-art stellar evolution calculations (Heger et al. 2005) tell us that iron cores of stars generally lose most of their angular momentum during their evolution due to magnetic torques. Therefore, they cannot become subject to strong rotationally-induced aspherities. However, anisotropic neutrino emission, convection and standing accretion shock instability (SASI, see e.g. Blondin et al. (2003)) - driven deviations from spherical symmetry which can lead to the emission of GWs of sizable amplitudes are likely to occur inside the PNS and the post-shock gain region and last for probably hundreds of ms (Müller & Janka 1997; Müller et al. 2004; Kotake et al. 2007; Marek et al. 2009; Murphy et al. 2009; Kotake et al. 2009a, b). Although there is qualitative consensus among the different core-collapse supernova groups about the aforementioned features being emitters of stochastic broad-band signals, their detailed quantitative character still remains quite uncertain, since self-consistent 3D simulations with proper long-term neutrino transport were not carried out so far, but would in principle be required. The most recent, very elaborate 2D simulations in that context were performed by Marek et al. (2009). Their numerical setup includes long-term multi-flavor neutrino transport, an effective relativistic potential and two different EoS. The first 100ms of after bounce they observed GWs from early prompt postbounce convection, peaking somewhat around or below 100Hz, depending on the employed EoS. After this early episode, the GW emission in their models is dominated by a growing negative amplitude related to anisotropic neutrino emission at frequencies Hz, while the GW signal associated with non-radial mass motions stems from density regimes gcm and peaks in the frequency range of Hz, highly sensitive to the nuclear EoS.
In the core-collapse scenario, PNS pulsations can provide another mechanism for GW emission. Ott et al. (2006a) pointed out that in the context of the acoustic mechanism (Burrows et al. 2006) excited core g-mode oscillations might emit very strong GWs. The GW signal is due to the nonlinear quadrupole components of the pulsations that, at least initially, are of g-mode character (Ott et al. 2006a).
In this paper, we present the gravitational wave analysis of a comprehensive set of three-dimensional MHD core-collapse supernova simulations in order to investigate the dependencies of the resulting GW signal on the progenitor’s initial conditions. Our calculations encompass presupernova models from stellar evolution calculations, a finite-temperature nuclear EoS and a computationally efficient treatment of the deleptonisation and neutrino emission during the collapse. General relativistic corrections to the spherically symmetric Newtonian gravitational potential are taken into account. Moreover, several models incorporate long-term neutrino physics by means of a leakage scheme. As for the progenitor configuration, we systematically consider not only variations in the precollapse rotation rate, but also in the nuclear EoS and the magnetic field topology. This study extends the work of Scheidegger et al. (2008), who investigated only two models with similar input physics. In this way, we carried out the so far largest parameter study of 3D MHD stellar collapse with respect to the prediction of GWs.
The structure of the paper is as follows. In section 2 we briefly describe the numerical methods, input physics and progenitor configurations applied in our core-collapse simulations. Furthermore, we describe how we extract GWs from our model set. In sec. 3 we discuss results of 25 three-dimensional MHD simulations with respect to the GW signal from convection, rotational core bounce, nonaxisymmetric rotational instabilities and very strong magnetic fields. In sec. 4, we summarise and present conclusions.
2 Numerical Methods
2.1 The 3D MHD code and its input physics
The algorithm used to solve the time-dependent, Newtonian MHD equations in our 3D simulations is based on a simple and fast cosmological MHD code of Pen et al. (2003); Liebendörfer et al. (2006), which has been parallelised with a hybrid combination of MPI and openMP and improved and adapted to the requirements of core-collapse supernova simulations (Käppeli et al. 2009). The 3D computational domain consists of a central cube of km volume, treated in equidistant Cartesian coordinates with a grid spacing of 1km. It is, as explained in detail in Scheidegger et al. (2008), embedded in a larger spherically symmetric computational domain that is treated by the time-implicit hydrodynamics code ’Agile’ (Liebendörfer et al. 2002). With this setup, the code scales nicely to at least 8000 parallel processes (Käppeli et al. 2009). The ideal MHD equations read
expressing the conservation of mass, momentum, energy and magnetic flux, respectively. stands for mass density, v is the velocity vector and the total energy (the sum of internal, kinetic and magnetic energy). The magnetic field is given by and the total pressure by (the sum of gas and magnetic pressure). The MHD equations are evolved from a condition
We employ the constrained transport method (Evans & Hawley 1988) to guarantee the divergence-free time evolution of the magnetic field. The right hand side of Eq. 2 and Eq. 3 take into account the effect of gravitational forces on the magnetohydrodynamical variables. The gravitational potential obeys the Poisson equation
We only implement the monopole term of the gravitational potential by a spherically symmetric mass integration that includes general relativistic corrections (Marek et al. 2006). In this approach, the Newtonian gravitational potential is replaced by an effective potential (‘case A’ potential, see Marek et al. (2006))
with being the gas pressure, the neutrino pressure and the internal energy density. The effective mass is given by
where is the neutrino energy. The metric function is given by
where is the radial fluid velocity. Note that Müller et al. (2008) recently proposed an effective relativistic gravitational potential for rapidly rotating configurations. However, since it would have added another degree of freedom to our parameter set, we did not consider it in this study.
The system of the MHD equations must be closed by a finite-temperature EoS. In core-collapse supernova simulations, an EoS has to handle several different regimes. For temperatures below 0.5MeV, the presence of nuclei and time dependent nuclear processes dominate the internal energy evolution. For simplicity, an ideal gas of Si-nuclei is assumed, which, for our GW study, is sufficiently accurate to describe the low baryonic pressure contribution at low densities in the outer core. At higher temperatures, where matter is in nuclear statistical equilibrium, we employ two alternative EoS: The Lattimer & Swesty EoS (Lattimer & Swesty (1991), LS EoS) and the one by Shen et al. (1998a). The first EoS is based on a phenomenological compressible liquid drop model; it also includes surface effects as well as electron-positron and photon contributions. The LS EoS assumes a nuclear symmetry energy of 29.3MeV and we will perform simulations with three choices of the nuclear compressibility modulus K (180,220,375MeV) as provided by Lattimer & Swesty (1991). Since variations in K affect the stiffness of the nuclear component of the EoS, this enables us to probe the effects of variations in stiffness while keeping the general EoS model fixed. The pure baryon EoS from Shen et al. (1998a) has a compressibility of K=281MeV and a symmetry energy of 36.9MeV. It is based on relativistic mean field theory and the Thomas-Fermi approximation. For matter in non-NSE (TMeV), the Shen EoS is coupled to the baryonic and electron-positron EoS given in Timmes & Arnett (1999) and Timmes & Swesty (2000). It employs an ideal gas for the nuclei and additionally includes contributions from ion-ion-correlations and photons.
2.2 The treatment of neutrino physics
The treatment of neutrino physics is an essential ingredient of core-collapse supernova simulations (Mezzacappa 2005). Multidimensional core-collapse supernova simulations therefore must rely on more or less severe approximations of the neutrino physics. In state-of-the-art 2D simulations, one way in which Boltzmann transport is approximated is the so-called ‘ray-by-ray plus scheme’ (Buras et al. 2006). It solves the full transport in separate 1D angular segments, where the neighbouring rays are coupled. Other groups rely on multi-group flux limited diffusion (MGFLD, Ott et al. (2008); Swesty & Myra (2009)). MGFLD treats all neutrinos in seperate energy groups, drops the momentum space angular dependence of the radiation field and evolves the zeroth moment of the specific intensity instead of the specific intensity itself. Ott et al. (2008) also compared MGFLD with angle-dependent (i.e. partly Boltzmann) transport in 2D. In three dimensions, simulations have been performed using ‘grey’ flux-limited diffusion (Fryer & Warren 2004), which oversimplyfies the important neutrino spectrum. It is important to resolve the neutrino spectrum, since the charge-current interaction rates go with the square of the neutrino energy.
In our 3D MHD simulations, we apply a parametrised deleptonisation scheme (Liebendörfer 2005). Detailed spherically-symmetric collapse calculations with Boltzmann neutrino transport show that the electron fraction in different layers in the homologously collapsing core follow a similar deleptonisation trajectory with density. The local electron fraction of a fluid element can be parametrised as a function of density . For our 3D simulations, we apply tabulated profiles which were obtained from detailed general relativistic, spherically symmetric three-flavour Boltzmann neutrino transport. Exemplary profiles are shown in Fig. 1. The offset between the profiles from the LS- and the Shen EoS is consistent with the different asymmetry energies of the two EoS. This in turn is reflected in different neutrino reaction rates and thus a different at a given density. These effects have been discussed in Sumiyoshi et al. (2008) and Fischer et al. (2009) for massive progenitor stars in the range of 40-50. However, this parametrisation scheme is only valid until a few milliseconds after bounce, since it cannot account for the neutronisation burst, as explained in Liebendörfer (2005) and Scheidegger et al. (2008).
In order to give an estimate of how long the above approximation is able to predict quantitatively reliable postbounce (pb) GW signals for times ms after bounce, we carried out three representative simulations (R1E1CA, R3E1AC and R4E1FC, see Tab. 1) that incorporate a postbounce treatment for neutrino transport and deleptonisation. The scheme we apply is based on a partial implementation of the isotropic diffusion source approximation (IDSA; Liebendörfer et al. (2009)). The IDSA splits the distribution function of the neutrinos into two components, a trapped component and a streaming component , representing neutrinos of a given species and energy which find the local zone opaque or transparent, respectively. The total distribution function is the sum of the two components, . The two components are evolved using separate numerical techniques, coupled by a diffusion source term. The trapped component transports neutrinos by diffusion to adjacent fluid elements, while in this paper the streaming component is discarded (), implying that these neutrinos are lost from the simulation immediately. The part of the IDSA for trapped neutrinos is implemented in three dimensions, including both electron neutrinos and electron anti-neutrinos (Whitehouse & Liebendörfer 2010, in preparation). The use of this partial implementation of the IDSA enables us to capture the neutrino burst and to continue our simulations into the postbounce regime, as shown in Fig. 2.
The leakage scheme is switched on at core bounce, when the central density of the core reaches its first maximum value. At present, we neglect the emission of and neutrinos and their antineutrinos, which in principle play an important role in the cooling of the PNS. However, since our leakage schemes neglects any absorption of transported neutrinos, it already overestimates the cooling of the PNS, even without the treatment of the and neutrinos. This is shown in 3.
2.3 Initial model configurations
We construct the initial conditions of our simulations by a parametric approach. All our models are launched from a progenitor star from the stellar evolutionary calculations of Woosley & Weaver (1995). Angular momentum was added to the presupernova model according to a shell-type rotation law (Eriguchi & Müller 1985)
where we define as spherical radius, and (x,y,z) as Cartesian coordinates. The constant A is the degree of differential rotation that controls the steepness of the angular velocity profile, the initial central rotation rate, and is the distance from the origin. In our whole model set we choose km as degree of differential rotation.
In order to guarantee a divergence-free initial state, the initial magnetic field configuration was set up employing its definition via the vector potential . The components of we chose are
In order to mimic a dipole-like field, we scale the vector potential with density according to
Finally, the magnetic field is derived from this vector potential
The initial toroidal- and poloidal components of the magnetic field are specified at a reference density of gcm according to Heger et al. (2005).
We compute a total of 25 models, changing the combination of total angular momentum, the EoS, and toroidal- and poloidal magnetic fields. The model parameters are summarised in table 1. The models are named after the combinations of initial central rotation rate, the EoS, and toroidal- and poloidal magnetic fields. The first two letters of the model name represent the initial central rotation rate [rads] according to
the second two letters stand for the applied EoS
|(K = 180MeV)||(K = 220MeV)||(K =325 MeV)|
, while the last two letters assign the order of magnitude of the toroidal-, and poloidal field strength [G] according to
The final subscript signs simulations which were carried out with the leakage scheme. Some of the values which we adopt as rotation rate and magnetic fields correspond to the values suggested in Heger et al. (2005). However, we point out that some of the initial rotation rates (rads) and magnetic fields (G) are larger and stronger compared to current predictions from stellar evolution calculations. However, we computed these models in order to cover a wide parameter space for our three-dimensional models with neutrino transport approximations.
2.4 Gravitational Wave extraction
The two independent polarisations of the dimensionless gravitational wave field in the transverse traceless gauge are given by
The spatial indices run from 1 to 3 (Misner et al. 1973). is the distance from the source to the observer and the unit polarisation tensors and in spherical coordinates are represented as
In the slow-motion limit (Misner et al. 1973; Finn & Evans 1990) the amplitudes and are linear combinations of the second time derivative of the transverse traceless mass quadrupole tensor . In Cartesian coordinates, the quadrupole tensor is expressed as
For convenience we evaluate below the GW amplitudes along the polar axis (, denoted as subscript I)
and in the equatorial plane (, , denoted as II)
Since a direct evaluation of Eq. (17) is numerically problematic, as discussed in Finn & Evans (1990), we apply alternative reformulations of the standard quadrupole formula, in which one or both time derivatives are replaced by hydrodynamic variables using the continuity and momentum equations. In the first moment of momentum density formula (see Finn & Evans (1990)), the first time derivative of the quadrupole moment yields
In the stress formulation Blanchet et al. (1990), is given by
where is the gravitational potential. We apply both ways of computing the GW signal and compare them by computing the overlap (Giacomazzo et al. 2009) of the resulting waveforms via
where [Hz] is the power spectral density of the strain noise from a given detector, e.g. LIGO, which were kindly provided by Shoemaker (2007, private communication) The results show good agreement with numerical deviations of a few percent, as displayed in Table 2. There, we computed the overlap of GW trains from the representative models R1E1CA, R3E1CA and R4E1FC, extracted by using Eq. 22 and Eq. 23. We assume optimal orientation of detector and source as well as the GW to be emitted along the polar axis.
Below, we generally apply Eq. 22 if not stated otherwise. Given our spherically-symmetric effective GR approach to the solution of the Poisson equation, this expression is physically best motivated, since it does not depend on spatial derivatives of the gravitational potential.
We also point out that the Newtonian quadrupole formalism to extract the gravitational radiation is not gauge invariant and only valid in the Newtonian slow-motion limit. However, it was shown by Shibata & Sekiguchi (2003) that the above method seems to be sufficiently accurate compared to more elaborate techniques, as it preserves phase while being off in amplitude by in neutron star pulsations. The energy carried away by gravitational radiation can be calculated by the following expression:
where . The equivalent frequency integral yields
where is the Fourier transform of the quadrupole amplitude .
In order to calculate the contribution to the GW signal due to magnetic stresses, we generalised Eq. (23), taking into account contributions from the magnetic field. Following the derivations of Kotake et al. (2004a) and Obergaulinger et al. (2006), the Cartesian quadrupole gravitational wave amplitude in the MHD case yields
where , and . Below, the GW source is assumed to be located at the Galactic centre at kpc. We estimate the signal-to-noise ratios (SNR) for optimal filtering searches according to Flanagan & Hughes (1998):
where is the Fourier transform of the gravitational wave amplitude, and is the power spectral density of strain noise in the detector. We assume optimal orientation of detector and source. and the Fourier spectra were normalised according to Parseval’s theorem.
In the subsequent four subsections we will discuss the GW signature of our 25 models. While presenting the resulting GW patterns, we will pay special attention to possible imprints of different finite temperature EoS, rotation rates, nonaxisymmetric instabilities, magnetic fields and a postbounce leakage scheme on the predicted 3D GW signals. The models’ initial conditions and similar relevant quantities are summarised in Table 1, whilst the GW data is listed in Table 3.
3.1 Non- or slowly rotating core collapse
Non- and slowly rotating progenitors (rads in our model set) all undergo quasi-spherically symmetric core collapse. As the emission of GWs intrinsically depends on dynamical processes that deviate from spherical symmetry, the collapse phase in our models (that neglect inhomogenities in the progenitor star) therefore does not provide any kind of signal, as shown in Fig. 3 for . However, subsequent pressure-dominated core bounce, where the collapse is halted due to the stiffening of the EoS at nuclear density gcm, launches a shock wave that plows through the infalling material, leaving behind a negative entropy gradient that induces so-called ‘prompt’ convective activity (e.g. Müller et al. (2004); Buras et al. (2006); Dimmelmeier et al. (2008); Marek et al. (2009); Ott (2009b)). The GW burst which is accompanied by such aspherities starts several ms after bounce when convective overturn starts to be effective. The criterion for convective instability (the ‘Ledoux condition’) is generally expressed as (Landau & Lifshitz 1959; Wilson & Mayle 1988)
where , , and are electron fraction, density, stellar radius and entropy per baryon respectively. Thermodynamic consistency requires , which implies that a negative entropy gradient always acts in a destabilising manner. Additionally, the neutronisation burst, occurring some ms after bounce, causes a negative lepton gradient at the edge of the PNS which further drives convection (cf. Fig. 13 of Sumiyoshi et al. (2004), Buras et al. (2006); Dessart et al. (2006)). Note that convection is be weakened in the rotational plane by positive specific angular momentum gradients in rotating cores (Endal & Sofia 1978).
Models without deleptonisation in the postbounce phase: Effects of the EoS and magnetic fields on the GW signature
During the early postbounce stage (ms), prompt convective motion is predominantly driven by a negative entropy gradient, as pointed out e.g. by Scheidegger et al. (2009) (see also Dimmelmeier et al. (2008); Marek et al. (2009); Ott (2009b)). The GW burst to be associated with prompt convection sets in ms after bounce in models based on the LS EoS, generally ms before the same feature occurs in the corresponding simulations using the Shen EoS, as indicated in Fig. 3. The reason for this behaviour is that the shock wave in the ‘Shen’-models carries more energy compared to those in models using a LS EoS. Therefore, the shock wave stalls at slightly later times at larger radii (see Fig. 4), and the conditions for convective activity are delayed compared to the LS runs. Note that the convective overturn causes a smoothing of the negative entropy gradient. As a result, the GW amplitude quickly decays (ms after bounce) and is not revived during the later evolution of the models without deleptonisation in the postbounce phase, as displayed in Figs. 3 and 5. Simulations that incorporate the Shen EoS return up to a factor of 2 smaller maximum amplitudes compared to their counterparts, as can be deduced from Fig. 3 and Table 3.
We also find that the GWs from simulations which were carried out with the stiff E3 show no significant deviations from those computed with E1 (see Table 3 and Fig. 5). The minor deviations of the GWs are entirely due to the stochastic nature of convection. The waveform spectra from the LS models cover a broad frequency band ranging from Hz, as displayed in Fig. 6. The spectral peak of the Shen models is shifted somewhat to lower frequencies, covering a range from Hz (see Fig. 6). When comparing the energy which is emitted during the first 30ms after bounce from the LS models to that of the corresponding Shen models, we find the latter models emit less energy (see Table 3). This discrepancy is due to the lower emission at higher frequency Hz in the Shen models and the fact that the emitted energy is proportional to . Moreover we find that already slow rotation (rotation rate R1) leads to a deformation of the PNS. This can be quantified by considering e.g. a density cut at gcm, which is just inside the convectively unstable region. For the slowly rotating model R1E1CA, this point is located at a radial distance of 80km along the polar axis at 20ms after bounce. However, in the equatorial plane, the same position is reached one radial grid zone further out relative to the origin due to the action of centrifugal forces. The short time variation in the quadrupole due to rotation combines with that of prompt convection and together they lead to somewhat stronger GW emission in the slowly rotating case compared to the non-rotating model set. This effect is strongest for the GW amplitude A, which is the ‘axisymmetric’ (l=2, m=0) component of the wave field. Despite this feature, the frequency content of models which only differ in rotation rate, stays practically unaltered.
We found the key controlling factors that govern
the GW emission from prompt convection
to be the i) radial location
of the convectively unstable zones and
ii) the related characteristic dynamical
for which we use as rough estimate
Computing the SNRs we find that all the simulations discussed above lie just below the detector limits of LIGO if we assume them to be located at a Galactic distance of 10kpc. Their single-detector optimal-orientation SNR is just a little above unity. However, for a successful detection, at least a SNR of 7 to 8 is necessary. Note that the SNRs of models with different EoS do not differ much at this stage since all share a similar spectral energy distribution within the window of LIGO’s maximum sensitivity. The current detector sensitivity does not allow for the detection of the high frequency tail of the LS models. Note, however, that for planned future detectors such as the Advanced LIGO facility, things change dramatically. As a direct consequence, these new detectors would permit the distinction between the prompt convection GW signal from the LS and Shen EoS, since the full spectral information would be available. However, we find it impossible to discriminate between the different LS EoS variants. Hence, our simulations indicate that the GW signature depends more strongly on the asymmetry energy than the compressibility parameter of the EoS.
These results are partly different from those previously published. Recently, Marek et al. (2009) as well as Ott (2009b) reported to observe GW from prompt convection in state-of-the-art 2D simulations which were launched from similar initial conditions as ours, namely the same progenitor star and the soft variant of the LS EoS (K=180MeV). Whilst the extracted GW amplitudes from early prompt convection of Marek et al. (2009) (cf. their model M15LS-2D) are in rough agreement with our results, the spectrum of their wave train peaks at considerably lower frequencies, namely around about Hz. We suppose that this discrepancy is due mainly to different radial locations of the unstable regions and the consequently encompassed amount of overturning matter. Ott (2009b) computed two models with different resolution. While one of their models (s15WW95) in particular fits our results well in all characteristic GW features, namely the size of amplitudes, band of emission and the amount of emitted energy, the better resolved model (s15WW95HR) showed that convection is much weaker due to less seed perturbations and hence the GW signal and the total amount of emitted energy considerably lower. We recently also tested this issue with a better resolved model (cf. model R1 of Scheidegger et al. (2009), with a grid spacing of 0.6km). This better resolved model showed considerably smaller seed perturbations around , as grid alignment effects are better suppressed at core bounce; hence prompt convection then is much weaker and a smaller GW amplitude (%) is emitted. However, better numerical resolution also leads to less numerical dissipation in the system, which eases the dynamical effects that follow. Thus, Scheidegger et al. (2009) found for ms considerably stronger GW emission from early prompt convection compared to the 1km resolved models.
Model with deleptonisation in the postbounce phase
Due to the absence of accurate postbounce neutrino transport, it is unclear how reliably the models discussed in the previous subsection predict the GW signals from the early postbounce period. In order to investigate this question, we carried out one computationally expensive simulation, model R1E1CA, that includes the emission of neutrinos after bounce but neglects the neutrino heating, which becomes relevant at ms. The absence of the latter makes it impossible for this model to track the long-term postbounce neutrino-driven convection. However, in comparison with purely MHD models, the inclusion of neutrino leakage makes is possible to distinguish effects on the GW signal due to entropy and lepton-gradient driven convection.
A detailed comparison of model R1E1CA with its purely magnetohydrodynamical counterpart R1E1CA shows that both follow a similar dynamical behaviour until about 20ms after bounce. Asphericities leading to GW emission are predominantly driven by entropy- and not lepton-induced convection in this supernova stage. Consequently, the wave trains emitted within this early period fit each other qualitatively (cf. Fig. 3 and Fig. 8). However, we find some quantitative deviations: The GWs of R1E1CA reach lower maximum values (see Table 3), as the ‘Ledoux’ unstable region encompasses less mass (Fig. 9), and the presence of neutrino cooling leads to a more rapid smoothing of the entropy gradient compared to the models discussed in the previous subsection. However, since the overturning matter in the top layers of the PNS has the same radial position, densities and dynamical timescales as in model R1E1CA, these models have similar GW spectra, peaking between Hz (cf. Fig. 6 and Fig. 10). Therefore the physically simpler models still provide reasonably accurate GW predictions in frequency space until about 20ms after bounce, although the amplitudes are overestimated a few % (cf. Figs. 3 and 8). Model R1E1CA’s later postbounce evolution (ms) differs strongly compared to its purely hydrodynamical counterpart R1E1CA. A negative radial lepton gradient, caused by the neutronisation burst and subsequent deleptonisation, drives convection inside the lower layers of the PNS (Dessart et al. 2006) at a radial position of 10-30km and a density range of gcm and therefore causes now the entire GW emission. The cooling PNS contracts with time, which causes the convective zones to migrate towards smaller radii and shrink. However, we point out that our leakage scheme overestimates the neutrino cooling processes, as shown in Fig. 2. Hence, this mechanism proceeds too quickly for model R1E1CA (Ott 2009, private communication). The PNS convection exhibits GW emission of roughly cm amplitude, as can be seen in Fig. 8 (ms). The corresponding spectral distribution is shown in Fig. 10. A broad peak rises between Hz and reflects the dynamical timescale of the violent overturn activity of a millisecond scale inside the PNS. The model’s SNR for LIGO at 10kpc is again around unity. The high frequency tail of the spectrum (Hz), which is present due to PNS convection, cannot contribute to the SNR as it lies below the current detector sensitivity. Our computed GW strains for PNS convection agree roughly in amplitude with the ones found in Ott (2009b) for axisymmetric MGFLD models. However, the amount of released energy emitted is found to be about one order of magnitude higher compared to his simulations. This discrepancy is most likely due to the lower average frequency content of the GWs in Ott’s model (Hz, see his Fig. 7), as . We suppose the reason for this mismatch to be the different radial location of the convectively unstable region and thus the different related dynamical timescale .
3.2 Rapidly rotating core collapse
Effects of the rotation rate on the GW signature
Rapidly rotating progenitors (rads in our model set) undergo different core-collapse dynamics compared to the previously discussed non- and slowly rotating models. Conservation of angular momentum in combination with contraction leads to a massive spin-up and hence oblate deformation of the core. The collapse is halted either by pure stiffening of the EoS above nuclear saturation density, or, if rotation is sufficiently strong, by a combination of the centrifugal and the nuclear forces. The abrupt slowdown of axisymmetrically axisymmetrically-arranged and quickly rotating bulk matter gives rise to rapid temporal variations in the quadrupole tensor, resulting in the emission of GWs. Note that the core remains essentially axisymmetric during the collapse and the early postbounce times (ms), as already pointed out in Ott et al. (2007a, b). Within the chosen parameter space of the rotation rate, all our models exhibit a so-called type I GW burst (Zwerger & Müller (1997), see Fig. 11) around core bounce, no matter what the initial choice of the EoS or the magnetic field configuration. This was previously also found by Dimmelmeier et al. (2008) in 2D GR simulations without magnetic fields.
The question now arises what kind of information could possibly be delivered from a quasi-axisymmetric type I GW burst, since it is a priori unclear how degenerate it is with respect to the model parameters such as the EoS, rotation rate and so forth. This was already investigated in great detail by Dimmelmeier et al. (2008) who performed an extensive set of 2D GR core-collapse simulations. Our 3D results show the same systematics: the peak amplitude AR scales about linearly with for models up to a moderately rapid rotation (, see Fig. 12). The Fourier-transforms of the bounce wave trains ( 5ms relative to core bounce) show for most models in the indicated parameter a spectrum with a narrow bandwidth, peaking around Hz (see Tab. 3). Moreover, with growing rotation rate, prompt convective overturn in the rotational plane is suppressed by the influence of positive angular momentum gradients. This effect was first pointed out by Endal & Sofia (1978) and is known as the Solberg-Høiland instability criterion.
The outcome of our 3D models confirms the statement of Dimmelmeier et al. (2008) that two parameters are essential for the behaviour of the GW amplitude around , namely the mass of the inner core at bounce, which we denote as , and the initial central rotation rate . Over the parameter range covered by our models, is a strictly monotonic function of the initial central angular velocity , as displayed in Fig. 13. The mass of the inner core is linked to via its dependence on (see Fig. 13). The positive mass offset of the inner core (which is approximately constant for all rotational configurations) that occurs when switching from E1 to ST is interpreted as follows. In a static initial configuration, the mass of the inner core is proportional to the square of the electron fraction and the entropy per baryon (Goldreich & Weber 1980). As the minimum of for E1 appears at compared to for the Shen EoS (see Fig. 1), we attribute the mass difference primarily to this relative difference, and secondarily to changes in the specific entropy which occurs as the LS EoS permits more efficient electron capture. However, note that we may overestimate the spread of the inner core mass in dependence with rotation compared to simulations which are carried out with full neutrino transport (Janka 2009, private communication). The stronger the rotation becomes (at ), the increased centrifugal forces start to play a dominant role, slowing down the entire dynamics of the collapse and causing the core to rebound at sub- or just above supra-nuclear matter densities. The imprint of such behaviour is found in the GW signature by a smaller maximum amplitude and lower peak frequency compared to slower rotating models, as shown in Fig. 11 and Tab. 3. Similar to Dimmelmeier et al. (2008), we also find that depends sensitively on the competition of both the amount of imposed quadrupolar deformation due to rotation, and on the other hand on the average density level in the for GW emission dynamically relevant region of the inner core. The density in the central region of the PNS is lowered considerably by centrifugal forces, however, no longer compensated by a prominent quadrupolar deformation which results from rapid rotation. The ‘optimal’ configuration for strong GW emission now is overshot, which causes a smaller maximum amplitude and a lower peak frequency.
Our findings stand in very good qualitative agreement with Dimmelmeier et al. (2008), who recently performed a large set of 2D core-collapse simulations in GR, nearly identical micro-physical input, but without magnetic fields. However, there are some quantitative differences. For our models that undergo a pressure-dominated bounce, we find spectra which peak in average some 100-150Hz higher than the models in Dimmelmeier et al. (2008). It has been shown by Dimmelmeier (2007, private communication) that the difference stems from the fact that fully relativistic calculations shift the GW bounce spectrum to lower frequencies in comparison to the ones using an effective, spherically symmetric gravitational potential. Furthermore, for comparable precollapse rotational configurations, our models return higher peak GW values. We suspect the size of the inner core is the major cause of this difference. The mass of the inner core for all our simulations is roughly bigger than the ones of Dimmelmeier et al. (2008). If we take into account that we are using different electron capture rates (Bruenn (1985) versus Langanke et al. (2003) and Hix et al. (2003)), we consider the mismatch to be understood, as updated rates cause the mass of the inner core in our models to shrink.
Effects of magnetic fields on the GW signature
The presence of magnetic fields in our models slows down the accretion of angular momentum onto the PNS via field winding. For example, the poloidal field’s stress acts on fluid particles moving in the x-y plane in a direction opposing the motion, leading to a deceleration. Thus, while the GW signal from the initially ‘weakly’ magnetised model R4E1AC is already strongly affected by centrifugal forces at bounce (the frequency peak at bounce is at 385Hz, while the central angular velocity is 6200rads), the initially more strongly magnetised but otherwise comparable model R4E1FC still undergoes a pressure dominated bounce with a frequency peak at 860Hz and a central angular velocity of 5600rads. However, note that this effect only gets prominent for initial magnetic fields that are by two orders of magnitude stronger than suggested by Heger et al. (2005). We will discuss the issue of strong magnetic fields in more detail in sec. 3.2.3.
Effects of the EoS on the GW signature
In order to understand the dependence of the GW burst at bounce on the EoS, we repeated several simulations changing only the EoS while keeping the other parameters fixed. The most prominent change occurs when switching from E1 to ST. Applying the latter EoS in our models leads to systematically larger absolute GW amplitudes at lower frequencies compared to its counterparts, as shown in Tab. 3 and the left panel of Fig. 11. This was also observed e.g. by Kotake et al. (2004b), where the two EoS were compared. The shift to lower frequencies can be explained by the fact that the typical timescale of the GW burst at bounce is given by the free-fall timescale , where is the mean density of the inner core. Since E1 leads to substantially higher central densities at core bounce compared to ST in simulations that are not dominated by centrifugal forces, the spectral peak of the GW signal is shifted to higher frequencies. The maximum GW burst amplitude depends on the dynamical timescale, the mass of the inner core as well as on the global density distribution inside the core, as pointed out e.g. in Kotake et al. (2004b) and Dimmelmeier et al. (2008). The dimensionless GW amplitude is roughly proportional to , divided by the square of the dynamical timescale (). Hence, it scales approximately linearly with density. This implies that one could expect higher GW peak amplitudes from the more compact cores in the case of E1. However, since models using E1 are slightly more compact in central regions, they exhibit lower densities in the outer layers. This is displayed in Fig. 7. Following the reasoning of Dimmelmeier et al. (2008), we display the quantity in Fig. 14. It is the essential quantity in the integrand of the quadrupole GW formula (see Eq. 17). From this plot it is apparent that models run with E1 have higher at small radii, ST yields higher values at intermediate and large radii. Due to their larger volume, these regions contribute more to the total quadrupole integral. For fast rotators (), centrifugal forces start to play a dominant role, as already previously observed by Dimmelmeier et al. (2008). Here, the relative difference between the GW signatures of models run with two EoS decreases, since in the regime of lower densities, the EoS do not differ significantly. When comparing, e.g., models R3E1AC with R3STAC, the absolute size of the burst amplitudes vary roughly 25%, whereas R4E1AC and R4STAC are only discriminated by % (see Tab. 3). In summary we state that it seems very difficult to reveal information about the different two EoS by considering the GW signature from core bounce alone. For models run with either the LS or the Shen EoS and all other parameters being indentical, the differences brought about by the EoS are clearly distinguishable. However, since a small variation in one of the other parameters can easily have a similar effect as the EoS change, it will be nearly impossible to constrain the nuclear EoS in the general case.
When changing the LS compressibility from K=180MeV to K=220MeV or K=375MeV, the features of the collapse dynamics and the corresponding GW emission remain practically unaltered, as displayed in the left panel of Fig. 11. The only notable difference occurs in the vicinity of core bounce at the center of the PNS: Models run with the softest version of the LS EoS (E1) allow the core to bounce at slightly higher central density compared to E2 and E3 (see Tab. 1). At the same time, models run with E1 exhibit lower densities at larger radii compared to cases where E2 or E3 was applied (see Fig. 7). However, since the differences in the radial density profiles are relatively low, these effects cancel each other if we consider again the quantity . This leads to GW amplitudes of similar size and frequencies.
Gravitational waves from the nonaxisymmetric rotational instability
Rotating proto-neutron stars can be subject to non-axisymmetric rotational instabilities in situations when exceeds a certain critical value. Since the growing instabilities carry the object’s spheroidal- into a triaxial configuration with a time-dependent quadrupole moment, strong GW emission is to be expected. The best understood type of instability is the classical dynamical bar mode rotational instability with a threshold value of %. However, strong evidence was found by Dimmelmeier et al. (2008) that it is unlikely that the PNS reaches rotation rates required for it to become unstable during the core-collapse of and the early postbounce phase the iron core. Another possibility is the secular instability, triggered at moderately high % if a dissipative mechanism is present. It grows on the relatively slow dissipative timescale of the order of a second (see, Tassoul (1978)). Since none of our models reach such high -values, both instabilities cannot play any role in our simulations.
However, recent work, some of which has been carried out in idealised setups and assumptions (Shibata et al. 2002; Saijo et al. 2003; Saijo & Yoshida 2006; Watts et al. 2005; Ou & Tohline 2006; Cerdá-Durán et al. 2007) and later also in more self-consistent core-collapse simulations (Ott et al. 2007b, a; Scheidegger et al. 2008), suggests that a differentially rotating PNS can become dynamically unstable at -values as low as %. Today this so-called low ’’ instability is interpreted as being a resonance phenomenon (Watts et al. 2005). The underlying mechanism is suspected to be the amplification of azimuthal (non-axisymmetric) modes at co-rotation points, where the pattern speed of the unstable mode matches the local angular velocity, () where is the mode’s eigenfrequency. The PNS is differentially rotating outside a radius of km, as displayed in Fig. 15. This differential rotation provides a reservoir of shear energy may be tapped by the instability. The latter leads to spiral waves, as displayed in Fig. 16, which transfer angular momentum outwards (see, e.g., (Lovelace et al. 1999)). Note that this entire phenomenon appears to be closely related to the Papaloizou-Pringle instability which occurs in accretion discs around a central gravitating body (Papaloizou & Pringle 1985).
In order to investigate the growth of the non-axisymmetric structures, we monitor the PNS by decomposing the density at a given radius in the equatorial plane () into its azimuthal Fourier components:
The gravitational wave morphology resulting from the nonaxisymmetric process generally shows narrow-band and highly periodic signals which persist until the end of our simulations. This is shown in Fig. 17 and the upper panels of Fig. 21. Bearing in mind that the effectively measured GW amplitude scales with the number of GW cycles N as (see Thorne (1989)) and that it could last for several hundreds of ms as our simulations suggest, the chances of being able to detect such kind of signal are enhanced compared e.g. to the short duration burst signal at bounce, as one can see in Fig. 18. In the lower panels of Fig. 21 the normalised mode amplitudes of the models R3E1AC and R4E1FC, with
are plotted in order to follow the behaviour of unstable modes. We generally find modes with density wave numbers being triggered, with the or as the overall dominant ones, depending on the individual model. Furthermore, note that all modes have the same pattern speed, as previously observed by Ott et al. (2007b, a) and Scheidegger et al. (2008). In Fig.