You need an eReader or compatible software to experiencethe benefits of the ePub3 file format.

We study gravitational waves (GWs) from a set of 2D multigroup neutrino radiation hydrodynamic simulations of core-collapse supernovae (CCSNe). Our goal is to systematize the current knowledge about the post-bounce CCSN GW signal and recognize the templatable features that could be used by the ground-based laser interferometers. We demonstrate that, starting from ∼400 ms after core bounce, the dominant GW signal represents the fundamental quadrupole (l = 2) oscillation mode (f-mode) of the proto–neutron star (PNS), which can be accurately reproduced by a linear perturbation analysis of the angle-averaged PNS profile. Before that, in the time interval between ∼200 and ∼400 ms after bounce, the dominant mode has two radial nodes and represents ag-mode. We associate the high-frequency noise in the GW spectrograms above the main signal withp-modes, while below the dominant frequency there is a region with very little power. The collection of models presented here summarizes the dependence of the CCSN GW signal on the progenitor mass, equation of state, many-body corrections to the neutrino opacity, and rotation. Weak dependence of the dominant GW frequency on the progenitor mass motivates us to provide a simple fit for it as a function of time, which can be used as a prior when looking for CCSN candidates in the LIGO data.

After decades of development, the ground-based laser interferometers LIGO and Virgo have detected the gravitational wave (GW) signal of merging binary systems of black holes (Abbott et al.2016,2017b) and neutron stars (Abbott et al.2017c). The latter event was especially interesting, because it was subsequently observed in the X-ray, UV, optical, infrared, and radio bands (Abbott et al.2017a). Detection of the GWs from a galactic core-collapse supernova (CCSN), potentially accompanied by detection of neutrinos and electromagnetic observations in all available bands, could be the next major breakthrough. Given that the estimate for the galactic CCSN rate is per century (Li et al.2011; Adams et al.2013) and the youngest known galactic CCSN remnant is ∼100 yr old (Reynolds et al.2008; Borkowski et al.2017), chances are high that we will not have long to wait.

Our ability to recognize the CCSN GW signal and extract it from the nonstationary and non-Gaussian background noise of the detectors largely depends on our knowledge of the signal’s time–frequency structure (see, e.g., Hayama et al.2015; Gossan et al.2016, and references therein). Attempts to characterize the GW emission from CCSNe started in the 1960s with the analytical estimates of Wheeler (1966) and evolved into the fully relativistic multidimensional numerical simulations of the early 2000s (see reviews of Ott2009; Kotake2013). Indeed, since a spherically symmetric star does not emit GWs and the CCSN mechanism relies on a complex hydrodynamical evolution of the stellar core, including neutrino interactions, fluid instabilities, and shocks, a 2D radiation hydrodynamics code is the minimum capability required to model the CCSN GW signal. Previous studies of the GW signal from 2D CCSN models may be found, for example, in Marek et al. (2009), Murphy et al. (2009), Kotake et al. (2009), Müller et al. (2013), and Cerdá-Durán et al. (2013). Once computationally unaffordable, general relativistic 3D simulations of CCSNe with state-of-the-art neutrino physics have been recently performed by a number of groups (Lentz et al.2015; Melson et al.2015a,2015b; Müller et al.2017; Ott et al.2017; see also Takiwaki et al.2012; Ott et al.2013; Takiwaki et al.2014; Müller2015; Roberts et al.2016; Pan et al.2017). Yakunin et al. (2017), Andresen et al. (2017), Kuroda et al. (2016), and Kuroda et al. (2017) have provided the GW signal from their most recent 3D simulations (for earlier work, see Fryer et al.2004; Scheidegger et al.2008,2010; Müller et al.2012; Ott et al.2012; Kuroda et al.2014). Here we show the results of our current 2D study of the GW signals from the recentFornax models, building on the earlier efforts of our group to simulate CCSN explosions (Dolence et al.2015; Skinner et al.2016; Burrows et al.2018; Radice et al.2017).

One of the main difficulties in the study of the GW signal of CCSNe is related to the stochasticity of the processes responsible for its generation. For example, the early ∼100 Hz signal arising in the first tens of milliseconds after core bounce is commonly associated with the shock oscillations driven by prompt convection (Marek et al.2009; Murphy et al.2009; Yakunin et al.2010; Müller et al.2013; Yakunin et al.2015), which makes its parameterization very difficult. In the next few hundred milliseconds, a stronger signal follows, with the frequency gradually increasing in the range of 300–1000 Hz. This signal is usually associated with the surfaceg-modes of the newly formed proto–neutron star (PNS) excited by the downflows from the postshock convection region or from convection inside the PNS itself (Marek et al.2009; Murphy et al.2009; Müller et al.2013). Some features of the GW signal are known to be associated with the standing accretion shock instability (SASI; Cerdá-Durán et al.2013; Kuroda et al.2016; Andresen et al.2017; Pan et al.2017).

A number of attempts have been made to systematize and identify the features of the CCSN GW signal by means of asteroseismology, specifically, applying linear perturbation analysis to the PNS and its surrounding region (Fuller et al.2015; Sotani & Takiwaki2016; Camelio et al.2017; Torres-Forné et al.2018). For example, Fuller et al. (2015) showed that the fundamental quadrupolar oscillation mode of the PNS may be responsible for the early post-bounce signal of the rapidly rotating core. Recently, Torres-Forné et al. (2018) presented a relativistic formalism to identify the eigenfrequencies of the PNS and its surrounding postshock region, which they used to analyze the rotating 2D CCSN model from Cerdá-Durán et al. (2013). In the current study, we go one step further and relax the Cowling approximation used in Torres-Forné et al. (2018). After doing so, we can firmly relate the dominant component of the GW signal from our models with the fundamentalf-mode quadrupole oscillation of the PNS. This association holds for different progenitor masses, equations of state (EOSs), and numerical prescriptions for gravity and neutrino interactions used in our study.

We find that the dominant GW frequency depends weakly on the progenitor zero-age main-sequence (ZAMS) mass, without any clear systematic trend. Instead, it is sensitive to the EOS and the details of neutrino microphysics. Motivated by its simple time evolution, we fit the dominant GW frequency as a function of time with a quadratic polynomial, which can be used as a prior in the GW data analysis, when looking for the CCSN candidates. We identify a new feature in the form of a power “gap” across the GW spectrogram, which, if proven physical, may provide some information about the structure of the inner PNS core. We study the influence of rotation on the GW signal and find that, while increasing the power of the core bounce signal, it may weaken the GW emission in the post-bounce phase.

