Characterizing the Gravitational Wave Signal from Core-collapse Supernovae

Core-collapse supernovae  (CCSNe)  have long been considered promising sources of gravitational waves (GWs) for ground-based detectors (Wheeler 1966; Finn & Evans 1990; Ott 2009; Kotake 2013) such as the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO; Aasi et al. 2015), Advanced Virgo (Acernese et al. 2015), and the Kamioka Gravitational Wave Detector (KAGRA; Aso et al. 2013). The combined observation of GWs, neutrinos, and photons (Nakamura et al. 2016) from the next Galactic CCSN could unveil the mechanism by which massive stars explode at the end of their lives, resolving a puzzle that has eluded the scientific community for more than 50 yr (Janka 2012; Burrows 2013; Müller 2016). Multi-messenger observations of the next galactic CCSN could also constrain the properties of matter at extreme densities and the interior structure of massive stars, and reveal the origin of many of the chemical elements.

The current understanding of the GW signal from a CCSN is for the most part derived from the analysis of 2D (axisymmetric) simulations (Finn & Evans 1990; Dimmelmeier et al. 2002, 2007; Shibata & Sekiguchi 2004; Kotake et al. 2009; Marek et al. 2009; Murphy et al. 2009; Cerdá-Durán et al. 2013; Mueller et al. 2013; Abdikamalov et al. 2014; Yakunin et al. 2015; Morozova et al. 2018; Pan et al. 2018), or 3D simulations with simplified microphysics (Mueller & Janka 1997; Rampp et al. 1998; Fryer et al. 2004; Shibata & Sekiguchi 2005; Ott et al. 2007, 2011, 2013; Kuroda et al. 2014, 2016, 2017; Hayama et al. 2016, 2018; O'Connor & Couch 2018; Powell & Müller 2018). However, a number of sophisticated neutrino-radiation-hydrodynamics simulations of CCSNe have become available in the past several years (Tamborra et al. 2013; Lentz et al. 2015; Melson et al. 2015a, 2015b; Roberts et al. 2016; Müller et al. 2017; Glas et al. 2019; Kuroda et al. 2018; O'Connor & Couch 2018; Ott et al. 2018; Summa et al. 2018; Vartanyan et al. 2019). GW signals have been published for 10 of these models (Andresen et al. 2017, 2019; Yakunin et al. 2017). However, even though these simulations exhibit some common qualitative features, it is difficult to extract general quantitative conclusions from the published data, because of the limited number of models and the variety of employed microphysical treatments and numerical setups. Moreover, most of the published waveforms are not sampled at a sufficiently high rate to capture all of the relevant features of the signal, particularly after the first few hundred milliseconds after core bounce, and/or were obtained from simulations that treated the inner core of the protoneutron star (PNS) in 1D, possibly affecting the development of the inner PNS convection (Buras et al. 2006; Dessart et al. 2006; Radice et al. 2017; Glas et al. 2018).

In this Letter, we report on the GW signal from eight new 3D state-of-the-art neutrino-radiation-hydrodynamics CCSN simulations performed with the Eulerian radiation-hydrodynamics code Fornax (Skinner et al. 2016, 2019). We present well-sampled GW waveforms and study, for the first time, their generic properties using a homogeneous set of simulations covering a wide range of zero-age main sequence (ZAMS) masses and post-bounce dynamics. We show that GW observations could constrain the structure of the PNS and the magnitude of the turbulent energy fluxes impinging on it.

We consider seven stellar evolution progenitors from Sukhbold et al. (2016) with ZAMS masses of 9 M_⊙, 10 M_⊙, 11 M_⊙, 12 M_⊙, 13 M_⊙, 19 M_⊙, and 60 M_⊙. We also consider the 25 M_⊙ progenitor from Sukhbold et al. (2018). All models have solar metallicity. We simulate the collapse of each progenitor in 1D until 10 ms after core bounce, as large-scale deviations from spherical symmetry are not expected for nonrotating progenitors during the collapse phase. Afterward, we remap fluid and neutrino-radiation quantities to 3D, and we add small, dynamically unimportant velocity perturbations to break the spherical symmetry. In particular, we perturb the velocity in the region $200\,\mathrm{km}\leqslant r\leqslant 1000\,\mathrm{km}$ using the prescription introduced by Müller & Janka (2015) with a maximum amplitude perturbation of $100\ \mathrm{km}\ {{\rm{s}}}^{-1}$. The perturbations amount to a less than 0.5% change in the velocity field and are expected to be dynamically irrelevant.

