A MODEL FOR GRAVITATIONAL WAVE EMISSION FROM NEUTRINO-DRIVEN CORE-COLLAPSE SUPERNOVAE

For the 2D simulations presented in this paper, we use BETHE-hydro (Murphy & Burrows 2008a, 2008b), an Arbitrary Lagrangian–Eulerian (ALE) hydrodynamics code. The basic equations of hydrodynamics are the conservation of mass, momentum, and energy:

$$\begin{equation} \frac{d \rho }{d t} = - \rho \nabla \cdot \bf {v}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \rho \frac{d \bf {v}}{d t} = - \rho \nabla \Phi - \nabla P, \end{equation} \tag{ 2 }$$

and

$$\begin{equation} \rho \frac{d \varepsilon }{d t} = - P \nabla \cdot \bf {v} + \rho (\mathcal {H} - \mathcal {C}), \end{equation} \tag{ 3 }$$

where ρ is the mass density, v is the fluid velocity, Φ is the gravitational potential, P is the isotropic pressure, ε is the specific internal energy, and $d/dt = \partial / \partial t + \bf {v} \cdot \nabla$ is the Lagrangian time derivative. In this work, the neutrino heating, $\mathcal {H}$, and cooling, $\mathcal {C}$, terms in Equation (3) are assumed to be

$$\begin{equation} \mathcal {H} = 1.544 \times 10^{20} L_{\nu _e} \left(\frac{100\,{\rm km}}{r} \right)^2 \left(\frac{T_{\nu _e}}{4 {\rm \, MeV}} \right)^2 \left[ \frac{\rm erg}{{\rm g} \, {\rm s}} \right], \end{equation} \tag{ 4 }$$

in which we assume the free-streaming limit,⁷ and

$$\begin{equation} \mathcal {C} = 1.399 \times 10^{20} \left(\frac{T}{2 {\rm \, MeV}} \right)^6 \left[ \frac{\rm erg}{{\rm g} \, {\rm s}} \right] \,. \end{equation} \tag{ 5 }$$

Note that these approximations for heating and cooling by neutrinos (Bethe & Wilson 1985; Janka 2001) depend upon local quantities and predefined parameters. They are ρ, temperature (T), the distance from the center (r), the electron–neutrino temperature ($T_{\nu _ e}$), and the electron–neutrino luminosity ($L_{\nu _e}$), which we give in units of 10⁵² erg s⁻¹. By using Equations (4) and (5), we gain considerable time savings by approximating the effects of detailed neutrino transport. For all simulations, we set $T_{\nu _e} = 4$ MeV. In Equation (4), it has been assumed that $L_{\nu _e} = L_{\bar{\nu }_e}$ and that the mass fractions of protons and neutrons sum to one. Therefore, the sum of the electron- and anti-electron–neutrino luminosities is $L_{\nu _e \bar{\nu }_e} = 2 L_{\nu _e}$. Closure for Equations (1)–(3) is obtained with the finite-temperature nuclear EOS of Shen et al. (1998) to which we have added the effects of photons, electrons, and positrons. As such, the EOS has the following dependencies:

$$\begin{equation} P = P(\rho,\varepsilon,Y_e), \end{equation} \tag{ 6 }$$

where Y_e is the electron fraction. Given this Y_e dependence, we also solve the equation:

$$\begin{equation} \frac{d Y_e}{dt} = \Gamma _e, \end{equation} \tag{ 7 }$$

where Γ_e is the net rate of Y_e change.

The heating and cooling functions, Equations (4) and (5), are valid for optically thin regions only. To mimic the reduction in heating and cooling at depth where the PNS is opaque to neutrinos, we multiply by $e^{-\tau _{\rm eff}}$. In 1D simulations, the effective optical depth for electron- and anti-electron-type neutrinos is calculated by

$$\begin{equation} \tau _{\rm eff} = \int ^{\infty }_r \kappa _{\rm eff}(r) dr, \end{equation} \tag{ 8 }$$

where the effective opacity, κ_eff is given by Equation (15) or (16) of Janka (2001), which we reproduce without derivation:

$$\begin{equation} \kappa _{\rm eff} (r) \approx 1.5 \times 10^{-7} X_{n,p} \left(\frac{\rho }{10^{10} \, {\rm g} \, {\rm cm}^{-3}} \right) \left(\frac{T_{\nu _e}}{4 {\rm \, MeV}} \right)^2 {\rm cm}^{-1}, \end{equation} \tag{ 9 }$$

where X_n,p, the composition weighting, is ${\sim} \frac{1}{2}(Y_n + Y_p)$ and Y_n and Y_p are the number fractions of free neutrons and protons. Because the radial density profile after bounce is monotonic and, to good approximation, a series of power laws, τ_eff is roughly a function of local density only. We use our 1D simulations to calibrate this τ–ρ relationship and parameterize τ as function of ρ for the 2D simulations.

The heating and cooling terms are set to zero during collapse and are turned on after bounce. We have experimented with several methods to do this, but find that they give similar results provided that these terms "turn on" within a few milliseconds of bounce. For convenience, we include heating and cooling once the maximum value of τ is greater than 100.

Using BETHE-hydro (Murphy & Burrows 2008a), we solve Equations (1)–(3) in two dimensions by the ALE method. To advance the discrete equations of hydrodynamics by one time step, ALE methods generally use two operations, a Lagrangian hydrodynamic step followed by a remap. The structure of BETHE-hydro's hydrodynamic solver is designed for arbitrary-unstructured grids, and the remapping component offers control of the time evolution of the grid. Taken together, these features enable the use of time-dependent arbitrary grids to avoid some unwanted features of traditional grids.

For the calculations presented in this paper, we use this flexibility to avoid the coordinate singularity of spherical grids in two dimensions. While spherical grids are generally useful for core-collapse simulations, the convergence of grid lines near the center places extreme constraints on the time step via the Courant–Friedrichs–Levy condition. A common remedy is to simulate the inner ∼10 km in 1D or to use an inner boundary condition. Another approach, which has been used in VULCAN/2D simulations (Livne 1993; Livne et al. 2004; Burrows et al. 2007c), is to avoid the singularity with a grid that is pseudo-Cartesian near the center and smoothly transitions to a spherical grid at larger radii.

In this paper, we use a similar grid, the butterfly mesh (Murphy & Burrows 2008a, 2008b) interior to 34 km and a spherical grid exterior to this radius. See Murphy & Burrows (2008a) for an example. For all 180° simulations, the innermost pseudo-rectangular region is 50 by 100 zones, and the region that transitions from Cartesian to spherical geometry has 50 radial zones and 200 angular zones. The outermost spherical region has 200 angular zones and 500 radial zones that are logarithmically spaced between 34 km and 8000 km in most cases. Due to EOS constraints, the outer boundary goes out to 6000 km for the 12 M_☉ model. The butterfly portion has an effective radial resolution of 0.34 km with the shortest cell edge being 0.28 km. In 2D simulations, we do not remap every time step. Since the time step is limited by the large sound speeds in the PNS core, we can afford to perform several Lagrangian hydrodynamic solves before remapping back to the original grid. By remapping rarely, we gain considerable savings in computational time.