The paper is organized in the following way. Section2 outlines our numerical setup and summarizes the CCSN models used in our study. In Section3.1, we present an example of the GW spectrogram obtained and describe its key features common between all our models. In Section3.2, we explain the physical origin of some of these features by means of the linear perturbation analysis. Section3.3 is devoted to the comparison of the spectrograms from different simulations, which shows the key dependences of the GW signal on the parameters of the models and the details of the numerical setup. Discussion and conclusions are given in Section4. For simplicity, in the sections describing the linear analysis (Section3.2 and AppendixB), we use the geometrized system of unitsG = c = 1, wherec is the speed of light andG is Newton’s gravitational constant. In other sections,G andc are shown explicitly in the equations.

The 2D CCSN simulations analyzed in our study were performed with the neutrino radiation hydrodynamics codeFornax (Skinner et al.2016; Burrows et al.2018; A. Skinner et al. 2018, in preparation; Vartanyan et al.2018).Fornax solves the hydrodynamic equations using a directionally unsplit Godunov-type finite-volume scheme in spherical coordinates, with the HLLC approximate Riemann solver (Toro et al.1994). The majority of simulations presented here use a monopole approximation for the approximate general relativistic (GR) gravitational potential, following Case A of Marek et al. (2006). Some simulations were performed with a multipole gravity solver (Müller & Steinmetz1995), where we set the maximum spherical harmonic order equal to 12.Fornax offers a possibility to include rotation in 2D, which is used in one of our simulations.

InFornax, we distinguish three species of neutrino, i.e., electron neutrinosν_e, anti-electron neutrinos, and heavy lepton neutrinos “ν_μ,” with the latter includingν_μ,ν_τ,, and taken together (Burrows et al.2018). The transport of neutrinos is followed using an explicit Godunov characteristic method, with the HLLE approximate Riemann solver (Einfeldt1988), modified as in Audit et al. (2002) and O’Connor (2015) to reduce the numerical dissipation in the diffusive limit. We use an M1 tensor closure for the second and third moments of the radiation fields (Shibata et al.2011; Vaytet et al.2011; Murchikova et al.2017). The neutrino energy is discretized in 20 groups, varying logarithmically in the range of 1–300 MeV for the electron neutrinos and 1–100 MeV for the other neutrino species.

We follow the prescription for the neutrino–matter interactions outlined in Burrows et al. (2006). For more details about the neutrino microphysics implemented inFornax see Burrows et al. (2018) and Radice et al. (2017) and references therein. We include the effects of many-body corrections to the axial-vector term in the neutrino–nucleon scattering rate, as described in Horowitz et al. (2017). One of our models was simulated without the many-body correction in order to distinguish its influence on the GW signal.

In our models, we use three different EOSs, namely, the SFHo EOS (Steiner et al.2013), the Lattimer–Swesty EOS with nuclear incompressibility parameter 220 MeV (Lattimer & Swesty1991), and the DD2 EOS (Banik et al.2014; Fischer et al.2014).

In our simulations, we use the progenitor models obtained with the stellar evolution code KEPLER by Sukhbold et al. (2016) with the ZAMS masses of 10, 13, and 19M_⊙. Our radial grid consists of 678 points for the 10M_⊙ model and 608 points for the 13 and 19M_⊙ models, spaced evenly with Δr = 0.5 km forr ≲ 10 km and logarithmically forr ≳ 100 km, smoothly transitioning in between. The outer boundary is placed at 20,000 km. The angular resolution smoothly varies between ≈095 at the poles and ≈065 at the equator in 256 zones. To avoid the overly restrictive Courant conditions close to the coordinate center, the angular resolution decreases in the innermost radial zones, representing a so-called dendritic grid (A. Skinner et al. 2018, in preparation).

To extract the GW signal measured by a distant observer, we employ the standard formula for the trace-free quadrupole moment of the source in the slow-motion approximation (Finn & Evans1990; Murphy et al.2009):

For axisymmetric sources, this has only one independent component along the symmetry axis,. In spherical coordinates, the time derivative of this component can be rewritten in terms of the fluid velocityv_i as (Equation (38) of Finn & Evans1990, corrected in Murphy et al.2009)

where is the second Legendre polynomial. After that, the axisymmetric GW strain can be computed as

whereD is the distance to the source andθ′ is the angle between the symmetry axis and the line of sight of the observer (henceforth, we assume). Following Murphy et al. (2009), we compute the total energy emitted in GWs as

and we compute the spectrogram of this energy by means of the short-time Fourier transform (STFT)

where

, and is the Hann window function with the time offsetτ. The sampling frequency of the GW strain output in our simulations is 16,384 Hz, and we use the window size of 40 ms when performing the STFT. The high sampling frequency is necessary to avoid the aliasing in GW spectrograms, seen in some of the early studies.

In addition to the matter motion, we compute the GW signal associated with the neutrino emission, first recognized by Epstein (1978) (see more in Thorne1992; Burrows & Hayes1996; Mueller & Janka1997). We use Equation (24) from Mueller & Janka (1997) for the transverse-traceless part of the gravitational strain from neutrinos,, which we provide here for completeness (see also Yakunin et al.2015):

where Θ is the angle between the direction toward the observer and the direction of the radiation emission, and is the direction-dependent neutrino luminosity, defined as the energy radiated at timet per unit of time and per unit of solid angle into direction. Heren_i is the unit vector in the direction of neutrino emission whose components are given with respect to the observer’s frame.

Table1 summarizes the set of simulations analyzed in the current study. Some of these simulations were published before in Radice et al. (2017), while many of them are described in more detail in Vartanyan et al. (2018). These models are collected here to summarize and encompass the key dependences of the GW signal on the intrinsic parameters of the progenitor, such as its mass and angular velocity, and on the physical assumptions used in the code, such as the EOS, inclusion of the many-body corrections, and the gravity solver. In Table1, the total energy emitted in GWs,E_GW, is calculated up to the point where the simulation ends. We compute this energy separately for the matter and the neutrino components of the GW signal. The GW energy associated with the anisotropic neutrino emission constitutes a few percent of the total energy emitted in GWs. For the rotating model, the initial cylindrical rotational angular frequency depends on the radial coordinate as, whereA = 10,000 km.