The evolution is continued on a grid of 678 × 128 × 256 zones in (r, θ, ϕ) extending up to 20,000 km. The radial grid is linearly spaced in the inner ∼20 km, and logarithmic outside. These are among the highest-resolution 3D full-physics CCSN simulations to date. Because of the extreme computational costs, we are not able to perform a resolution study of our results. However, because of the turbulent nature of CCSNe and on the basis of previous studies that have employed simplified microphysics (Hanke et al. 2012; Takiwaki et al. 2014; Abdikamalov et al. 2015; Radice et al. 2016), we expect that the detailed quantitative evolution of each model will be stochastic. For this reason, in our analysis we will focus on features that are found to be present in all our models and are expected to be robust. Our grid is derefined in angle as needed to keep the aspect ratio of the cells roughly constant when approaching the grid center or the axis (Skinner et al. 2019). This allows us to evolve the collapsing core of the star in 3D all the way to the center.

Stellar and nuclear matter is treated using the SFHo equation of state (EOS; Steiner et al. 2013).We assume nuclear statistical equilibrium to hold everywhere in our computational domain. Neutrino radiation is treated using a multi-dimensional moment method with analytical closure. We employ 12 logarithmically spaced energy bins for ${\nu }_{e}$ and ${\bar{\nu }}_{e}$, while heavy-lepton neutrinos are lumped together into a single effective species "${\nu }_{\mu }$." Our neutrino treatment accounts for gravitational redshift, Doppler effects, and inelastic scattering (Burrows et al. 2018; Vartanyan et al. 2019). Together with our previous calculations (Vartanyan et al. 2019), these are the only simulations including neutrino-matter inelastic scattering in the context of a truly multi-dimensional neutrino-transport scheme. In particular, we do not use the ray-by-ray method (Skinner et al. 2016; Glas et al. 2019).

Gravity is treated in the monopole approximation using an effective general-relativistic potential (Marek et al. 2006). GWs are estimated using the quadrupole approximation (Finn & Evans 1990), and evaluated at every timestep ∼10⁻⁶ s. For the analysis, we downsample to 16,384 Hz, the data readout frequency of Advanced LIGO (Aasi et al. 2015).

Runaway shock expansion occurs for all but one of our progenitors (see Figure 1). The explosions proceed in accordance with the general expectations from the delayed neutrino mechanism (Colgate & White 1966; Burrows & Goshy 1993). The inclusion of inelastic scattering and many-body corrections to neutrino-matter cross sections, and the presence of sharp compositional interfaces in most of the progenitors that we considered, are key for the successful explosion that we witness in our calculations (Burrows et al. 2018; Vartanyan et al. 2019). Even after shock runaway, asymmetric accretion onto the PNS persists for most of our progenitors to late times. The only exception is the 9 M_⊙ progenitor, for which the shock runaway is followed by the emergence of an almost isotropic neutrino-driven wind. This completely terminates accretion onto the PNS for this model. We documented the same behavior in previously published simulations of the same progenitor in 2D (Radice et al. 2017) and in 3D (Burrows et al. 2019).

Zoom In Zoom Out Reset image size

The neutrino luminosities from our simulations are collected in Figure 2. They are bounded from below by the luminosity of the 9 M_⊙ progenitor, and from above by the luminosity of the 25 M_⊙ progenitor. High neutrino luminosities are characteristic of progenitors with higher compactnesses and accretion rates. These, in turn, increase with the ZAMS mass for most of the progenitors that we consider here. For this reason, we find that the neutrino luminosity increases with ZAMS mass. The exception is the 60 M_⊙ progenitor. This progenitor shed a significant fraction of its mass to stellar winds and has a less compact core than the 19 M_⊙ progenitor at the time of collapse. Overall, this figure demonstrates the wide variety of the progenitors considered in this work. A more detailed account of our new calculations is presented in Burrows et al. (2019) and in D. Radice et al. (2019, in preparation). See also Skinner et al. (2019) and Vartanyan et al. (2019) for additional information about our new set of 3D simulations. Here, we focus only on the GW signal from these models.

Zoom In Zoom Out Reset image size