We make several approximations to expedite the calculations. For one, gravity is calculated via $\vec{g} = -{\it G M}_{\rm int}/r^2$, where M_int is the mass interior to the radius r. Furthermore, since we do not solve the neutrino transport equations, we do not obtain electron/positron capture and emission rates, which are necessary for self-consistent Y_e profiles and evolution. Instead, we use a Y_e versus ρ parameterization that is quantitatively accurate during collapse and gives the correct qualitative trends after bounce. Liebendörfer (2005) observed that 1D simulations including neutrino transport produce Y_e values during collapse that are essentially a function of density alone. This allows for a parameterization of Y_e as a function of ρ. To change Y_e, we use results of 1D SESAME (Burrows et al. 2000; Thompson et al. 2003) simulations to define the function Y_e(ρ), and we employ the prescription of Liebendörfer (2005) to calculate local values of Γ_e. Though this prescription is most appropriate during collapse, we also apply it after bounce to continue deleptonizing postshock material at later times. While this is not entirely consistent with our heating and cooling prescription, the most important effects of deleptonization are included without the need for expensive, detailed neutrino transport.

The approximations for neutrino–matter interactions, electron capture rates, and gravity give simulations that reproduce the primary features of the core-collapse problem. These include a PNS, a stalled shock at ∼200 km, a gain radius at ∼100 km (that divides the gain region above with net heating from the cooling region below), prompt convection, PNS convection from ∼20 to ∼45 km, postshock convection, and the SASI. While the Y_e–ρ parameterization is designed to give the proper Y_e profile during collapse, we do not reproduce the Y_e trough at later times that extends from ∼40 km to the shock (e.g., Dessart et al. 2006). Dessart et al. (2006) conclude that the negative lepton gradient from ∼15 to ∼30 km in simulations employing consistent neutrino transport drives the PNS convection. In our simulations, PNS convection is driven by a negative entropy gradient. Although 1D simulations that solve neutrino transport more consistently show a similar negative entropy gradient, the depth and width are smaller in those simulations (compare Figures 13 and 6 of Murphy & Burrows 2008b; Liebendörfer et al. 2005). Furthermore, our PNS convection velocities are roughly a factor of 2 larger than those obtained by Dessart et al. (2006). In summary, while we obtain a PNS convective region, our PNS convection is predominantly driven by an entropy gradient rather than by a negative lepton gradient and is more extended and vigorous compared to more realistic calculations. However, we show in later sections that though the GW signal due to PNS convection is important, it is not an overwhelmingly dominant feature of GW emission from CCSNe. Therefore, the differences in PNS convection should not change the qualitative conclusions of this paper.

As initial conditions, we use the 12, 15, 20, and 40 M_☉ core-collapse progenitor models of Woosley & Heger (2007), where these masses correspond to the zero-age-main-sequence (ZAMS) masses. In Figure 1, we plot the mass accretion rate versus time after bounce (solid black lines) for each of these models. In general, because the more massive progenitors have extended core structures, the accretion rates for them are higher at a given time. Since higher accretion rates require larger neutrino luminosities for successful explosions (Murphy & Burrows 2008b), this translates into higher critical neutrino luminosities for the more massive progenitors. In fact, if we assume that only the neutrino mechanism leads to explosion, the required neutrino luminosity for successful explosions of the 40 M_☉ progenitor is so high, ≳10⁵³ erg s⁻¹, we doubt that it will explode at all and suspect that the non-explosive simulation using $L_{\nu _e} = 6 \times 10^{52}$ erg s⁻¹ represents the likely outcome of core collapse for the 40 M_☉ progenitor. In general, these trends with progenitor mass will translate to trends in the GW signal, which we discuss in forthcoming sections.

Zoom In Zoom Out Reset image size

Although core structure follows a general trend with the progenitor's ZAMS mass, the correlation is known to be somewhat chaotic (Woosley et al. 2002). For example, a substantial drop in the 20 M_☉ model's accretion rate corresponding to the interface between the Si and O burning shells temporarily reduces the accretion below even the 12 and 15 M_☉ progenitors at ∼300 ms past bounce. This has consequences for the development of the SASI, explosion, and the GW signal.

We extract the GW signal from our hydrodynamic simulations via the slow-motion, weak-field formalism (e.g., Misner et al. 1973) and consider the dominant mass-quadrupole only. In this approximation, the GW field is directly proportional to the transverse-traceless (TT) part of the second time derivative of the reduced mass-quadrupole tensor $I\hspace{-5pt}{\scriptscriptstyle \stackrel{\hbox{--}}{}}_{jk}$. Specifically,

$$\begin{equation} h^{TT}_{jk} (t,{\bf x}) = \frac{2}{c^4}\frac{G}{D} \bigg [\frac{d^2}{dt^2}\, I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{jk} (t - D/c)\bigg ]^{TT}, \end{equation} \tag{ 10 }$$

and

$$\begin{equation} I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{jk} = \int \rho \bigg (x^j x^k - \frac{1}{3}\delta ^{jk} x_i x^i\bigg) d^3 x\,. \end{equation} \tag{ 11 }$$

In the above, D is the observer distance to the source of emission. All other variables and constants have their usual meaning and the Einstein sum convention is assumed. For numerical convenience, we employ the so-called first-moment-of-momentum-divergence (FMD) recast of Equation (10) proposed by Finn & Evans (1990), which employs the continuity equation to remove one time derivative. In axisymmetry, $I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{ij}$ reduces to one independent component ($I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{zz}$) and, adapted to spherical coordinates, the FMD formula for its first time derivative reads

$$\begin{eqnarray} \frac{d}{dt}\, I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{zz} &=& \frac{8\pi }{3} \int _{-1}^{1} d\cos \theta \int _{r_1}^{r_2} dr\, r^3 \rho \, \, \cdot \nonumber \\ &&\bigg [ P_2(\cos \theta)\, v_r + \frac{1}{2} \frac{\partial }{\partial \theta } P_2(\cos \theta)\, v_\theta \bigg ], \end{eqnarray} \tag{ 12 }$$

where P₂(cos θ) is the second Legendre polynomial in cos θ, and v_r and v_θ are the fluid velocities in the radial and lateral direction, respectively. The axisymmetric⁸ GW strain h₊ = h_θθ(dropping the TT superscript) is then given by

$$\begin{equation} h_+ = \frac{3}{2}\frac{G}{Dc^4} \sin ^2 \alpha \frac{d^2}{dt^2}\,\, I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{zz}, \end{equation} \tag{ 13 }$$