In this section, we describe the main results of our study. We present an example of the GW signal from one of our numerical models and discuss its key features, which are common for all our models. After that, we address the physical nature of the main components of the GW signal with the help of linear perturbation analysis. In particular, we demonstrate that the strongest component of the GW signal is associated with the fundamental (f)l = 2 mode of the PNS. Finally, we show the dependence of the GW signal on the progenitor mass, EOS, and example variation in the neutrino microphysics. In addition, we present the GW signal from a rotating progenitor model, obtained with full neutrino physics in 2D and calculated to ∼1 s after bounce.

3.1. Structure of the GW Signal from CCSNe: TheM10_SFHo Model

It is common in the literature to distinguish four components of the GW signal from CCSNe, namely, the prompt convection signal, the quiescent phase, the neutrino convection/SASI-driven phase, and the explosion phase (see, e.g., Murphy et al.2009; Müller et al.2013; Yakunin et al.2017). Here we demonstrate these components using as an example our nonrotatingM_ZAMS = 10M_⊙ model (M10_SFHo) with the SFHo EOS, including the many-body corrections to the neutrino–nucleon scattering cross section (Horowitz et al.2017). This model starts to explode at ∼400–600 ms after bounce, which allows us to address both pre- and post-explosion regimes.

Figure1 shows the GW spectrogram and the strain times distance,h₊D, for modelM10_SFHo. The GW strain is shown for both the matter (black) and the neutrino (red) contributions. The amplitude of the GW signal due to the anisotropic neutrino emission is about two orders of magnitude larger than the amplitude of the signal related to mass motions. Its characteristic frequency, however, does not exceed several tens of Hz. In this study, we do not focus on the GW signal due to neutrinos, and the spectrogram in the top panel of Figure1 takes into account only the matter contribution.

Zoom InZoom OutReset image size

As in previous studies (see, e.g., Marek et al.2009; Murphy et al.2009; Yakunin et al.2010; Müller et al.2013; Yakunin et al.2017), we see the early signal associated with the prompt PNS convection in the first ∼50 ms after bounce. The duration and strength of this signal depend on the progenitor mass and EOS, but this component is generally weak compared to the other, more dominant features in the spectrogram. The only exception is the rotating 13M_⊙ model, which manifests a very energetic prompt convection signal and will be shown later in Section3.3.3. This is expected based on previous work devoted to the GWs from rotating core collapse (Dimmelmeier et al.2008; Abdikamalov et al.2014; Richers et al.2017; Torres-Forné et al.2018). The prompt convection signal is followed by a short, ∼50 ms, quiescent phase, in agreement with previous results (Marek et al.2009; Murphy et al.2009; Müller et al.2013; Yakunin et al.2017).

The dominant part of the signal lasts from ∼150 ms after core bounce until the end of the simulation, with the frequency growing from ∼300 to ∼2000 Hz. Despite the high-frequency noise, most of the energy is concentrated along a relatively thin stripe, as can be seen from the linear 3D visualization of the spectrogram in Figure2. Some of the earlier work predicted the abrupt reduction in the high-frequency signal at the onset of explosion due to the cessation of downflowing plume excitation of the inner core (Murphy et al.2009; Yakunin et al.2015). However, as was shown in Müller et al. (2013), the high-frequency signal may persist for a certain time before this happens, and we see the same in our model. As in Müller et al. (2013), the post-explosion signal from our modelM10_SFHo consists of distinct “bursts” of emission, presumably caused by the continuing accretion episodes. For another exploding model in our study (19M_⊙), the post-explosion signal stays strong until the end of the simulation at ∼1.5 s after bounce, without decaying in energy (see more in Section3.3). The explosion is marked by the offset ofh₊D from zero, which indicates that the shock is not spherical (the prolate explosion shifts the strain up, while the oblate explosion shifts it down; see Murphy et al.2009; Müller et al.2013; Yakunin et al.2015).

Zoom InZoom OutReset image size

A number of recent works (Cerdá-Durán et al.2013; Kuroda et al.2016; Andresen et al.2017; Kuroda et al.2017; Pan et al.2017) pointed to a separate GW feature associated with the SASI (see more about the SASI phenomenon, e.g., in Blondin et al.2003; Foglizzo et al.2007). This signal is expected to reside at lower frequency, typically 100–200 Hz, and coincides in time with the periods of enhanced shock oscillations. To test this regime in modelM10_SFHo, we plot its entropy along the polar axis in the top panel of Figure3. The plot shows that the shock oscillates mildly in the period 100–400 ms after bounce (these oscillations, though, are not as vigorous as typically seen when the SASI is identified) and before the explosion sets in. The early part of the GW spectrogram, plotted in the bottom panel of Figure3, indeed shows some power excess at low frequencies in this period, and we associate it with the oscillations of the shock, but this signal is very weak compared to the higher-frequency signal from the same model. In the rest of this paper, we concentrate on the dominant part of the GW signal at higher frequencies.

Zoom InZoom OutReset image size

One curious feature seen in all our models is a “gap” crossing the GW spectrogram at ∼1300 Hz. We checked the dependence of this feature, which at first glance looked like a numerical artifact, on simulation parameters, such as time step, resolution, and output frequency. For example, the dependence on the GW signal on the grid resolution for modelM10_LS220_no_manybody is given in AppendixA. The “gap” persisted at exactly the same location for all combinations of numerical parameters we considered for a given model, varying only slightly between the different models. We do not exclude the possibility that the “gap” is physical and attempt to explain it in part with the trappedg-mode of the PNS inner core in the next subsection.

3.2. Analytical Explanation of the Key Features of GW Signal from theM10_SFHo Model

In this subsection we focus on explaining the dominant features of the GW spectrogram, using modelM10_SFHo as an example (shown in Figure1), by means of a linear perturbation analysis.

The system of equations we solve combines the linearized equations of general relativistic hydrodynamics in a spherically symmetric conformally flat background metric (Banyuls et al.1997; Torres-Forné et al.2018) together with the Poisson equation. It can be summarized in the form

and

Here scalar functionsη_r = η_r(r) and represent the amplitudes of the decomposition of radial (ξ^r) and polar (ξ^θ) Lagrangian displacements of a fluid element with respect to its equilibrium position in terms of spherical harmonics:

whereσ is the mode frequency. The scalar function is the amplitude of the lapse function perturbation