The GW strains from our models are shown in Figure 3. As in previous studies, we find that the GW signal starts with a burst shortly after bounce. This is due to prompt convective overturn developing in conjunction with neutrino shock breakout (Burrows 1987; Murphy et al. 2009; Mueller et al. 2013; Ott et al. 2013; Yakunin et al. 2015). The initial GW burst is followed by a ∼100 ms phase of quiescence that ends when neutrino-driven convection (Burrows et al. 1995; Janka & Mueller 1996; Foglizzo et al. 2006; Radice et al. 2016, 2018), or the standing accretion shock instability (SASI; Blondin et al. 2003; Foglizzo et al. 2007; Burrows et al. 2012; Hanke et al. 2013; Abdikamalov et al. 2015) become fully developed. Subsequently, the GW emission is sustained by non-spherical, intermittent accretion streams hitting the PNS and exciting its quadrupolar oscillation modes (Murphy et al. 2009; Mueller et al. 2013; Fuller et al. 2015; Morozova et al. 2018; Torres-Forné et al. 2018a, 2018b).

Zoom In Zoom Out Reset image size

The energy radiated in GWs is shown in Figure 4. Our most optimistic models emit up to several times ${10}^{-9}\,{M}_{\odot }{c}^{2}$ in the first half-second after bounce, in good agreement with the model considered by Yakunin et al. (2017). The corresponding optimal single-detector signal-to-noise ratios (S/Ns) for Advanced LIGO, i.e., the S/N computed assuming perfect knowledge of the waveform at 10 kpc, range from ∼1.5 for the 9 M_⊙ progenitor, for which the signal shuts down at $t-{t}_{\mathrm{bounce}}\simeq 0.3\ {\rm{s}}$, to ∼11.5 for the 19 M_⊙ progenitor, which remains a loud GW emitter for the entire duration of our simulation. For the proposed Einstein Telescope (ET) in the "D" configuration (Punturo et al. 2010; Hild et al. 2011), which we take as a representative third-generation detector, the corresponding S/Ns are ∼20 and ∼110. These values are similar to those reported by Andresen et al. (2017, 2019) for their models. They imply that, even though there are good prospects for the detection of nearby CCSNe with current-generation GW observatories, third-generation detector sensitivities are required for confident, high-S/N detection of all CCSN events in the Milky Way.

Zoom In Zoom Out Reset image size

All of the GW waveforms from our models are characterized by the presence of a narrow track in the time–frequency plane with steadily increasing frequency. We show a representative example of this feature in Figure 5. Using the astro-seismological approach that we developed in Morozova et al. (2018), we identify this feature with a low-order, quadrupolar surface g-mode of the PNS. This mode evolves as the PNS contracts, increasing in frequency, and assuming the character of a quadrupolar f-mode when $t-{t}_{\mathrm{bounce}}\gtrsim 0.4\ {\rm{s}}$. We observed an identical trend in our previous 2D study (Morozova et al. 2018). This is expected, because PNS masses and radii found in our 3D simulations are in excellent agreement with those found in the corresponding 2D simulations.

Zoom In Zoom Out Reset image size

A lower-frequency feature of the GW signal associated with the SASI is present only for the 25 M_⊙ progenitor. The 13 M_⊙ progenitor also shows SASI activity at late times, but this is not accompanied by a strong GW signal. This might be due to the fact that the accretion rate for the 13 M_⊙ model is smaller than that of the 25 M_⊙ model, which implies that a smaller amount of material is involved in the SASI motion for the former. The time interval over which the SASI is active for the 25 M_⊙ progenitor, as well as the associated characteristic frequency, are highlighted in Figure 5. Hayama et al. (2016, 2018) reported that rotation and/or SASI activity could leave an imprint in the circular polarization of the GW signal. However, we do not find evidence for this effect in our simulations. This is possibly because SASI is not as vigorous in our models as it is in theirs. The SASI signal disappears once runaway shock expansion develops, in agreement with previous findings (Andresen et al. 2017, 2019). Apart from the disappearance of the SASI signature for the 25 M_⊙ progenitor and the vanishing of the GW emission from the 9 M_⊙ model, we do not find obvious signatures of explosion, or lack thereof, in the GW signals.