where α is the angle between the symmetry axis and the line of sight of the observer and the second time derivative is obtained via straightforward differentiation.⁹ For our simulations, this approach yields good results and we do not find it necessary to employ the modification of Equation (10) proposed by Blanchet et al. (1990) that removes both time derivatives, but introduces a sensitive dependence on derivatives of the Newtonian gravitational potential, which we treat only in spherical fashion.

In our analysis of the GW signature, we not only consider the dimensionless GW strain h₊, but also compute the total emitted energy in GWs, given by

$$\begin{equation} E_\mathrm{GW} = \frac{3}{10} \frac{G}{c^5} \int _0^t \bigg (\frac{d^3}{dt^3}\, I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{zz}\bigg)^2\,dt, \end{equation} \tag{ 14 }$$

and the GW spectral energy density via

$$\begin{equation} \frac{d E_\mathrm{GW}}{df} = \frac{3}{5}\frac{G}{c^5} (2\pi f)^2 |\tilde{A}|^2, \end{equation} \tag{ 15 }$$

where $\tilde{A}$ denotes the Fourier transform of $A \equiv \frac{d^2}{dt^2}\, I\hspace{-5pt}{\scriptscriptstyle \stackrel{ \hbox{--}}{}}_{zz}$ computed via

$$\begin{equation} \tilde{A}(f) = \int _{-\infty }^{\infty } A(t)\, e^{-2\pi i f t} dt\,\,. \end{equation} \tag{ 16 }$$

At this point, it is useful to define for future reference the dimensionless characteristic GW strain (Flanagan & Hughes 1998), in terms of the GW spectral energy density,

$$\begin{equation} h_\mathrm{char} = \sqrt{\frac{2}{\pi ^2} \frac{G}{c^3} \frac{1}{D^2} \frac{{\it dE}_\mathrm{GW}}{df}}\,\,. \end{equation} \tag{ 17 }$$

For signals with relatively stable frequencies and amplitudes, Fourier transforms and their energy spectra are adequate frequency analysis tools. However, for signals with time-varying amplitudes and frequencies, a short-time Fourier transform (STFT) is more appropriate. The STFT of A(t) is

$$\begin{equation} \tilde{S}(f,\tau) = \int _{-\infty }^{\infty } A(t)\, H(t - \tau) \, e^{-2\pi i f t} dt, \end{equation} \tag{ 18 }$$

where τ is the time offset of the window function, H(t − τ). We use the Hann window function:

$$\begin{equation} H(t-\tau) = \left\lbrace \begin{array}{@{}lcl} \frac{1}{2} \left(1 + \cos \left(\frac{\pi (t - \tau)}{\delta t}\right) \right) & \mbox{for} & |t-\tau | \le \displaystyle\frac{\delta t}{2}\\ 0 & \mbox{for} & |t-\tau | > \displaystyle\frac{\delta t}{2}\\ \end{array} \right., \end{equation} \tag{ 19 }$$

where δt is the width of the window function. The analog of the energy spectrum of the Fourier transform is the spectrogram, $|\tilde{S}(f,\tau)|^2$. Using the spectrogram, we define an analog to the energy emission per frequency interval (Equation (15)):

$$\begin{equation} \frac{d E^*_\mathrm{GW}}{df}(f,\tau) = \frac{3}{5}\frac{G}{c^5} (2\pi f)^2 |\tilde{S}(f,\tau)|^2 \,. \end{equation} \tag{ 20 }$$

We emphasize that the GW strains reported in this paper are based upon matter motions alone and do not include the low-frequency signal that results from asymmetric neutrino emission (Burrows & Hayes 1996; Müller & Janka 1997). Accurate calculations of asymmetric neutrino emission require multi-dimensional, multi-angle neutrino transport to capture the true asymmetry of the neutrino radiation field (see, e.g., Ott et al. 2008). Our choice to parameterize the effects of neutrino transport by local heating and cooling algorithms is based upon assumptions of transparency, which ignore diffusive effects and would exaggerate the asymmetries and resulting GWs. For example, Kotake et al. (2007) estimated the neutrino GW signal using a similar heating and cooling parameterization and obtained GW strain amplitudes that are ∼100 times the matter GW signal. However, with an improved ray-tracing-based method, the same authors find much smaller amplitudes that are larger than those due to matter motions by only a factor of a few (Kotake et al. 2009). This is in agreement with the GW estimates of Marek et al. (2009) who used 1D ray-by-ray neutrino transport and coupled neighboring rays in 2D hydrodynamic simulations.

Studying the matter GW signal alone is worthwhile. Although the neutrino GW strain amplitudes can be as large or even larger than the contribution by matter (Burrows & Hayes 1996; Müller & Janka 1997; Müller et al. 2004; Marek et al. 2009), the typical frequencies, f, of the neutrino GW signal (∼10 Hz or less) are typically much lower than the frequencies of the matter signal (≳100 Hz). Consequently, the GW power emitted, which is proportional to f², can be much higher for the matter GW signal. Furthermore, although future GW detectors (e.g., Advanced LIGO) will have improved sensitivity at low frequencies, current detectors have response curves that are not sensitive to the lower frequencies of the neutrino GW signal.

In Figure 1, we plot the GW strain (Equation (13)) times the distance to a 10 kpc source, h₊D, versus time after bounce for all simulations. Though there is some diversity in amplitude and timescale among these GW strains, there are several recurring features that exhibit systematic trends with mass and neutrino luminosity. We illustrate these features in Figure 2 with the GW strain of the simulation using the 15 M_☉ progenitor and $L_{\nu _e} = 3.7 \times 10^{52}$ erg s⁻¹. Before bounce, spherical collapse results in zero GW strain. Just after bounce the prompt shock loses energy and stalls, leaving a negative entropy gradient that is unstable to convection. Because the speeds of this prompt convection are larger than those of steady-state postshock or PNS convection afterward, the GW strain amplitude rises to h₊D ∼ 5 cm during prompt convection and settles down to ∼1 cm roughly 50 ms later, which is consistent with the results of Ott (2009b) and Marek et al. (2009). Later in this section, we show that during both phases, convective motions in postshock convection above the neutrinosphere and PNS convection below it contribute to the GW strain. Since nonlinear SASI oscillation amplitudes increase around 550 ms past bounce, the GW signal strengthens from h₊D ∼ 1 to 10 cm and is punctuated by spikes that are coincident in time with narrow plumes striking the PNS "surface" (at ∼50 km). Marek et al. (2009) also noted this correlation.

Zoom In Zoom Out Reset image size