and we define. The conformal factor of the metric,ψ, is equal to 1 in our numerical setup. The details of the derivation of Equations (8)–(11) are given in AppendixB. In the limitδα = 0 (the Cowling approximation), Equations (8)–(9) coincide with Equations (31)–(32) of Torres-Forné et al. (2018).

In Equations (8)–(11),P is the pressure,ρ is the rest-mass density of the matter,h is the specific enthalpy,c_s is the relativistic speed of sound, Γ₁ is the adiabatic index, is the radial component of the gravitational acceleration,, is the relativistic Brunt–Väisälä frequency, which in our case is equal to (see also Müller et al.2013)

and is the relativistic Lamb frequency

These quantities describe the spherically symmetric equilibrium background configuration, which we find by averaging the hydrodynamical output of our 2D simulations over polar angle.⁴

Figure4 shows the Brunt–Väisälä frequency for the averaged profile of theM10_SFHo model (middle panel). Black lines show the radial coordinates where the density is equal to 5.0 × 10⁹ g cm⁻³, 10¹⁰ g cm⁻³, and 10¹¹ g cm⁻³, with the latter density surface commonly used as a definition of the PNS boundary. Colored gray are the regions where the Brunt–Väisälä frequency is imaginary (), which means that they are convectively unstable. To further emphasize the convection, in the bottom panel of Figure4 we plot the anisotopic velocity defined as (Takiwaki et al.2012; Pan et al.2017)

where denotes spherical averaging. From this plot, one can see a convective layer inside the PNS, between ∼10 and ∼20 km, but the convective velocities there are much smaller than the ones above the PNS surface. The outer convective zone between the PNS boundary radius and the shock radius (shown in the top panel of Figure4) is recognized as the main driving region for the GW signal (Murphy et al.2009). The middle panel of Figure4 shows the imaginary Brunt–Väisälä frequency in the center of the PNS; however, the bottom panel shows no convection in that region.

Zoom InZoom OutReset image size

To solve Equations (8)–(11), we place the outer boundary condition at the radial coordinate whereρ = 10¹⁰ g cm⁻³ (solid black line in Figure4). There, we impose the condition ΔP = 0 on the Lagrangian perturbation of the pressure, which physically corresponds to a free surface of the PNS (see, e.g., Reisenegger & Goldreich1992). Mathematically, this boundary condition can be written as

Our treatment of the outer boundary condition is, therefore, different from the one in Torres-Forné et al. (2018), whereη_r = 0 at the shock position is imposed instead. At the innermost point, we impose a small radial displacement, use the regularity condition (Reisenegger & Goldreich1992)

and assume. As in Torres-Forné et al. (2018), we apply a trapezoidal rule to discretize the radial derivatives in Equations (8)–(11). Starting from the innermost point, we integrate the equations outward, inverting a 4 × 4 matrix of coefficients at every step. We use the bisection method to find the solutions satisfying Equation (17) at the outer boundary. The frequenciesσ/2π corresponding to these solutions are the eigenfrequencies of our model.

Figure5 shows the eigenfrequencies ofl = 2 (quadrupolar) modes overplotted on the GW spectrogram for modelM10_SFHo. Each eigenfrequency is represented by a number of nodes in the corresponding mode, i.e., the number of times the radial displacement functionη_r changes its sign along the radial coordinate. To avoid crowding the numbers, we show only the modes with the number of nodes <7 above the frequency 700 Hz and <4 below that frequency. Since the GW signal itself was obtained from the numerical simulations using a quadrupole formula (Finn & Evans1990), we primarily focus onl = 2 modes in this study. At the same time, we cannot exclude the case of nonlinear coupling between thel = 2 modes and the modes of differentl values, which can explain certain features of the GW signal (see, e.g., Torres-Forné et al.2018). For the interested reader, thel = 3 andl = 4 modes are shown in AppendixC.

Zoom InZoom OutReset image size

The left panel of Figure5 shows the results for the modes obtained under the Cowling approximation ( andf = 0 in Equations (8)–(11)). Starting from ∼0.4 s after bounce, a fundamental mode (thef-mode, with zero radial nodes) can be clearly identified. Above this mode, one can seep-modes (acoustic), for which the frequencies increase with the number of nodes, while below it there areg-modes, for which the frequencies decrease with the number of nodes. Before ∼0.4 s after bounce, as in Torres-Forné et al. (2018), we see the mixing and crossing between the different modes, during which they change the number of nodes.

The right panel of Figure5 shows the full solution of Equations (8)–(11), for. As in the left panel, the fundamental (f) mode clearly stands out after ∼0.4 s post-bounce time, but in this case it agrees very well with the strongest component of the GW radiation. This result is expected from physical grounds, and it shows that the Cowling approximation can indeed affect the analysis and should be used with caution when interpreting the GW signal from numerical simulations of CCSNe. Thep- andg-modes can be also identified in the left panel of Figure5. Interestingly, the GW spectrogram from our simulations shows almost no power below thef-mode, suggesting that the higher-orderg-modes of the PNS are not excited. Aside from the possible SASI and neutrino signal, which are expected to operate at the frequencies ≲100 Hz (see Kuroda et al.2016; Andresen et al.2017), there is no other apparent mechanism that could fill this “excluded region” of the spectrogram.

Figure6 shows the dependence of the obtained results on the position of the outer boundary, placed at the radial coordinate where the density reaches a given value. We remind the reader that the three choices of boundary density correspond to the three black lines in Figure4, with the middle value,ρ = 10¹⁰ g cm⁻³, being our default choice. From Figure6, it is seen that our approach does not let us capture the outerp-mode frequencies very accurately, because the result is very sensitive to the position of the outer boundary. This is probably related to the fact that thep-modes represent the sound waves propagating between the PNS surface and the shock position, a region that is not taken into account in our analysis. This, however, does not affect the qualitative conclusion that the high-frequency noise on the GW spectrogram above the dominant feature is at least partially associated with these modes. Another possible source of this noise is the turbulent convection between the PNS and the shock front, which is chaotic and does not necessarily represent any simple eigenmode of the system.

Zoom InZoom OutReset image size