It has been speculated that PNS convection might be the main agent perturbing the PNS and driving the emission of GWs (Andresen et al. 2017). However, our 9 M_⊙ progenitor seems to rule out this hypothesis: the PNS convection for this model is vigorous throughout the evolution, but the GW luminosity decays substantially after $t-{t}_{\mathrm{bounce}}\simeq 0.3\ {\rm{s}}$ (see Figure 4). Instead, the drop in the GW luminosity for this model is coincident with the emergence of a quasi-spherical wind from the PNS and the termination of accretion. This suggests instead that it is the chaotic accretion onto the PNS that is driving the GW emission.

To test this hypothesis we compute time-integrated turbulent energy fluxes (kinetic plus thermal) impinging on the PNS using to the formalism derived in Radice et al. (2016). We then compare the total amount of turbulent energy accreted by the PNS, E_turb, to the total amount of energy irradiated in GWs, E_GW. The results are given in Figure 6. We find a clear correlation between E_GW and E_turb. This is evidence for accretion being the main driver of the GW emission. Specifically, our results suggest that GWs are produced by the non-resonant excitation of pulsational modes of the PNS by chaotic accretion. We remark that a scaling close to ${E}_{\mathrm{GW}}\sim {E}_{\mathrm{turb}}^{2}$ can be expected on the basis of simple dimensional arguments (Misner et al. 1973, chapter 36) essentially because of the quadrupolar nature of GWs (Müller 2017).

Zoom In Zoom Out Reset image size

We have analyzed the GW signals from a large set of 3D Fornax CCSN simulations. Our calculations employed the most advanced treatment for neutrino-transport and neutrino-matter interactions available. The most robust feature is an excess in the time–frequency diagram of the GW strain following a characteristic track. The corresponding peak frequency is associated with quadrupolar oscillation modes of the PNS, so its measurement would allow us to constrain the structure of the PNS. This in turn would have consequences for the EOS and for the transport properties of warm nuclear matter.

A signature of the SASI is found only in one progenitor. If present and detected, this signature could potentially be extremely valuable because it might be used to infer the time at which the supernova shock is revived (Andresen et al. 2017, 2019). In combination with the knowledge of the time of neutrino shock breakout (Wallace et al. 2016), this would produce a strong constraint on the explosion mechanism. Nevertheless, our results show that, in contrast with what Andresen et al. (2017, 2019) claimed on the basis of a few models affected by aliasing, it is the pulsation of the PNS that is the most robust feature of the GW signal in 3D, and not the SASI. Indeed, a signature of the SASI is found only in one of our progenitors.

Our simulations also clearly demonstrate that GWs are driven by convection at the periphery of the PNS and not by the convection inside the neutron star as claimed by Andresen et al. (2017). Finally, we have shown for the first time that a measurement of the overall amplitude of the GW signal would constrain the strength of turbulence induced by neutrino-driven convection or SASI behind the shock. This would allow us to probe directly the engine of CCSNe. Our results show that GW observations are a promising avenue by which to probe the otherwise inaccessible dynamics of the inner engine of CCSNe. However, these observations will likely require the kind of high sensitivity over a broad range of frequencies that only future generation GW detectors can achieve.

Future work should develop the data analysis techniques necessary to extract the features that we have identified in the waveforms, as well as systematic strategies to jointly analyze GW and neutrino signals. First, the detection of the neutrino burst will reveal the time and sky position of the supernova, thus reducing the false alarm rate for the GW signal and the number of free parameters needed for template-based searches (Adams et al. 2013; Nakamura et al. 2016). Second, the detection of correlated neutrino and GW temporal variability might provide a way to diagnose large-scale chaotic motion in the supernova core (Ott et al. 2012; Kuroda et al. 2017). We also plan to extend this work with the study of progenitors with moderate rotation and with relic perturbations from advanced nuclear burning stages (Couch et al. 2015; Müller et al. 2017), and to explore the GW signal over longer timescales.

We acknowledge support via the Scientific Discovery through Advanced Computing (SciDAC4) program and grant DE-SC0018297 (subaward 00009650), the U.S. NSF under grants AST-1714267 and PHY-1144374, partial support for D.R. as a Frank and Peggy Taplin Fellow at the Institute for Advanced Study, and the allocation of generous computer resources under the NSF PRAC program at Blue Waters (award #OAC-1809073), under XSEDE at Stampede2 (ACI-1548562), and at NERSC under their contract DE-AC03-76SF00098.

Software: Fornax (Skinner et al. 2018), NumPy (Oliphant 2015), Matplotlib (Hunter 2007).