The final feature after ∼800 ms is associated with explosion. The signatures of explosion are twofold. First, during explosion, postshock convection and the SASI subside in strength and the higher frequency (∼300–400 Hz) oscillations in h₊D diminish. Second, global asymmetries in mass ejection result in long-term and large deviations of the GW strain. In Figure 2, a monotonic rise of h₊D to nonzero, specifically positive, values corresponds to a prolate explosion in this simulation. This is similar to the "memory" in the GW signal of asymmetric neutrino emission (Burrows & Hayes 1996; Müller & Janka 1997). When the explosion is spherical, the strain drops to zero and remains there, and when it is oblate, the strain maintains negative values. Examples of h₊D curves showing prolate, oblate, and spherical explosions are shown in Figure 3. The simulation using a 20 M_☉ progenitor and $L_{\nu _e} = 3.4 \times 10^{52}$ erg s⁻¹ (gray line) exploded with an oblate structure, the simulation with M = 12 M_☉ and $L_{\nu _e} = 3.2 \times 10^{52}$ erg s⁻¹ (light gray line) has a prolate explosion, and the simulation with M = 12 M_☉ and $L_{\nu _e} = 2.2 \times 10^{52}$ erg s⁻¹ explodes almost spherically.

Zoom In Zoom Out Reset image size

Figure 4 shows snapshots of the entropy (in units of k_B/baryon) distribution during explosion for the three models highlighted in Figure 3. Lighter shades (warmer colors in the online version) represent higher entropies, while darker shades (cooler colors in the online version) represent lower entropies. The general shapes of the matter interior to the shocks is oblate (M = 20 M_☉ and $L_{\nu _e} = 3.4 \times 10^{52}$ erg s⁻¹), prolate (M = 12 M_☉ and $L_{\nu _e} = 3.2 \times 10^{52}$ erg s⁻¹), and spherical (M = 12 M_☉, $L_{\nu _e} = 2.2 \times 10^{52}$ erg s⁻¹; and thin black line). We emphasize that the GW signal is sensitive to ℓ = 2 accelerations of matter and is somewhat blind to differences in composition or higher order asymmetries. Therefore, even though the GW signal may indicate that the explosion is in general "spherical," "oblate," or "prolate," the entropy (hence temperature and composition) distributions may not be.

Zoom In Zoom Out Reset image size

In Figures 5 and 6, we localize the source of GW emission. Figure 5 shows the spatial distribution of

$$\begin{equation} \left| \frac{d}{dt} \left[ r^3 \rho \left(P_2(\cos \theta)\, v_r + \frac{1}{2} \frac{\partial }{\partial \theta } P_2(\cos \theta)\, v_\theta \right) \right] \right|, \end{equation} \tag{ 21 }$$

which is the time derivative of the integrand of Equation (12) and determines the GW strain (Equation (13)). As in Figure 2, these results are for the 15 M_☉ progenitor and $L_{\nu _e} = 3.7 \times 10^{52}$ erg s⁻¹. Lighter shades (brighter colors in the online version) correspond to a larger contribution to the GW strain integral. The left panel is at 615 ms past bounce, and represents the phase with a stalled shock (∼250 km) and vigorous postshock-convection/SASI motions. This is our first indication that motions below the gain radius (∼100 km), in particular deceleration of plumes and PNS convection, are the strongest sources of GWs before explosion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The right panel is at 1060 ms past bounce and shows that accelerations at an aspherical shock (∼4000 km at this time) are responsible for the "memory" signature at late times. Using the appropriate shock velocities and asymmetries and postshock densities and velocities in Equations (12) and (13), we obtain the correct order of magnitude for the memory signal. However, during explosion a high entropy wind emerges that encounters the swept-up material of the primary shock and produces secondary shocks. We find that these secondary shocks and their asymmetric structure can add a non-trivial contribution to the "memory" signal. The strength of these winds is a strong function of the neutrino luminosity (Qian & Woosley 1996; Thompson et al. 2001). Since we employ a constant neutrino luminosity, an accurate characterization of the late-time "memory" (such as saturation) are beyond the scope of this paper.

Figure 6 quantifies the sources of GW radiation and their spatial distribution as a function of radius. The model shown is the same as in Figures 2 and 5. For reference, the entire GW signal is shown in both panels (solid-black line). The signals originating from >50 km (orange, online version) and <50 km (blue, online version) are shown in the top panel. Below 50 km, PNS convection dominates the motions, and above this radius postshock convection and the SASI dominate. As expected, most of the signal associated with prompt convection originates in the outer convective zone, though motions below 50 km and presumably from PNS convection account for a fair fraction (∼20%). Afterward, the contributions to the GW amplitude from below and above 50 km are comparable. This suggests that motions associated with both PNS convection and the postshock-convection/SASI region contribute significantly to the GW signal. The "memory" signature of the GW strain during explosion clearly originates from the outer (exploding) regions. Interestingly, once the model explodes, and the nonlinear postshock-convection/SASI motions subside, the GW signal from below 50 km diminishes as well, though PNS convective motions do not.

To further understand the origin of GW emission, the bottom panel of Figure 6 refines the spatial origin of the GW emission into five regions. We plot the GW strain from five overlapping regions, and each extends from the center to different outer radii (30, 40, 50, 60, and 100 km). PNS convection extends from ∼20 to ∼40 km and the gain radius is ∼100 km. However, as we explain in section Section 3.4, the turbulent motions of postshock convection/SASI penetrate below the gain radius to radii of ∼60 km during the most vigorous phases. Therefore, the partial GW signals with outer radii of 30 and 40 km contain signals only due to PNS convection; 50 km encompasses PNS convection and gravity waves that are excited by the overlying convection/SASI; 60 km encompasses PNS convection, gravity waves, and the deepest penetration of the convection/SASI plumes (see Section 3.4); and 100 km encompasses all of these contributions and the gain radius. First, we note that extending the outer radius from 60 to 100 km, the gain radius, adds very little signal during the most vigorous postshock-convection/SASI phase from 550 to 800 ms past bounce. Therefore, the most relevant motions are those at ∼60 km and below, which include the postshock-convection/SASI plume decelerations, excited gravity waves, and PNS convection. Furthermore, the strength of the GW signal from all radii increases and decreases in synchrony with the vigor of postshock-convection/SASI motions, which are restricted to radii above ∼50 km. This suggests that the influence of the convection/SASI diminishes only gradually with depth. Even though the mechanism for convection/SASI and their nonlinear motions are above 50 km, they influence PNS convective motions below 50 km, which results in GW emission whose vigor is ultimately due to the postshock-convection/SASI motions.

Hence, we conclude that features of the SASI and postshock convection, in particular the plumes striking the PNS "surface," are ultimately responsible for the strong GW signal from ∼550 ms until explosion at ∼800 ms after bounce. Some of the signal comes from the accelerations at this interface. The rest comes from gravity waves and motions within the PNS (<50 km) that are excited by the plumes striking this layer. Whether these motions are excited g-waves, enhanced PNS convection, or a combination of both is difficult to distinguish in our simulations. A focused investigation on this aspect is warranted.