At the same time, the frequency of the fundamental mode in Figure6 is almost insensitive to the position of the outer boundary, and the low-orderg-modes depend weakly on it. Importantly, this shows that the dominant GW frequency is not just proportional to the Brunt–Väisälä frequency at the surface of the PNS, as was suggested in earlier work. Indeed, Figure4 shows that the three black lines corresponding to the different outer boundary locations pass through very different values of the Brunt–Väisälä frequency. The fact that the fundamental quadrupolar eigenfrequency in Figure6 is nearly independent of the position of the outer boundary tells us that the dominant frequency of the GW signal is defined by the entire structure of the PNS, rather than by its surface characteristics alone.

The left panel of Figure7 illustrates the time evolution of the radial eigenfunctionη_r for thel = 2 modes associated with the dominant frequency of the GW signal. The eigenfunctions are normalized to 1 and plotted as a function of radial coordinate from the innermost grid point up to the location of the outer boundary. In Figure7, they are shifted along they-axis according to the time after bounce at which they are calculated (the time is indicated on the left side of the panel and directed downward). As we already mentioned, starting from ∼400 ms after bounce and until the end of the simulation, the main signal is represented by thef-mode, which has the largest amplitude at the PNS boundary surface and gradually decreases toward the center. Before that, in the time interval between ∼200 and ∼400 ms, this mode is smoothly connected to ag-mode having two radial nodes (see also the left panel of Figure5). The right panel of Figure7 shows the energy density defined as (Torres-Forné et al.2018)

for the corresponding eigenfunctions of the left panel. The figure shows that the shape of the fundamental eigenfunction is very similar in the case of the Cowling approximation (black lines) and in the case when (red lines). The energy density of the modes shows less agreement. Note that the definition of contains the mass density, which is larger in the inner region than at the the surface of the PNS. Therefore, even a barely visible disagreement between the eigenfunctions in the inner region may lead to a large disagreement between the energy density distributions (see, e.g., the 0.48 s snapshot in Figure7).

Zoom InZoom OutReset image size

Finally, in Figure8 we attempt to address the nature of the “gap” seen in our spectrograms by performing the linear perturbation analysis of the PNS inner core only. For that, we place the outer boundary at the inner maximum of the Brunt–Väisälä frequency, which roughly corresponds to the radial coordinate of 10 km (see Figure4), and solve the system of Equations (8)–(11) in that inner region, using boundary condition (17). This approach is not strictly accurate, but it gives us an idea about the eigenfrequencies of the inner core. In Figure8, the digits show the number of radial nodes in the corresponding modes. All modes with the number of nodes less than 7 are shown, without selection. The resolution of the PNS core in our simulations does not exceed a few tens of grid points, which leads to the spurious nodes and the step-like behavior of the eigenfrequencies (for the same reason, it does not make sense to plot the modes with a larger number of nodes). Nevertheless, a part of the core eigenfrequencies lies very close to the “gap” position in the spectrogram and roughly reproduces its morphology.

Zoom InZoom OutReset image size

We speculate that the “gap” may be the result of interaction between the trapped PNS core mode and the other (p- orf-) modes of the system, probably by means of an avoided crossing (Christensen-Dalsgaard1981; Yoshida & Eriguchi2001; Stergioulas2003). One of the simplest examples of the avoided crossing phenomenon is the case of two coupled classical oscillators, where the eigenfrequencies demonstrate characteristic splitting in the strong coupling regime (Novotny2010). Analogously, one may view the PNS as a coupled system of the inner core and the outer convectively stable shell, mediated by the inner PNS convection region (Dessart et al.2006). In this picture, the modes of the inner core may repel the modes of the shell, leading to an empty region in the frequency space with the width related to the strength of the coupling between the two. At the same time, the inner mode itself is most likely not excited, because it is shielded from the downfalling plumes of the postshock convection region by the PNS surface (however, the inner PNS convection itself may be a source of mode excitation; see, e.g., Müller et al.2013). To clarify the nature of the “gap,” higher-resolution simulations are necessary. If this spectrogram feature is real, it could serve as an interesting analysis tool to probe the structure of the inner PNS.

3.3. Dependence of the GW Signal on Parameters

In this section, we outline the key dependences of the GW signal on the progenitor mass and rotational angular velocity, on the EOS, and on the details of the microphysics, such as the inclusion of the many-body corrections to the neutrino–nucleon scattering rate and the implementation of the gravity solver. While qualitatively the GW signals from all our models are very similar, the frequency of the dominant feature is sensitive to the EOS and the neutrino–nucleon opacities and almost insensitive to the progenitor mass. We apply the analysis of Section3.2 to all nonrotating models from our set, and we confirm the association between the dominant GW feature and the fundamentall = 2 mode in each case.

3.3.1. Dependence of the GW Signal on the Progenitor Mass

Figure9 shows the GW spectrograms and waveforms for modelsM10_SFHo,M13_SFHo, andM19_SFHo, which are simulated with the identical numerical setup and differ only in the progenitor mass. Two of the models,M10_SFHo andM19_SFHo, explode at ∼400 and ∼350 ms after bounce, respectively, and have a characteristic explosion “tail” in their waveforms (Murphy et al.2009; Yakunin et al.2010; Müller et al.2013). ModelM13_SFHo does not explode. White markers on the spectrograms indicate the eigenfrequencies of the fundamentall = 2 modes, found as described in Section3.2 for each model. In general, we see good agreement between the analytical eigenfrequencies and the dominant GW signal, with the largest deviation seen in the post-explosion phase of the 19M_⊙ model.

Zoom InZoom OutReset image size

For comparison, the red lines in the top panel of Figure9 show the peak GW frequencyf_peak computed as suggested in Murphy et al. (2009) and Müller et al. (2013), where it is associated with the surface value of the Brunt–Väisälä frequency divided by 2π. In Figure9, we use Equation (17) of Müller et al. (2013), omitting the factor, whereM_PNS andR_PNS are the mass and radius of the PNS, respectively. As in Pan et al. (2017), we find that removing this factor results in better agreement betweenf_peak and the dominant feature of the GW spectrogram.⁵ The three lines correspond to the three different density isosurfaces ofρ = 5.0 × 10⁹ g cm⁻³, 10¹⁰ g cm⁻³, and 10¹¹ g cm⁻³, which can represent the PNS surface. The plot shows that using the conventional definition for the PNS surface,ρ = 10¹¹ g cm⁻³, the analytical formula forf_peak may provide a good fit to the dominant GW signal. At the same time, the value off_peak is sensitive to the location of the PNS surface, which currently lacks strict physical definition.