Müller et al. (2004) conclude that the GW signal from PNS convection is a few factors smaller in amplitude than the GWs from postshock convection. At first sight, this conclusion appears at odds with our results that the signals above and below 50 km are roughly equal. However, to calculate the GW signal from PNS convection, Müller et al. (2004) used results from a simulation that extended out to only 60 km (Keil et al. 1996), which ignores any influence that postshock-convective or SASI motions have on PNS convection. This is consistent with the GW signal we obtain after explosion, when postshock convection and SASI go away, but PNS convection remains. These results emphasize the importance of calculating PNS convection in the full context of core-collapse dynamics.

We now return to Figure 1 and analyze the trends with progenitor mass and $L_{\nu _e}$. Each panel represents a single progenitor model and the curves are labeled by $L_{\nu _e}$. The accretion rates (in M_☉ s⁻¹) at the stalled shock of non-exploding models are shown for comparison. The prominent features described in Figure 2 are apparent in all models, except the run using M = 40 M_☉ and $L_{\nu _e} = 6.0 \times 10^{52}$ erg s⁻¹ that does not explode and manifests the larger amplitudes associated with the nonlinear SASI at much later times. Otherwise, lower luminosity simulations take longer to develop the nonlinear SASI GW signal and explode at later times, but all runs seem to saturate at roughly similar amplitudes. This is consistent with the results of Murphy & Burrows (2008b), in which they note that lower luminosity runs take longer to reach saturated shock oscillations, but that the saturated amplitudes are similar for the different neutrino luminosities. Similarly, Kotake et al. (2009) conclude that the GW amplitudes are independent of neutrino luminosity. Finally, we note that the significant drop in accretion rate in simulations using the 20 M_☉ progenitor, instigates strong nonlinear SASI motions and associated GW signatures. This might suggest using GWs as a diagnostic for the location of shell burning or compositional interfaces. However, very few of the other simulations show such a convincing correlation.

The neutrino luminosities required to explode the 40 M_☉ progenitor model are much larger than simulations obtained that consistently calculate neutrino transport (Dessart et al. 2006; Ott et al. 2008; Marek & Janka 2009). Therefore, if other mechanisms fail to explode higher mass stars, we suspect that the simulation using $L_{\nu _e} = 6 \times 10^{52}$ erg s⁻¹, which does not explode by the neutrino mechanism, gives an upper limit on the most likely GW emission strength for massive progenitors that collapse to form black holes without explosion. Taken at face value, the low convective and SASI motions of this model imply that, in the absence of rapid rotation, the GW power emitted by the most massive progenitors that form black holes could be quite low.

If rotation rates are low and PNS core g-modes are weak, then postshock-convection/SASI motions are the most probable source of GWs. In this scenario, the GW emission of Galactic CCSNe (∼10 kpc) is probably too weak to detect the GW strain waveform directly as a time series. Instead, GW observers will search for excess spectral power above the detector noise in time–frequency maps (e.g., Abbott et al. 2009). The energy spectra versus frequency for all simulations is plotted in Figure 7 assuming a distance to the source, D, of 10 kpc. Rather than showing dE_GW/df, we plot h_char (Equation (17)), which is proportional to $\sqrt{{\it dE}_{\mathrm{GW}}/df}$. As in Figure 1, each panel shows the spectra from simulations for the same progenitor model, and within each panel the spectra are color-coded and labeled by $L_{\nu _e}$. For reference, the noise thresholds for Initial LIGO (solid black line, Gustafson et al. 1999), Enhanced LIGO (dot-dashed-black line; R. Adhikari 2009, private communication), and Advanced LIGO (burst mode, dashed-black line; D. Shoemaker 2006, private communication) are included.

Zoom In Zoom Out Reset image size

Characteristically, the spectra show broad peaks with maximum power at ∼300 to ∼400 Hz. The lower and higher frequencies correspond to lower and higher progenitors masses, respectively. Later in this section and Section 3.4, in analyzing the spectra as a function of postbounce time, we show in Figure 8 that simulations with lower $L_{\nu _e}$ explode later and obtain higher frequencies at later times. However, the spectra in Figure 7, which are calculated using the entire time sequence, do not clearly show such a correlation. The spectra associated with 12, 15, and 20 M_☉ progenitors show a second, weaker peak at ∼100 Hz that we associate with prompt convection. Generally, all spectra are only marginally above or just below the design noise threshold for Initial LIGO. On the other hand, all simulations have power above the design noise threshold for Enhanced and Advanced LIGOs, with maximum spectral signal-to-noise ratios of ∼5 and ∼10–20, respectively.

Zoom In Zoom Out Reset image size

Characteristic GW amplitudes, frequencies, optimal theoretical single-detector signal-to-noise ratios, and energies for all simulations presented in this paper are listed in Table 1. The first two columns identify the progenitor mass and neutrino luminosity of the simulation. Following are the maximum amplitude of the GW strain (|h_+,rms|), the signal-to-noise ratios with respect to the Initial, Enhanced, and Advanced LIGO sensitivity curves (S/N_LIGO, S/N_eLIGO, and S/N_advLIGO). We also provide values for the emitted GW energy (E_GW) and h_char (Equation (17)) for the maxima of the characteristic strain spectra and the frequencies f_char,max at which they are located. The values in this table do not present clean, monotonic trends with progenitor mass or neutrino luminosity. However, we note some general trends. More massive progenitors tend to produce higher frequencies.¹⁰ On the other hand, for a given progenitor model, higher neutrino luminosities give lower frequencies. In Section 3.4, we explain that the characteristic frequency of GWs is connected with the buoyancy frequency in the postshock region, which is higher for more massive progenitors, increases with time, and is largely independent of neutrino luminosity. This explains the first correlation with progenitor mass, but to explain the anti-correlation with neutrino luminosity, we note that the higher luminosity simulations explode earlier when the buoyancy frequency is lower.

Table 1. Quantitative Summary of GW Characteristics^a

Massb	$L_{\nu _e}$c	|h+,max|d	S/NLIGOe	S/NeLIGOf	S/NadvLIGOg	hchar,maxh	fchar,maxi	EGWj
12	1.8	0.27	0.8	2.0	7.9	1.04	253	6.1
12	2.2	0.20	0.7	1.7	6.7	0.61	200	3.5
12	2.8	0.46	1.2	2.7	10.8	1.35	273	8.5
12	3.2	1.13	1.6	4.0	15.7	2.23	294	20.4
15	3.2	0.37	1.1	2.6	11.4	2.14	345	17.5
15	3.4	0.40	1.0	2.4	10.2	1.49	406	14.4
15	3.7	0.37	1.1	2.8	11.5	1.83	365	12.8
15	4.0	1.53	2.1	5.3	21.2	3.10	347	46.1
20	3.2	0.32	1.4	3.4	14.5	3.00	348	29.5
20	3.4	0.70	1.8	4.5	19.5	4.68	347	57.8
20	3.6	0.56	1.5	3.9	16.1	2.67	429	34.5
20	3.8	0.48	1.5	3.8	15.7	3.14	369	33.8
40	6.0	0.33	1.0	2.6	12.2	3.23	423	36.2
40	10.0	0.68	1.4	3.6	17.3	3.93	359	47.1
40	12.0	0.82	1.4	3.4	14.0	2.04	323	16.6
40	13.0	4.53	0.8	1.8	9.4	1.08	1k	4.1

Notes. ^aThis table lists the integrated GW characteristics of the 2D simulations. These simulations represent a two-dimensional parameterization that investigates the dependence of GW emission on progenitor mass (Column 1) and neutrino luminosity (Column 2). ^bProgenitor model (M_☉). ^cNeutrino luminosity (10⁵² erg s⁻¹). ^dMaximum GW strain (10⁻²¹ at 10 kpc). ^eOptimal theoretical single-detector signal-to-noise ratio using the Initial LIGO sensitivity curve (Gustafson et al. 1999). ^fOptimal theoretical signal-to-noise ratio using the Enhanced LIGO sensitivity curve (R. Adhikari 2009, private communication). ^gOptimal theoretical single-detector signal-to-noise ratio using the burst-mode Advanced LIGO sensitivity curve (D. Shoemaker 2006, private communication). ^hMaximum of the characteristic strain spectrum defined in Equation (17) (10⁻²¹ at 10 kpc). ⁱFrequency location of h_char,max (Hz). ^jGW energy emitted (10⁻¹¹ M_☉ c²). ^kBecause this simulation explodes early and with large asymmetry, the low-frequency "memory" signature in the GW strain dominates the energy spectrum.

Download table as:  ASCIITypeset image

In general, the typical GW strain for a source at 10 kpc is a few times 10⁻²² to 10⁻²¹, and typical values of E_GW range from a few times 10⁻¹¹ to a few times 10⁻¹⁰ M_☉ c². Kotake et al. (2009) conclude that the stochastic properties of the SASI preclude any relationship between the neutrino luminosity and the GW amplitude. While there appears to be no trends, we note several competing dependencies on neutrino luminosity that result in a non-monotonic and complicated relationship. In Section 3.4, we connect the GW amplitude to the speed of plumes, v_p, at the base of the SASI region. In addition, we note that although lower luminosity simulations take longer to saturate the SASI motions, the saturation value eventually achieved is similar for all runs. Since most of the GW power is emitted during the phase of strongest postshock-convection/SASI motions, the total energy emitted is determined by the duration of this phase, which starts earlier for the higher neutrino luminosity, but also ends much earlier at the higher luminosities, so that the total duration of the stronger SASI phase, and therefore, E_GW is a non-monotonic function of neutrino luminosity.

Focusing on the energy spectrum for an entire GW strain data set is most appropriate when the signal is regular in both amplitude and frequency. However, Figures 1 and 2 show that the amplitude and frequency change substantially in these simulations. In this case, a time–frequency analysis, or spectrogram, is more appropriate (see Section 3.1). To obtain energy spectra versus time we take the Fourier transform of the GW strain convolved with a Hann window function with width ∼20 ms and sample the time domain, τ in Equation (18), at intervals of 2 ms. In Figure 8, we show color maps of dE*_GW/df (Equation (20)) versus frequency and time. The warmer colors in Figure 8 reflect higher power, and vice versa for cooler colors. In general, simulations that explode later (i.e., with higher mass progenitors and/or lower neutrino luminosities) have GW power at higher frequencies. During explosion, the frequency at peak power drops to 10s of Hz, owing to asymmetric mass ejection.

Unlike in the spectra, which are calculated using the entire time domain (Figure 7), prompt convection, nonlinear SASI, SASI plumes, and explosion are quite apparent in these time–frequency plots. Figure 9 labels these features in the spectrogram for the simulation using M = 15 M_☉ and $L_{\nu _e} = 3.7 \times 10^{52}$ erg s⁻¹. The three features showing the largest power are prompt convection just after bounce, nonlinear SASI motions, and explosion. Though the power declines between prompt convection and the start of the nonlinear SASI, the frequency at peak power increases monotonically from ∼100 Hz at bounce to ∼300–400 Hz just before explosion. In the next section, we propose that the characteristic timescale for the deceleration of plumes striking the PNS "surface" determines the GW peak frequencies. The solid black line in Figure 9 shows our analytic estimate for this characteristic frequency, f_p, and that it agrees with the frequency and time evolution of the GW signal from simulation.

Zoom In Zoom Out Reset image size

Figure 6 shows that motions associated with the broad region, which includes PNS convection, postshock convection, and the SASI, contribute to the overall GW signal. Marek et al. (2009) come to a similar conclusion. However, in Sections 3.2 and 3.3, we noted that the strength of the GW signal waxes and wanes as the vigor of postshock-convection and SASI motions rises and falls. Moreover, spikes in the GW strain often correlate in time with downdrafts (plumes) striking the PNS "surface" (also noted by Marek et al. 2009). Hence, we argue that the strength and characteristic frequencies of the GW emission are determined predominantly by the deceleration of the plumes striking the PNS "surface."

In other contexts, the "surface" of the PNS is defined as the position of the neutrinosphere, the location of a steepening density gradient, or the narrow region over which the average radial velocity of the accreting matter approaches zero. These definitions roughly coincide in space and are in fact intimately related, but they neither explain nor quantify the strong decelerations of the downdrafts in an otherwise continuous medium. We note, however, that buoyancy forces and their radial profile at the boundary between convective and stably stratified regions explain the deceleration of the plumes and, hence, the GW amplitudes and frequencies. We, therefore, define the PNS "surface" as the region where the negatively buoyant downdrafts become positively buoyant and reverse their downward motions.

An important quantity in analyzing convection and buoyancy is the square of the Brunt–Väisälä (or buoyancy) frequency

$$\begin{equation} N^2 = \left(\frac{G M_r}{r^3}\right) \left(\frac{1}{\Gamma _1} \frac{d \ln P}{d \ln r} - \frac{d\ln \rho }{d \ln r}\right), \end{equation} \tag{ 22 }$$

where M_r is the mass interior to the radius, r, and the thermodynamic derivative, Γ₁ = (∂ln P/∂ln ρ)_S, is evaluated at constant entropy, S. In the local dispersion relation for waves, this corresponds to the characteristic frequency associated with gravity waves (not to be confused with GWs), whose dominant restoring force is buoyancy, which sets the frequency scale. Alternatively, N² < 0 indicates linear instability and is in fact the Ledoux (and Rayleigh–Taylor) condition for convection. Naively, one might expect the radius where N² changes sign from negative to positive to mark the lower boundary of the convection zone. However, the Ledoux condition is derived from a linear analysis and ignores important nonlinear effects. In practice, convective plumes have momentum when they reach this boundary and penetrate beyond the Ledoux condition boundaries. This is called overshoot and has a long history and many proposed prescriptions, though many of these are either inappropriate for the dynamics in core-collapse SNe or lack adequate physical motivation, resulting in a free parameter (see Meakin & Arnett 2007b; Arnett et al. 2009; and the references therein). For this paper, we adopt an analysis based upon buoyancy and the bulk Richardson number (Meakin & Arnett 2007b), which is physically motivated and contains no free parameters.