The agreement between the analytic eigenfrequencies and the GW spectrograms allows us to compare the spectrograms by comparing the frequencies. Figure10 shows the full results of the linear perturbation analysis for the three considered models, performed as in Section3.2. In the bottom panel of Figure10, the filled symbols show the modes with larger than zero number of nodes, while the open symbols show thef-mode frequencies. This plot demonstrates that the PNS eigenfrequencies in general, and the frequencies of the fundamental quadrupolar mode in particular, are strikingly similar between the models, despite the large difference in their progenitor masses and even in the waveforms themselves.

Zoom InZoom OutReset image size

This may reflect the fact that the evolution of the PNS radius is very similar between the models with different progenitor masses, which is shown in the top panel of Figure10 and was already noticed in the literature for a wide range of progenitors differing not only by the ZAMS mass but also in the metallicity (see, e.g., Figure 7 of Bruenn et al.2016; Figure 10 of Summa et al.2016; Figure 15 of Radice et al.2017).⁶ Indeed, if the GW signal from CCSNe is so tightly related to the PNS eigenmodes, the structure of the PNS should be the main factor defining the time–frequency structure of this signal.

The same argument cannot be applied to the amplitude of the GW signal, which, instead, must depend on the mechanism of excitation of the PNS modes. It was shown in many previous studies that the GW signal from CCSNe experiences sudden increases in amplitude at the moments when the PNS surface is hit by the downfalling accretion “plumes” (Murphy et al.2009; Müller et al.2013; Yakunin et al.2015). It is, therefore, natural to expect that the GW power will depend on the details of the postshock accretion, which takes place above the PNS surface and is largely determined by the core structure of the progenitor. Figure11 shows the energy emitted in GWs due to the matter motions alone as a function of time for modelsM10_SFHo,M13_SFHo, andM19_SFHo. In these models, we do not see a monotonic dependence of the GW power on the progenitor mass, with modelM13_SFHo producing the weakest signal among the three. In fact, it is hard to expect such a monotonic dependence, because the dependence of the progenitor core structure itself on the progenitor ZAMS mass is not monotonic (Sukhbold et al.2016) and, moreover, may be intrinsically chaotic (Sukhbold et al.2017). For this reason, we advise using caution when deducing the dominant signal frequency based on the total GW energy spectrum, especially if it is done for the purpose of comparing models with different progenitor masses. Accretion downflows hitting the PNS surface at random moments of time may give more weight to the system eigenfrequencies in those moments, complicating the overall picture. Instead, the comparison of the time–frequency spectrograms serves this purpose best.

Zoom InZoom OutReset image size

The strongest signal among all our models is produced by modelM19_SFHo. In Figure12, we present the linear 3D representation of the GW spectrogram from this model. This figure emphasizes the point made in Section3.1, that the GW signal may stay strong for a long time after the explosion (more than a second in the case ofM19_SFHo). The large offset from zero seen in the GW strain of this model at late times (bottom panel of Figure9) suggests a very asymmetric character for its explosion. This is indeed the case, as demonstrated in Figure 6 of Vartanyan et al. (2018), which shows snapshots of the electron fraction and entropy of this model at different moments of time. For the analogous snapshots of the 10M_⊙ model, we refer the reader to Radice et al. (2017).

Zoom InZoom OutReset image size

The dependence of the CCSN GW signal on the progenitor mass was previously studied in a number of works (Murphy et al.2009; Müller et al.2013; Yakunin et al.2015). For example, Müller et al. (2013) report ∼30% differences in the typical emission frequencies between their 11.2 and 25M_⊙ models, which they admit to be small for this large a mass difference. Our comparison, however, shows even smaller scatter, no more than ∼5%–10% in frequency across the considered mass range, without any systematic trend. On one hand, we cannot exclude that at least part of the difference between the 11.2 and 25M_⊙ models of Müller et al. (2013) may come from the fact that they were simulated with a slightly different EOS (we discuss the dependence of the signal on EOS in the next subsection). On the other hand, the models of Müller et al. (2013) treat general relativity more accurately by solving the relativistic equations of hydrodynamics in the conformally flat approximation, while our work uses the effective potential approach, which may also affect the dominant frequency of the signal (see Müller et al.2013).

3.3.2. Dependence of the GW Signal on the Equation of State

Figure13 shows the GW spectrograms and waveforms for modelsM10_LS220,M10_SFHo, andM10_DD2, which were simulated with three different EOSs. All other numerical parameters and the details of microphysics are the same between these models. We find that the EOS has a large impact on the amplitude of the GW signal, its dominant frequency, the total emitted energyE_GW (see Table1), and even the qualitative outcome of the simulation (modelM10_SFHo explodes, unlike the other two). Similar EOS sensitivity of the simulation outcome was recently reported by Pan et al. (2017) in the context of the GW signal from black hole formation in failed SNe.

Zoom InZoom OutReset image size

In order to emphasize the dependence of the dominant GW frequency on the EOS, we overplot the GW spectrograms of these models in the bottom panel of Figure14. Open markers represent thef-mode eigenfrequencies found from the linear perturbation analysis of the models, as described in Section3.2. The top panel of Figure14 shows the evolution of the PNS radii taken at a density of 10¹⁰ g cm⁻³. Compared to the top panel of Figure10, the difference between the PNS radii in Figure14 is slightly larger and more systematic, which translates into the systematic ∼10%–15% difference in the dominant frequencies of the GW signal, which, in turn, is well captured by our analysis (the largest disagreement is seen in theM10_DD2 model). Interestingly, among the three EOSs used in our study, SFHo is the “softest” one, while DD2 is the “hardest.” Nevertheless, the smallest PNS radius and the largest GW frequency are produced by the LS220 EOS. This suggests that the EOS dependence of the GW signal, as well as the overall core evolution, may not necessarily be described in terms of a single stiffness parameter defined at zero temperature.

Zoom InZoom OutReset image size