Beyond the boundaries defined by the Ledoux condition, N² is positive and, therefore, coherent plumes that penetrate this boundary experience a buoyancy force back toward the central regions of convection. In this context, the buoyancy acceleration, b(r), felt by a Lagrangian mass element displaced from r_i to r is related to the buoyancy frequency by

$$\begin{equation} b(r) = \int ^{r_i}_{r} N^2(r) dr \,. \end{equation} \tag{ 23 }$$

For plumes approaching the lower boundary of postshock convection, the magnitude of the integrand is largest where N²>0. Therefore, we take as r_i the lower Ledoux condition boundary. With a plume velocity in the radial direction, v_p, at r_i as initial conditions and Equation (23) as the position dependent acceleration, we integrate the plume's equation of motion to determine the depth of penetration, D_p. Because the sources of GWs are time-changing quadrupole motions, the inverse of the characteristic time of the buoyancy impulse for each plume should correspond roughly to the peak frequencies in the GW spectra and spectrograms (Figures 7 and 8). This timescale, t_p, is defined as the half-width at half-maximum (HWHM) of the buoyancy acceleration pulse and the associated frequency is f_p = 1/(2πt_p).¹¹ It is this definition of f_p for which we give quantitative results in this paper.

Using the bulk Richardson number, a dimensionless measure of the boundary layer stiffness, we obtain an analytic description of f_p that is independent of, and corroborates, our method for calculating f_p. There is a long tradition of using this dimensionless parameter in atmospheric sciences to accurately describe the boundary layer in atmospheric convection (Fernando 1991; Fernando & Hunt 1997; Bretherton et al. 1999; Stevens 2002). More recently Meakin & Arnett (2007b) have successfully used similar arguments to describe the boundaries of convection in stellar evolution. The bulk Richardson number is

$$\begin{equation} R_b = \frac{\Delta b D_p}{v_p^2}, \end{equation} \tag{ 24 }$$

where Δb is the change in buoyancy acceleration across the boundary, D_p is a length scale such as the penetration depth, and v²_p is the typical velocity of the plume. This is not to be confused with the gradient Richardson number, R_g, which is derived from a linear analysis and is a condition for shear instabilities in a stratified medium. While R_g = 1/4 is the derived value that separates stability from instability in a linear analysis, R_b is a rough measure of the stiffness of a stable layer next to a convective layer in nonlinear flows. Approximately, it is the ratio of work per mass done by buoyancy (ΔbD_p) and the kinetic energy per mass of the plume (v²_p). Therefore, R_b ≲ 1 offers little resistance to penetration, R_b>1 represents a stiff boundary, and where R_b ∼ 1, the work done by buoyancy balances the plume's kinetic energy, causing it to turn around. Hence, R_b is appropriate in characterizing the penetration of downdrafts into the underlying stable layers.

To estimate scales and trends, we approximate the integral in Equation (23), giving

$$\begin{equation} \Delta b \sim N^2 D_p, \end{equation} \tag{ 25 }$$

and we assume that R_b ∼ 1, which we argued above is where buoyancy balances the plume's kinetic energy and gives the turn-around point. This gives (ND_p/v_p)² ∼ 1, or

$$\begin{equation} D_p \sim \frac{v_p}{N}. \end{equation} \tag{ 26 }$$

We note that this scaling is similar to the gravity wave displacement amplitude derived by Meakin & Arnett (2007a), which they derive by equating the gravity wave pressure fluctuations to the ram pressure of the convective eddies at the boundary. An estimate of the characteristic timescale is given by the velocity divided by the acceleration, t_p ∼ v_p/Δb. Substituting in Equations (25) and (26) gives¹¹

$$\begin{equation} f_p \sim \frac{N}{2 \pi }. \end{equation} \tag{ 27 }$$

Interestingly, this characteristic frequency is insensitive to the plume velocity and penetration depth. Rather, f_p is most sensitive to the buoyancy frequency at the turning point (N_turn).

Although the above analysis suggests that the characteristic frequency is independent of v_p, the amplitude of the GW strain is directly dependent upon these velocities. Very roughly, the GW strain is

$$\begin{equation} h_+ \sim \frac{4 \pi G}{Dc^4} \rho \, f_p \, v_p \, r^3 \Delta r \, \Delta \mu, \end{equation} \tag{ 28 }$$

where we use Equations (12) and (13), Δμ is the finite difference approximation of dcos θ, and we assume that ∂(ρv_p)/∂t ∼ ρ∂v_p/∂t and that ∂v_p/∂t ∼ f_pv_p. Figure 10 compares h₊D with the maximum downward plume speed below 120 km, v_p, as a function of time after bounce. The spikes in v_p correspond to the strongest plumes that strike the PNS "surface," and the baseline, below which v_p does not dip, shows the time-averaged accretion of material onto the PNS. As the PNS accumulates more mass, exerting stronger gravitational forces, the background and plume velocities evolve upward slightly until ∼550 ms after bounce. At this time, the motions in the postshock-convection/SASI region increase dramatically, leading to a significant rise in v_p. After the start of explosion (≳800 ms), it no longer makes sense to define a plume velocity, since the GW signal during explosion is dominated by explosion dynamics. Hence, we do not show v_p after ∼800 ms.

Zoom In Zoom Out Reset image size

Initially, prompt convection dominates the signal and v_p is not relevant. As prompt convection settles to steady-state convection, the plumes become more relevant. From ∼150 ms to ∼550 ms, both v_p and the GW amplitude are low, but as v_p grows so does the GW amplitude. Strikingly, the increase in GW wave strain at ∼600 ms coincides with the rise in v_p. In fact, many, but not all, spikes in h₊D coincide with rapid changes in v_p. Furthermore after ∼600 ms, substituting values for typical flows near the boundary (e.g., ρ ∼ 10¹¹ g cm⁻³, v_p ∼ 4 × 10⁹ cm s⁻¹, r ∼ 60 km, Δr ∼ D_p ∼ 20 km, f_p ∼ 300 Hz, and Δμ ∼ D_p/(πr)) into Equation (28) gives h₊D ∼ 5 cm, which is the order of magnitude of the GW strain.

Figures 11 and 12 depict the buoyancy frequency, N (solid blue line, online version), (GM_r/r³)^1/2 (dashed purple line, online version), and the radial velocity, v_r (black points), as a function of radius at 300 and 620 ms after bounce for the 2D simulation that uses the 15 M_☉ progenitor model and $L_{\nu _e} = 3.7 \times 10^{52}$ erg s⁻¹. The velocities clearly show (from large to small radii), spherical infall, the asymmetric shock, the postshock-convection/SASI region, a stable layer, and PNS convection. At 300 ms (Figure 11), the dynamics are dominated by postshock-convection and mild SASI oscillations, while at 620 ms, the SASI oscillations have increased substantially. Roughly, the buoyancy frequency scales with the average density, Equation (22). However, at the core, the spatial derivatives, and therefore N, approach zero. In the regions of PNS convection (∼20 to ∼40 km) and postshock convection (≳100 km), N² is negative. Green line segments in both figures indicate the regions that are unstable to convection by the Ledoux condition in 1D simulations, and the penetration depths, D_p, which are computed by integrating the plume's equation of motion, are labeled and indicated by an orange line. At 300 and 620 ms after bounce, D_p at the base of the region of postshock convection is ∼25 and ∼39 km, respectively. These analytic penetration depths agree with those obtained from the numerical simulations. Note that the boundaries of PNS convection in the 2D simulations coincide with the boundaries given by the Ledoux condition. This is due to the fact that the relatively large buoyancy frequency and small convective speeds imply very small overshoot depths (Equation (26)) at the edges of PNS convection. The plume frequencies, f_p, obtained from the HWHM timescale of the acceleration pulses, are ∼196 Hz at 300 ms and ∼272 Hz at 620 ms.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To better understand the contributions to f_p, we plot in Figure 13 the characteristic frequency (f_p), the buoyancy frequency at D_p (N_turn/(2π)), and (GM_r/r³)^1/2/(2π) at D_p as a function of time after bounce. Strikingly, f_p and N_turn/(2π) agree quantitatively. Since f_p is calculated using the equation of motion and is independent of the approximations and assumptions that led to Equation (27), the fact that f_p is nearly equal to N_turn/(2π) for all times provides implicit confirmation of our approximations and assumptions in deriving the analytic result. Furthermore, we calculated the ratio N_turnD_p/v_p (the square root of R_b after inserting the approximation that Δb ∼ N²D_p) for all models (Figure 14), found that its value is around 1.5–2, and that it changes less than 30% during the simulation. This validates both the assumption that R_b ∼ 1 characterizes where the plumes turn around and the approximation that Δb ∼ N²D_p.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Initially, all three measures of frequency in Figure 13 coincide. However, while N_turn/(2π) tracks the time evolution of f_p, (GM_r/r³)^1/2/(2π) diverges from the other two measures. By plotting both N_turn and (GM_r/r³)^1/2, which is the inverse free-fall time, we highlight global (inverse free-fall time) and local (difference of logarithmic derivatives) properties of the PNS that contribute to Equation (22). The inverse free-fall time is proportional to the square root of the average density (compactness), which is strongly dependent upon the dense-matter EOS. Marek et al. (2009) reported a change in peak GW frequency with stiff and soft dense-matter equations of state, though they make reference to only compact PNSs without explaining the timescales in the context of the buoyancy frequency. Figure 13 shows that both the global compactness and the local gradients at the penetration depth are important in determining N_turn. These local gradients are strongly dependent upon the details of deleptonization and cooling. Therefore, the time evolution of the characteristic frequencies is a strong function of both the dense-matter EOS in the PNS below the boundary layer and the local microphysics of the boundary layer itself.

Figure 9 compares f_p with the spectrogram, dE*_GW/df (Equation (20)), of the simulation using the 15 M_☉ progenitor and $L_{\nu _e} = 3.7 \times 10^{52}$ erg s⁻¹. From bounce until explosion, the frequency at peak power coincides with f_p. This strong correlation persists whether prompt convection, convection, or the SASI are the dominant instabilities and implies that the plume's buoyant deceleration determines the characteristic GW frequency at all times before explosion.

In Figure 15, we plot f_p as a function of time after bounce for all 2D simulations. Each shade (color in the online version) represents a set of simulations with the same progenitor model but different $L_{\nu _e}$: black for 12 M_☉, blue (online) for 15 M_☉, green (online) for 20 M_☉, and red (online) for 40 M_☉. As in Figure 9, careful inspection of these f_p–time curves show that they agree with the trends of peak frequency in the GW spectrograms, Figure 8. For simulations that use the same progenitor model, their f_p evolutions are practically indistinguishable. The only differences are that the higher luminosity simulations explode earlier (the end of f_p–time curves mark the approximate time of explosion) and consequently have lower frequencies at explosion. In general, the higher mass progenitors produce higher frequencies, with the exception that the 12 and 15 M_☉ models have very similar f_p–time curves.

Zoom In Zoom Out Reset image size

Yoshida et al. (2007) also investigated the interaction between the plumes of postshock-convection/SASI region and the stable layers below it. Their motivation was to investigate the excitation of PNS core g-modes, which are vital for the success of the acoustic mechanism (Burrows et al. 2006, 2007b), by these plumes. They concluded that the driving pressure perturbations had frequencies that are much too low to cause excitation of the PNS core g-modes. As a result, they suggest that core g-mode excitation is inefficient, rather the core g-modes are forced oscillations. However, their analysis relied upon the power spectrum of pressure perturbations at the base of 2D simulations with inner boundaries at ∼50 km. While they attempt to handle the subsonic outflow through this inner boundary consistently, they are forced to set the radial velocity at this inner boundary to maintain the subsonic structure above the boundary. This boundary condition certainly limits the motions at the base of their calculations, and their unfortunate choice of the location of their inner boundary is exactly where the plume's motions are most important for exciting core g-modes.

Furthermore, there is some confusion in the literature concerning the appropriate frequencies and excitation mechanism for the core g-modes. While pressure perturbations play some role in exciting core g-modes, it is well known that buoyancy is the dominant driving force in g-waves (e.g., Unno et al. 1989, Chapter III, Section 13). Therefore, velocity perturbations, which result from both pressure and buoyancy forces, are more relevant in characterizing the perturbation spectrum at the surface of the PNS. The power spectrum of the pressure perturbations is most useful at highlighting the sounds waves that are emitted by oscillating PNS. Furthermore, in their discussion of the relevant frequencies, Yoshida et al. (2007) reference the ∼30 to ∼80 Hz shock-oscillation frequencies in Table 1 of Burrows et al. (2007b) and conclude that motions associated with these frequencies are much too low to excite core g-modes which typically have frequencies around 300 Hz. However, these shock-oscillation frequencies are not relevant in driving oscillations at the surface of the PNS. Rather, we have shown that the characteristic frequencies associated with the plume impulses are most relevant and have frequencies that are very similar to the core g-mode frequencies. Though we address neither the acoustic mechanism nor the core g-modes in this paper, in light of our characteristic frequency analysis, it is possible that the plumes associated with postshock convection/SASI excite core g-mode oscillations (Burrows et al. 2007b).