We quantify the dependence of thef-mode eigenfrequency on time for modelsM10_LS220,M10_SFHo, andM10_DD2 by fitting it with a polynomial. We find that a simple quadratic function in the form, wheref is frequency in Hz,t is time in seconds, andA,B, andC are coefficients, adequately describes the dependence over the first ≳1.5 s after bounce, while the core keeps shrinking. Eventually, the PNS will cool down and deleptonize, which could lead to the flattening of the frequency–time curve. The quadratic fits are shown with the lines of corresponding color in Figure14 and explicitly written down in the lower right corner of the figure. These fits can be used as priors when looking for the CCSN GW signal in the data from ground-based laser interferometers, such as LIGO, Virgo, or KAGRA. At the same time, we emphasize that the accuracy of these fits may be affected by the details of the physics and microphysics used in our (and other) codes. For example, to demonstrate the sensitivity of the GW signal to the details of the neutrino opacity, we compare the spectrograms from modelsM10_LS220 andM10_LS220_no_manybody in Figure15. The many-body corrections to the neutrino–nucleon scattering cross section decrease the neutrino opacity, which leads to the faster contraction of the PNS, as shown in the top panel of Figure15. The bottom panel of Figure15 shows that neglecting these corrections results in a ∼10% shift in the dominant GW frequency. Another factor influencing the GW frequency is the description of the gravitational field (see, e.g., Müller et al.2013). Taking all these factors into account, we expect the accuracy of the fits from Figure14 to be not worse than ∼30%.

Zoom InZoom OutReset image size

On the other hand, the power of the GW signal does demonstrate monotonic dependence on the stiffness of the EOS, with the hardest EOS (DD2 in our case) producing the weakest signal. As we already mentioned in the previous subsection, the amplitude of the GW signal is largely determined by accretion and postshock convection, which act as driving forces for the excitation of the PNS oscillations. It was found in previous work (Marek et al.2009; Kuroda et al.2016,2017) that softer EOSs result in more vigorous SASI activity. While we do not clearly identify SASI in any of the three models, we also find that the shock oscillations are strongest in theM10_SFHo model and weakest in theM10_DD2 model. This leads to stronger excitation of the PNS modes and more powerful GW signals in case of the softest EOS.

3.3.3. Dependence of the GW Signal on Rotation

Simulations of rotating core collapse were the first to predict and study the GW emission from CCSNe (Dimmelmeier et al.2007,2008; Ott et al.2007,2012; Abdikamalov et al.2010). Because of the symmetry breaking introduced by rotation, these models produce strong GW signals already at the early stages of collapse and bounce, which makes even short (few tens of milliseconds) simulations very informative. Not very demanding in terms of the neutrino physics, these simulations progressed enough to establish the connection between the properties of the GW signal and the progenitor core parameters (Summerscales et al.2008; Röver et al.2009; Logue et al.2012; Abdikamalov et al.2014; Engels et al.2014; Fuller et al.2015; Powell et al.2016; Richers et al.2017). The main limitation of these papers is that fast rotating cores are not very common among CCSN progenitors (Heger et al.2005; Woosley & Heger2006). Here we focus on the GW signal from a moderately (Ω₀ = 0.2 rad s⁻¹) rotating progenitor and follow it for a full second after bounce, which, to the best of our knowledge, is currently the longest simulation of its kind, for which the GW signal has been extracted.

In the rotating model, we use the multipole gravity solver of Müller & Steinmetz (1995). For all other models shown before, we used a monopole approximation for the gravitational potential (Marek et al.2006). As an aside, to show how the gravity implementation alone influences the GW signal, we compare the spectrograms and waveforms from modelsM13_SFHo andM13_SFHo_multipole in Figure16. The difference between the models is noticeable, although not large, resulting in ∼10% shift in the dominant frequency by the end of theM13_SFHo_multipole simulation. This tells us that the full general relativistic approach to gravity (which is, strictly speaking, the only correct approach) is important for the accurate quantitative description of the GW signal.

Zoom InZoom OutReset image size

Figure17 shows the GW spectrograms and the waveforms from the nonrotating (M10_SFHo_multipole) and rotating (M10_SFHo_rotating) models, which have identical numerical setups, apart from the angular velocity. In agreement with the previous literature, the rotating model generates a strong GW signal at the core bounce, lasting for a few tens of milliseconds, followed by the short quiescent phase. At the same time, the main component of the GW signal is noticeably weaker for this model, though the dominant frequency does not seem to change much. Interestingly, the “gap” still persists in the GW spectrograms of both models.

Zoom InZoom OutReset image size

The main findings of our study can be briefly summarized as follows:

• 1.  
  We reproduce the dominant, long-lasting GW signal from CCSNe by means of linear perturbation analysis and associate it with ag-mode having two radial nodes at the early stage (∼200–400 ms after bounce) and with thef-mode later on (from ∼400 ms until more than a second after bounce). This finding presages future opportunities for the analytical study of the CCSN GW signal.
• 2.  
  We demonstrate a weak dependence of the dominant GW frequency on the progenitor ZAMS mass and provide a simple quadratic fit for it as a function of time for three different EOSs. This may help identify possible CCSN candidates in the GW data from ground-based laser interferometers.
• 3.  
  We identify a new feature in the GW spectrogram, which looks like a “gap” across the noisy GW emission in the first ∼200–700 ms after bounce. Our attempts to explain it as a numerical artifact failed. We explain the “gap” as the interaction between the outerp-modes andg-modes of the PNS inner (∼10 km) core, probably as a result of avoided crossing.
• 4.  
  We show the effect of moderate (0.2 rad s⁻¹) initial progenitor rotation on the GW signal. The rotation strengthens the bounce signal but weakens the dominant part of the post-bounce GW emission.

All simulations analyzed in our study are 2D, which raises the question of how our conclusions will change in the full 3D case. It is known from previous studies that the success of an explosion in the CCSN simulations largely depends on the hydrodynamical instabilities and the associated turbulent pressure behind the stalled shock (Burrows et al.1995; Murphy et al.2013; Couch & Ott2015; Müller & Janka2015; Abdikamalov et al.2016; Takahashi et al.2016; Müller et al.2017), which also increases the exposure of matter to neutrino heating (Buras et al.2006; Murphy & Burrows2008). However, it is known that turbulence has different properties in 2D and 3D (Kraichnan1967), and it has been shown that this difference artificially facilitates explosion (Hanke et al.2012; Dolence et al.2013; Couch & O’Connor2014; Takiwaki et al.2014; Abdikamalov et al.2015; Couch & Ott2015). Therefore, if the properties of turbulence in the gain region were directly reflected in the GW spectrogram, we would expect it to differ in 3D. Instead, our analysis suggests that the strongest component of the GW signal is associated with the fundamental mode of the PNS itself, which is expected to be nearly spherical even in the 3D case. Turbulence in this case acts only as a driving force exciting the mode oscillations. This makes us believe that the time–frequency structure of the GW signal shown here and its linear analysis will still be applicable for 3D models, while the amplitude may change (become smaller). The same was recently suggested in Yakunin et al. (2017), where the authors obtained similar behavior of the GW signal for a 2D and a 3D model. In their 3D case, convection was characterized by a larger number of relatively small scale structures, as opposed to the few massive accretion funnels in 2D. This led in 3D to smoother GW energy emission but caused only moderate changes in its spectral distribution, vis-à-vis their 2D results, during the first 450 ms of the signal. More about the comparison between the 2D and 3D GW signals from CCSNe may be found in Andresen et al. (2017).

Interestingly, our linear analysis presents an opportunity to predict the dominant frequency of the GWs from CCSNe based on 1D simulations. However, this approach should be applied with great caution, because, for example, the evolution of the PNS radius differs between the 1D and 2D simulations for the same models (Radice et al.2017). At the same time, such an analysis allows one to quickly cover large regions of parameter space related to the EOS and microphysics, in order to investigate which of the parameters has the strongest influence on the GW signal.

Pan et al. (2017) suggested that increasing the sensitivity of the next-generation GW detectors in the ∼1000 Hz window is very important for the study of the BH formation in failed SNe (see also Kuroda et al.2018). We add to this statement that high-frequency sensitivity is crucial for the detection of the GW signal from the successful SN explosions as well. Increasing the sensitivity of aLIGO and KAGRA in this band would help us to fully exploit the luck of the next nearby SN discovery and trace the high-frequency GW signal of a newborn NS.

We thank Aaron Skinner and James Stone for helpful discussions and feedback. We thank Pablo Cedrá-Durán and José Antonio Font for finding an error in the original calculations, and for other helpful suggestions. The authors would like to acknowledge support of the U.S. NSF under award AST-1714267, the Max-Planck/Princeton Center (MPPC) for Plasma Physics (under award NSF PHY-1144374), and the DOE SciDAC4 grant DE-SC0018297 (under subaward 00009650). The authors employed computational resources provided by the TIGRESS high-performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Princeton University Office of Information Technology. They also acknowledge a supercomputer allocation by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy (DOE) under contract DE-AC03-76SF00098. D.R. acknowledges support from a Frank and Peggy Taplin Membership at the Institute for Advanced Study and the Max-Planck/Princeton Center (MPPC) for Plasma Physics (NSF PHY-1523261).

Figure18 shows the GW spectrograms and waveforms of modelM10_LS220_no_manybody for three different grid resolutions, of which the lowest (“standard”) is used in all other models of this study. It is clear from the figure that the overall structure of the GW signal and its spectrogram depend weakly on resolution, demonstrating the robustness of our results.

Zoom InZoom OutReset image size

Here we derive the set of equations describing linear perturbations of a spherically symmetric background, including the perturbation of the lapse function. Generally, our calculations follow the same scheme as described in Torres-Forné et al. (2018), and for simplicity of comparison we use the same notation where possible. All equations are given in geometrized units. We start with a static spherically symmetric conformally flat spacetime metric in isotropic coordinates:

whereα is the lapse function,ψ is the conformal factor, and is the flat spatial 3-metric. In this metric, the equations of general relativistic hydrodynamics for a perfect fluid can be rendered in the form (Banyuls et al.1997; Torres-Forné et al.2018)

Here is the energy-momentum tensor of a perfect fluid, whereρ is its rest-mass density,P is the pressure,u^μ is the 4-velocity,h ≡ 1 +  + P/ρ is the specific enthalpy, and is the specific internal energy. denotes the Christoffel symbols, and is the determinant of the 3-metric,. The conserved rest-mass densityD, momentum density in thej-directionS_j, and total energy densityE are defined as

where is the Lorentz factor, andνⁱ and represent the Eulerian and “advective” velocities, in the spherically symmetric case equal touⁱ/W and, respectively.

As in Torres-Forné et al. (2018), we consider the linear perturbations of the system with respect to the equilibrium static background, for which the only nonzero radial component of Equation (22) is

where is the radial component of the gravitational acceleration. At the same time, in addition to the perturbation of density, pressure, and velocity, we introduce the nonzero perturbation of the lapse function,α. This addition does not fully relax the Cowling approximation, but it closely mimics the conditions of our numerical simulations, where the shift vectorβⁱ = 0 and the conformal factor is fixed atψ = 1. Following Torres-Forné et al. (2018), we denote the Eulerian perturbations of the quantities byδ and the Lagrangian perturbations by Δ, where the relation between the two for any quantity, e.g.,ρ, is

Hereξⁱ is the Lagrangian displacement of a fluid element, related to the advective velocity as

After perturbing the quantities by substituting, e.g.,, and leaving only terms of linear order, Equations (21) and (22) can be rewritten as

with the three components of Equation (29) taking the form

The condition of adiabaticity of the perturbations

wherec_s is the relativistic speed of sound and Γ₁ is the adiabatic index, allows one to write (Torres-Forné et al.2018)

where

is the relativistic version of the Schwarzschild discriminant ande ≡ ρ(1+). For a spherically symmetric background, the only nonzero component of is. Due to the adiabatic nature of perturbations, Equation (23) does not add any information.

To close the system of Equations (21)–(22), we use the Poisson equation

where Φ is the gravitational potential. Using the relationα = e^Φ, we rewrite it as

Following Torres-Forné et al. (2018), we consider only polar perturbations and expand them in terms of spherical harmonics as

where,,η_r, and are scalar functions depending only on radial coordinate. With thisAnsatz, and using adiabaticity condition (33), Equation (37) may be brought to the form

To conveniently find the numerical solution, we introduce the function and break this second-order equation into two first-order equations:

Equation (31) results in

where, after Torres-Forné et al. (2018), we have defined. Using Equations (33), (34), and (42) in Equations (28) and (30), we get

where is the relativistic Brunt–Väisälä frequency defined as

and is the relativistic Lamb frequency

Finally, using Equation (42), we bring Equations (40) and (41) to the form

To find the eigenfrequencies of the linear perturbation modes,f = σ/(2π), we numerically solve the system of first-order differential Equations (43), (44), (47), and (48).

Figure19 shows the results of the linear perturbation analysis of the M10_SFHo model, analogous to the one described in Section3.2, for the l = 3 and l = 4 modes.

Zoom InZoom OutReset image size