Equation-of-state Dependence of Gravitational Waves in Core-collapse Supernovae

Each stellar core-collapse simulation we perform uses a progenitor with a zero-age main-sequence mass of 20 M_⊙ from Woosley & Heger (2007) and is started at the onset of collapse in our 2D domain. We utilize the FLASH (v.4) multiscale, multiphysics adaptive mesh refinement framework (Fryxell et al. 2000) modified for core collapse in Couch (2013), Couch & O'Connor (2014), and O'Connor & Couch (2018b). The grid setup is a 2D cylindrical geometry that spans 1 × 10⁹ cm in the radial direction and ±1 × 10⁹ cm along the cylindrical axis. Our maximum block size is 2 × 10⁸ cm, and we allow a total of 10 levels of mesh refinement, giving a minimum block size of ∼3.9 km and minimum grid zone spacing of ∼326 m. Beyond 80 km, we decrease the highest level of refinement attainable by the adaptive mesh refinement routines so as to maintain an angular resolution of at least 082 everywhere. This corresponds to a maximum resolution aspect ratio, Δx_i /r, of about 0.014. We use standard Newtonian hydrodynamics and solve gravity by employing a modified general relativistic effective potential (case A) for the monopole contribution (Marek et al. 2006), retaining spherical harmonic orders up to 16. We note that the use of an effective potential leads to, in general, a quantitative overprediction of the GW frequency due to the inconsistency of the underlying Newtonian hydrodynamics (Müller et al. 2013; Zha et al. 2020).

We incorporate a multidimensional, multispecies, and energy-dependent neutrino transport following an M1 scheme outlined in O'Connor & Couch (2018b). As such, we distinguish three species of neutrino, i.e., electron neutrinos ν_e , electron antineutrinos ${\bar{\nu }}_{e}$, and the composite group of heavy lepton neutrinos/antineutrinos ν_x . For each EOS employed, a set of neutrino opacities is generated using the neutrino transport library NuLib (O'Connor 2015). Neutrino opacities are computed for 12 logarithmically spaced energy groups, the first group centered on 1 MeV and the last one centered on ∼315 MeV. For the emission of electron-type neutrinos and antineutrinos, we include electron and positron capture interactions on protons and neutrons, respectively, as well as electron capture on heavy nuclei following Bruenn (1985). For the emission of heavy lepton neutrino pairs we include emission from electron–positron annihilation to $\nu \bar{\nu }$ pairs following Bruenn (1985) and nucleon–nucleon bremsstrahlung following Burrows et al. (2006) and Hannestad & Raffelt (1998). We include these pair emissions in a parameterized way (Betranhandy & O'Connor 2020). We include isoenergetic scattering of neutrinos of all types on nucleons, nuclei, and alpha particles following Bruenn (1985). For the charged-current interactions on nucleons, as well as the neutral-current scattering interactions, we include corrections for weak magnetism following Horowitz (2002). We also include inelastic scattering of neutrinos on electrons following Bruenn (1985).

For the simulations that include rotation, we map the progenitor model onto the cylindrical grid via an artificial rotation profile,

$$\begin{eqnarray}&&{\rm{\Omega }}(r)={{\rm{\Omega }}}_{0}{\left[1+{\left(\displaystyle \frac{r}{A}\right)}^{2}\right]}^{-1},\end{eqnarray} \tag{ 1 }$$

where $r=\sqrt{{z}^{2}+{R}^{2}}$ is the spherical radius for a given cylindrical radius R and symmetry axis coordinate z, Ω₀ is the central angular speed of the star, and A is the differential rotation parameter. Using this same scheme, Pajkos et al. (2019) find that the parameter A strongly correlates with the stellar core compactness (O'Connor & Ott 2011) and calculate an optimal value of A = 1021 km for the compactness ξ_2.5 = 0.2785 of the 20 M_⊙ progenitor used here.

In this study we focus on nucleonic EOSs where neutrons and protons interact via a Skyrme-type force while immersed in a degenerate gas of noninteracting electrons, positrons, and photons. ⁵ Because leptons and photons are considered noninteracting, for an EOS at a given state their contribution to thermodynamics of the system will be the same regardless of the interactions between nucleons. Furthermore, nucleon contributions dominate the EOS only at densities near or above nuclear saturation density, i.e., number density n_sat ≃ 0.155 baryons fm⁻³ or mass density ρ_sat ≃ 2.7 × 10¹⁴ g cm⁻³, where matter is in a homogeneous state. Moreover, since GWs are mainly emitted from the PNS, which has for the most part densities in excess of one to a few n_sat, we discuss only the bulk nucleon contributions next. For details on the low-density part of the EOS, n ≲ n_sat, where nucleons may cluster to form nuclei, and the lepton and photon contributions, see Lattimer et al. (1985), Lattimer & Swesty (1991), and Schneider et al. (2017).

The FLASH simulation code is set up to solve the hydrodynamic equations of motion where the properties of matter are obtained from the state ζ = (n, y, T) of the system. Here n = n_n + n_p is the baryon number density, where n_n (n_p ) is the neutron (proton) number density, y = n_p /n is the proton fraction, and T is the temperature. We assume that the system is everywhere charge neutral, i.e., the number of protons is equal to the number of electrons minus positrons, and in thermal equilibrium. As an example of the thermodynamic quantity of the EOSs we study, we consider the finite-temperature nucleonic energy per baryon

$$\begin{eqnarray}&&{\epsilon }_{B}(n,y,T)={\epsilon }_{\mathrm{kin}}(n,y,T)+{\epsilon }_{\mathrm{pot}}(n,y).\end{eqnarray} \tag{ 2 }$$

Above, _kin (_pot) is the specific kinetic (potential) energy. In Skyrme-type models, _B (n, y, T) is a function parameterized to fit the properties of bulk nuclear matter and finite nuclei at zero temperature (Lattimer & Swesty 1991; Schneider et al. 2019b). Thus, in this class of models, the potential term is temperature independent and the temperature dependence is fully contained in the kinetic term, to be discussed later. With this in mind, we rewrite Equation (2) in the form

$$\begin{eqnarray}&&{\epsilon }_{B}(n,y,T)={\epsilon }_{\mathrm{thermal}}(n,y,T)+{\epsilon }_{\mathrm{cold}}(n,y),\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{thermal}}(n,y,T)={\epsilon }_{\mathrm{kin}}(n,y,T)-{\epsilon }_{\mathrm{kin}}(n,y,T=0),\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{cold}}(n,y)={\epsilon }_{\mathrm{pot}}(n,y)+{\epsilon }_{\mathrm{kin}}(n,y,T=0).\end{eqnarray} \tag{ 4b }$$

Expression (3) is useful, as _thermal contains only temperature-dependent terms, while _cold is the zero-temperature EOS. Therefore, _cold may be written as a series expansion of the specific energy around nuclear saturation density for isospin symmetric nuclear matter (n_n = n_p ), i.e.,

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{cold}}(n,y)\equiv {\epsilon }_{\mathrm{cold}}(x,\delta )={\epsilon }_{\mathrm{is}}(x)+{\delta }^{2}{\epsilon }_{\mathrm{iv}}(x).\end{eqnarray} \tag{ 5 }$$

As usual in the literature, $x=\tfrac{1}{3{n}_{\mathrm{sat}}}(n-{n}_{\mathrm{sat}})$ and $\delta =1-2y=\tfrac{1}{n}({n}_{n}-{n}_{p})$ is the isospin asymmetry. The expansion in Equation (5) assumes quadratic approximation for the isospin dependence and neglects nonanalytic terms in δ (see Kaiser 2015; Wen & Holt 2021). Up to second order, the isoscalar (is) and isovector (iv) expansion terms are (e.g., Piekarewicz & Centelles 2009; Margueron et al. 2018)

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{is}}(x)={\epsilon }_{\mathrm{sat}}+\displaystyle \frac{1}{2!}{K}_{\mathrm{sat}}{x}^{2}+\cdots ,\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{iv}}(x)={\epsilon }_{\mathrm{sym}}+{L}_{\mathrm{sym}}x+\displaystyle \frac{1}{2!}{K}_{\mathrm{sym}}{x}^{2}+\cdots ,\end{eqnarray} \tag{ 6b }$$

and we recall that, by definition of the nuclear saturation density n_sat, the linear term in x of _is(x) vanishes. The empirical parameters shown in the expansion Equation (6) are the energy per baryon at nuclear saturation density _sat, the isoscalar incompressibility modulus K_sat, the symmetry energy _sym and its slope L_sym, and the isovector incompressibility modulus K_sym. In principle, all of the expansion parameters in Equation (6) can be determined experimentally, with varying degrees of difficulty that usually increase with the order of the parameter. Often, isovector empirical parameters are harder to constrain than their isoscalar counterparts.

We now discuss the kinetic energy per baryon term

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{kin}}(n,y,T)=\displaystyle \frac{1}{n}\left(\displaystyle \frac{{{\hslash }}^{2}{\tau }_{n}}{2{m}_{n}^{\star }}+\displaystyle \frac{{{\hslash }}^{2}{\tau }_{p}}{2{m}_{p}^{\star }}\right).\end{eqnarray} \tag{ 7 }$$

Here τ_t ≡ τ_t (n, y, T) is the kinetic energy density, ${m}_{t}^{\star }\equiv {m}_{t}^{\star }(n,y)$ is the effective mass for a nucleon of type t (t =n, p), and, in the zero-temperature limit, ${\tau }_{t}(n,y,T=0)=\tfrac{3}{5}{\left(3{\pi }^{2}{n}_{t}\right)}^{2/3}{n}_{t}$. We highlight that the quantities τ_t do not depend on any model parameters besides the ones that determine the effective masses ${m}_{t}^{\star }$ (Lattimer & Swesty 1991; Schneider et al. 2017). Hence, the finite-temperature component of the EOS, represented here by _thermal, is fully determined by the effective masses of nucleons (Prakash et al. 1997; Constantinou et al. 2014).

In Skyrme-type models, the inverses of the nucleon effective masses have a simple linear relation with nucleon densities. They are computed from

$$\begin{eqnarray}&&\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{t}^{\star }}=\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{t}}+{\alpha }_{1}{n}_{t}+{\alpha }_{2}{n}_{-t},\end{eqnarray} \tag{ 8 }$$

where α_i are Skyrme parameters and m_t the nucleon vacuum masses. For n_t , if t = n, then −t = p, and vice versa. The parameters α_i are set to reproduce the desired effective mass behavior. Because there are only two such parameters, ${m}_{t}^{\star }$ is fully determined by choosing the value of ${m}_{t}^{\star }$ in two independent points in phase space (n, y) or, alternatively, (x, δ). Note that other forms, some considerably more complicated, for the effective masses are possible and, likely, more realistic (Prakash et al. 1997; Constantinou et al. 2015; Schneider et al. 2019a; Huth et al. 2021). Once the parameters that determine the effective mass of nucleons have been computed, the other parameters of the model, which regulate the potential energy term _pot but were not discussed here for brevity, may be obtained from the nuclear physics constraints in the expansions of Equation (6).

Based on the model discussed above, Schneider et al. (2019b) constructed 97 Skyrme-type EOSs using open-source code SROEOS (Schneider et al. 2017). They built a baseline EOS by fixing 10 EOS observables to their most likely values based on current nuclear physics constraints combined from Danielewicz et al. (2002), Antoniadis et al. (2013), Nättilä et al. (2016), and Margueron et al. (2018); see Table 1 for their values. From the baseline EOS, Schneider et al. (2019b) allowed 1σ and 2σ variations in pairwise combinations of empirically accessible parameters in four sets: s_M = {m^⋆, Δm^⋆}, s_S = {_sym, L_sym}, s_K = {K_sat, K_sym}, and ${s}_{P}=\{{P}_{\mathrm{SNM}}^{(4)},{P}_{\mathrm{PNM}}^{(4)}\}$. Here ${m}^{\star }={m}_{n}^{\star }(n={n}_{\mathrm{sat}},y=1/2)\simeq {m}_{p}^{\star }({n}_{\mathrm{sat}},1/2)$ is the effective mass of neutrons at nuclear saturation density for symmetric nuclear matter (SNM), ⁶ ${\rm{\Delta }}{m}^{\star }={m}_{n}^{\star }({n}_{\mathrm{sat}},0)-{m}_{p}^{\star }({n}_{\mathrm{sat}},0)$ is the neutron–proton mass difference in the limit of pure neutron matter (PNM). ${P}_{{\rm{SNM}}}^{(4)}$ and ${P}_{{\rm{PNM}}}^{(4)}$ is the pressure of SNM and PNM at 4n_sat, respectively. Two parameters were fixed in every EOS: the nuclear saturation density n_sat = 0.155 baryons fm⁻³ and _sat = 15.8 MeV baryon⁻¹, as their uncertainties are small enough to barely affect core-collapse physics. The choice of pressure at a large density as a constraint, which is computed from the relation $P(n,y)={n}^{2}{\left.\tfrac{\partial {\epsilon }_{\mathrm{cold}}(n^{\prime} ,y)}{\partial n^{\prime} }\right|}_{n^{\prime} =n}$, to determine the parameters of the Skyrme model is due to the fact that the pressure of nuclear matter at large densities is significantly better constrained through a combination of high-energy collision analysis and computational models (Danielewicz et al. 2002) and the maximum observed mass of NSs (e.g., Antoniadis et al. 2013) than the higher-order terms in the expansions of Equation (6), which can only be poorly constrained by a mixture of theoretical modeling and astrophysical constraints (Margueron et al. 2018; d'Étivaux et al. 2020).

Table 1. Nuclear Matter Properties of the Baseline EOS, Following Schneider et al. (2020)

Quantity	Baseline	Variation 2σ	Unit
m⋆	0.75	0.20	mn
Ksat	230	30	MeV baryon−1
nsat	0.155		fm−3
sat	−15.8		MeV baryon−1
sym	32		MeV baryon−1
Lsym	45		MeV baryon−1
Ksym	−100		MeV baryon−1
Δm⋆	0.10		mn
${P}_{\mathrm{SNM}}^{(4)}$	125		MeV fm−3
${P}_{\mathrm{PNM}}^{(4)}$	200		MeV fm−3

Note. See Schneider et al. (2019b) for further details. 2σ variations are shown for the two parameters we investigate in this work.

Download table as:  ASCIITypeset image

Finally, a comment on the uncertainties σ is warranted. The variations around the baseline values for m^⋆, Δm^⋆, and the expansion terms in Equation (6) were determined by Schneider et al. (2019b) combining analysis of experimental results and theoretical models by Margueron et al. (2018), astrophysical observations by Antoniadis et al. (2013) and Nättilä et al. (2016), and modeling of heavy ion collisions by Danielewicz et al. (2002). These uncertainties thus have to be taken with caution, as not only are they approximate, but the values of some observables may also be correlated by theory and observations (see, e.g., Tews et al. 2017; Margueron et al. 2018; d'Étivaux et al. 2020).

It is unreasonable, considering current computational resources, to perform multidimensional CCSN simulations for approximately 100 EOSs. Nonetheless, we can still learn a great deal from a smaller subset of simulations. We thus focus on variations in only two empirical parameters and determine their effects on GW signatures. We perform simulations only for the baseline EOS of Schneider et al. (2019b) and four EOSs where either m^⋆ or K_sat differs by 2σ from their baseline values; see Table 1. These choices are motivated from the observation that K_sat (m^⋆) has a strong (weak) effect on the cold NS structure, while it has a weak (strong) effect on the structure of a hot PNS. Since PNSs are hot soon after their birth, we expect the effect of varying m^⋆ while keeping all else fixed to be much more significant in our simulations than varying K_sat while keeping every other empirical parameter, including m^⋆, constant. This is because variations in m^⋆ modify the temperature dependence of the EOS, while modifications in K_sat only affect the cold component. Hence, we expect a similar correlation in our results. Moreover, the relative uncertainty in m^⋆ is approximately twice that of K_sat. Particularly, our baseline simulation has m^⋆ = 0.75 (specified as a fraction of the neutron mass vacuum value) and K_sat = 230 MeV baryon⁻¹. The 2σ value for m^⋆ (K_sat) below baseline is 0.55 (200 MeV baryon⁻¹) and above baseline is 0.95 (260 MeV baryon⁻¹).

Besides studying EOS variations, we perform two additional simulations for both m^⋆ = 0.55 and m^⋆ = 0.95 that include rotation with Ω₀ = 1 rad s⁻¹ and Ω₀ = 2 rad s⁻¹. See Table 2 for an overview of the nine simulations performed in this paper.

All of the EOS tables, and each of the corresponding NuLib tables used in this work, are available at https://doi.org/10.5281/zenodo.4973545.

To extract the GW signal measured by a distant observer, we adopt the standard formula utilizing the second time derivative of the trace-free quadrupole moment. This quadrupole moment in the slow-motion, weak-field formalism is given by (Blanchet et al. 1990; Finn & Evans 1990)

$$\begin{eqnarray}&&{I}_{{ij}}=\int {d}^{3}x({x}_{i}{x}_{j}-\displaystyle \frac{1}{3}{x}^{2}{\delta }_{{ij}})\rho (t,{\boldsymbol{x}}),\end{eqnarray} \tag{ 9 }$$

where the indices range over (x, y, z), δ_ij is the Kronecker delta, and ρ is the source energy density. For axisymmetric sources, the only independent component of the quadrupole moment is I_zz . Unless otherwise stated (for example, in Section 3.4), we compute the first derivative in terms of the fluid velocity, v_i , via the analytic expression (Finn & Evans 1990),

$$\begin{eqnarray}&&\displaystyle \frac{{{dI}}_{{ij}}}{{dt}}=\int {d}^{3}x({x}_{i}{v}_{j}+{v}_{i}{x}_{j}-2{/3\delta }_{{ij}}{x}_{k}{v}_{k})\rho (t,{\boldsymbol{x}}),\end{eqnarray} \tag{ 10 }$$

and the second time derivative is computed numerically. The axisymmetric GW strain for an observer located at the equator a distance D away from the source is then given as (Thorne 1980)

$$\begin{eqnarray}&&{h}_{+}(t,{\boldsymbol{x}})=\displaystyle \frac{3}{2}\displaystyle \frac{G}{{{Dc}}^{4}}\displaystyle \frac{{d}^{2}}{{{dt}}^{2}}{I}_{{zz}},\end{eqnarray} \tag{ 11 }$$

where G is the gravitational constant and c is the speed of light.

The GW spectrograms are extracted by short-time Fourier transforming (STFT) the GW signal,

$$\begin{eqnarray}&&\tilde{S}(f,\tau )={\int }_{-\infty }^{\infty }s(t)H(t-\tau ){e}^{-2\pi {ift}}{dt},\end{eqnarray} \tag{ 12 }$$

into its frequency domain f. τ is the time offset of the Hann window function H(t − τ) (see, e.g., Murphy et al. 2009). The sampling frequency of our simulation output is ∼ 2 × 10⁶ Hz, which is reduced in the post-processing to ∼20,000 Hz. Unless otherwise stated, we use a Hann window width of 40 ms and compute STFTs with a Δτ of 3.6 ms.

To investigate the hypothetical detectability of our simulations, we follow Murphy et al. (2009) to calculate the dimensionless characteristic strain (Flanagan 1998),

$$\begin{eqnarray}&&{h}_{\mathrm{char}}=\sqrt{\displaystyle \frac{2}{{\pi }^{2}}\displaystyle \frac{G}{{c}^{3}}\displaystyle \frac{1}{{D}^{2}}\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{df}}},\end{eqnarray} \tag{ 13 }$$

with the spectral energy density

$$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{GW}}}{{df}}=\displaystyle \frac{3}{5}\displaystyle \frac{G}{{c}^{5}}{\left(2\pi f\right)}^{2}{\tilde{| A| }}^{2},\end{eqnarray} \tag{ 14 }$$

where $\tilde{A}$ denotes the Fourier transform of d² I_zz /dt².

When investigating the effective mass dependence of GWs for nonrotating progenitors, the three simulations we compare are carried out for a minimum of 860 ms post-bounce. Each stellar core collapses until the EOS stiffens and bounce occurs (at 318.5 ± 0.1 ms after the onset of the simulation), followed by the formation and outward propagation of the shock wave that initially stalls at ∼170 km at ∼100 ms after bounce. The top panel of Figure 1 shows the evolution of the mean shock radius. Shock revival occurs for m0.95 at ∼300 ms post-bounce and for m0.75 at ∼900 ms post-bounce, whereas the simulation with the lowest value of the effective mass parameter, m0.55, does not explode in the simulated time. This is consistent with the generally increasing neutrino heating expected for higher effective masses (bottom panel of Figure 1), making conditions more favorable for shock runaway (Janka 2017; Schneider et al. 2019b; Yasin et al. 2020). The difference in neutrino heating between the simulations is a reflection of the neutrino luminosity and neutrino mean energy, both of which increase with the effective mass. We remark that this trend is seen for all neutrino/antineutrino species, with eventual deviations related to the times of explosion.

Zoom In Zoom Out Reset image size

Neutrino heating not only provides the shock with thermal support and aid in shock revival but also drives turbulence and convection in the gain region (Couch & Ott 2015) that imprints on the GWs. Figure 2 shows the GW signals. While we note a strong signal at epochs ≲70 ms that is typically attributed to prompt convection (Murphy et al. 2009), our analysis focuses on epochs after 0.1 s when the quiescent phase has ended and the excitation of PNS oscillations begins. During this phase, such PNS oscillations strongly dominate the GW signal. The amplitude of the GWs increases with higher effective mass. We attribute this, at least in part, to the increasingly (neutrino-driven) "violent" gain region in simulations with higher effective mass, correlating with the amplitude of the PNS oscillations via accretion plumes striking the convectively stable layer below (Murphy et al. 2009; Yakunin et al. 2015). Radice et al. (2019) show, albeit in 3D, that the energy radiated in GWs is proportional to the amount of turbulent energy accreted by the PNS. This relation is further supported via analytic arguments in Powell & Müller (2019). Here, we note (but do not show) that the total energy radiated increases with the effective mass throughout all evolutionary stages above 0.1 s. It is worth noting that the development of an explosion, as long as there is sustained accretion, can lead to GWs with an amplitude at a level higher than it would have been without an explosion (Powell & Müller 2019; Radice et al. 2019); therefore, the high GW amplitude seen in m0.95, and to some extent in m0.75, after the onset of explosion can be in part attributed to this as well. We return with an in-depth analysis and discussion regarding the source of the GWs in Section 3.4.

Zoom In Zoom Out Reset image size

As seen in previous studies (e.g., Murphy et al. 2009; Yakunin et al. 2015), we observe a prolate shape of the shock for the m0.95 model that is reflected in the GW signal by a shift in the positive direction at ∼0.6 s.

In Figure 3, we show the spectrograms where the typical, dominant PNS oscillation GW signal starts at ∼0.1 s and increases in frequency as the PNS cools and contracts. In all three spectrograms, we note a particular region that appears to lack emission, located at roughly constant frequency (∼1.2 kHz) throughout the simulation time. Morozova et al. (2018) were the first to address this "power gap," which they could not attribute to a numerical artifact, and thus do not rule out a physical origin. As the dominant oscillation mode evolves toward the gap region from below, an interaction appears to take place with another mode above the power gap, consistent with the avoided-crossing description in Morozova et al. (2018) and Sotani & Takiwaki (2020), resulting in the lower-frequency mode flattening out in frequency and dropping in amplitude, while the other higher-frequency mode rises in frequency and amplitude. The crossings between modes occur at ∼650, ∼850, and ∼1000 ms for the m0.95, m0.75, and m0.55 simulations, respectively. We report our further exploration of the power gap in Section 3.5.

Zoom In Zoom Out Reset image size

We show in Figure 4 the radius of the PNSs (top panel) and the three spectrograms overlaid on one another (bottom panel). We define the PNS radius as the maximum distance to where the density is 10¹¹ g cm⁻³, which leads to sporadic jumps in the evolution profile of the early-exploding run, m0.95. The dashed lines in the spectrogram plot are quadratic fits of the form At² + Bt + C to the collection of frequency points of most power (one extracted per STFT time step of 3.6 ms). The PNS radius and the frequency of the dominant oscillation mode are clearly effective mass dependent. It has been realized through asteroseismology studies (e.g., Morozova et al. 2018; Torres-Forné et al. 2019a) and works using semianalytical fits to the dominant frequency (e.g., Murphy et al. 2009; Cerdá-Durán et al. 2013; Müller et al. 2013; Pan et al. 2018; Pajkos et al. 2019; Powell & Müller 2019, 2020; Torres-Forné et al. 2019a; Sotani & Takiwaki 2020) that the mass and radius are the primary determinants of the dominant PNS oscillation frequency. In our simulations, the masses of the accreting PNSs remain equal until ∼0.5 s after bounce and differ by at most 4% by t_pb = 860 ms (which is the end of run m0.95) owing to explosion. However, the EOS dependence results in relatively large variations of the PNS radius, having evolved differently since bounce and differing by ∼15% between our two extreme cases at t_pb = 860 ms. Hence, the radius holds the dominant role in setting the frequency and accounts for the ∼36% difference in the peak frequency at 860 ms between runs m0.55 and m0.95. In summary, as expected owing to the role of the effective mass on the thermal properties outlined in Schneider et al. (2019b) and mentioned above, we find that a higher value of the effective mass parameter m^⋆ yields a more compact PNS, producing a higher oscillation frequency.

Zoom In Zoom Out Reset image size

Figure 5 shows the characteristic strain, h_char (Equation (13)), where we assume a distance of 10 kpc to the source. The frequency corresponding to the largest h_char, the so-called peak frequency, is shifted to higher frequencies and higher power when increasing the effective mass, i.e., when increasing the PNS compactness. The peak frequencies are located around 1–2 kHz, indicating that these amplitudes are due to the PNS oscillations. The power gap near 1.2 kHz is also clearly seen in this plot. Furthermore, because the peak frequency rises slower in the m0.55 run and the window over which h_char is calculated is limited to the maximum common run time, t_pb = 860 ms, the characteristic strain for this run above the power gap is lower than below the gap. This is not the case for the m0.75 and m0.95 runs, where the peak frequency crosses the power gap ≃200 ms earlier. Another feature we observe is a small shoulder of excess emission around ∼100 Hz, which could be explained by instabilities such as neutrino-driven convection and turbulent flows in the gain region, as well as the SASI, for which the GW radiation is typically emitted at such low frequencies (e.g., Cerdá-Durán et al. 2013; Pan et al. 2018; Pajkos et al. 2019). For the m0.95 run there is also an excess emission at frequencies ≲100 Hz due to the asymmetric explosion.

Zoom In Zoom Out Reset image size

The gray dashed line in Figure 5 represents a sensitivity curve of Advanced LIGO (Barsotti et al. 2018). The frequency and amplitude ranges are well within the bounds of the sensitivity curve. Caution should be taken when concluding their detectability since 2D simulations typically overpredict the signal strength compared to 3D simulations. Several 3D studies report on amplitudes decreasing by a factor of 10−20 in their 2D/3D comparisons (Andresen et al. 2017; O'Connor & Couch 2018b; Mezzacappa et al. 2020), which would make a detection at 10 kpc more marginal (see, e.g., Szczepanczyk et al. 2021). Although the amplitude and detection prospects may need revision in 3D simulations, we do not expect the evolution of the dominant PNS oscillation frequency to change significantly when going from 2D to 3D owing to an approximately preserved evolution of PNS mass and radius. 3D studies with 2D counterparts support this reasoning (e.g., O'Connor & Couch 2018a; Radice et al. 2019), meaning that the quadratic fits (recall Figure 4) and the EOS dependence are likely preserved for more realistic, 3D, CCSNe.

For our models with the highest and lowest values of the effective mass parameter, we run additional simulations using rotation profiles with central angular speeds of 1 and 2 rad s⁻¹ at the onset of collapse (see Equation (1)). Figure 6 shows the mean shock radius (top panel) and the neutrino heating (bottom panel). In contrast to the nonrotating m0.95 simulation, no shock revivals occur for either m0.95r1 or m0.95r2 during their ∼0.5 s simulation time.

Zoom In Zoom Out Reset image size

Due to centrifugal support, matter does not settle as deeply in the gravitational potential, releasing less gravitational binding energy during core collapse. This results in lower neutrino luminosity and slower contraction of the PNS (Summa et al. 2018; Westernacher-Schneider et al. 2019). When we increase the rotation rate, we observe a significant nonlinear decrease of the neutrino/antineutrino luminosities (as well as mean energies), which is reflected in the neutrino heating (Figure 6). Other processes may factor into the observed heating also: for the most rapidly rotating models, the centrifugally supported shocks are located at a large radius, increasing the region where neutrinos may deposit energy. On the other hand, when including rotation, convection is inhibited owing to a positive angular momentum gradient according to the Solberg–Høiland stability criterion (Endal & Sofia 1978; Fryer & Heger 2000), which is known to result in less prominent neutrino heating (Murphy et al. 2013).

Figure 7 shows the GW signal for the models varying in both the effective mass and rotation rate. Expectedly (e.g., Pajkos et al. 2019), the two bottom panels illustrate how the signal becomes increasingly muted with a higher rotation rate. As demonstrated by the two neighboring panels on top, and contrary to the nonrotating simulations, we do not find any effective mass dependence of the GW amplitude for our simulations that include rotation. This is of special note for the 1 rad s⁻¹ models, as they still exhibit an effective mass dependence in neutrino heating. The centrifugal support stabilizes the gain region to the extent that the different amounts of neutrino heating have little to no impact on the GW amplitude. Furthermore, it takes longer for the convection instabilities in the 2 rad s⁻¹ models to become fully nonlinear in the gain region, and thus these simulations exhibit an extended quiescent phase until ∼150–170 ms. Since these simulations are axisymmetric, the suppression of GW emission we see here due to rotation could be overestimated as compared to 3D. Furthermore, for sufficiently high rotation rates one expects nonaxisymmetric instabilities to give rise to potentially strong GWs (Scheidegger et al. 2010; Takiwaki & Kotake 2018; Andresen et al. 2019).

Zoom In Zoom Out Reset image size

Figure 8 shows the PNS radii (top panel) and the polynomial fits to the frequency points of most power (bottom panel). Correlated with the radius (which increases nonlinearly with higher rotation rate), the frequency of the dominant PNS oscillation mode decreases with increasing rotation. The effective mass dependence of this mode is still clear for the 1 rad s⁻¹ simulations, with rotation at this level only reducing the peak frequency by ∼10% from the respective nonrotating counterparts. For the 2 rad s⁻¹ models, due to the relatively little PNS excitation (and late onset thereof; see Figure 7), the frequency distribution of power is much broader and at a lower level. Nevertheless, using a linear fit to the frequency points of most power for these rapidly rotating models (see the bottom panel of Figure 8), we see a further reduction in the peak frequency and remark that the rapid rotation washes out any discernible frequency dependence that the slightly different PNS radii of the fast-rotating simulations would cause.

Zoom In Zoom Out Reset image size

The results above highlight how the compactness of the PNS at the time of formation and its cooling rate are the underlying factors of the GW signatures (under otherwise-similar conditions).

For the simulations where we investigate the effect of the incompressibility modulus, the PNS masses and radii are nearly identical across the models during their ∼0.6 s simulation time. This agrees with the spherically symmetric results of Schneider et al. (2019b). Expectedly, we see no systematic differences in the GW signatures between the three models, and any minor differences are likely stochastic. Figure 9 shows the characteristic strain as a function of frequency, which reflects the similarity in both amplitude and frequency between the three simulations. We note that the decrease in GW emission in the 1.0–1.2 kHz range is even more evident in Figure 9 now that the characteristic strains of the simulations h_char overlap. We discuss the incompressibility dependence of the gap further in Section 3.5.

Zoom In Zoom Out Reset image size

In a typical PNS there are two convectively stable regions, one in the inner core and one that composes the outer shell of the PNS that faces the turbulent post-shock region. According to canonical wave theory (Handler 2013; also see Gossan et al. 2020 for a direct application in CCSNe), apart from oscillatory p-modes (where pressure acts as a restoring force) at high frequency, these stable regions may also house g-modes at lower frequencies (where buoyancy acts as a restoring force). These excitations may emit GWs at select eigenmode frequencies when excited. Often, a buoyantly unstable region separates the core and the outer shell, the PNS convection zone. In 2D, many groups have attributed their GW signal to g-modes from the outer shell (Cerdá-Durán et al. 2013; Müller et al. 2013; Pan et al. 2018). A viable perturbing mechanism is accretion plumes impinging onto the PNS from the violent gain region above (Murphy et al. 2009), also recognized in 3D (O'Connor & Couch 2018b; Powell & Müller 2019; Radice et al. 2019; Andresen et al. 2021).

Adding to the possibilities of strong GW sources, Mezzacappa et al. (2020) conclude that after a few hundreds of milliseconds, PNS convection takes over as the dominant source of GW emission in their 3D study. Andresen et al. (2017) report that the dominant source of emission in their 3D runs originates from the bottom end of the convectively stable surface layer, driven by convective overshoot from the convection zone below rather than by accretion from above, and further show the same scenario for a 2D counterpart to a 3D run. It should also be mentioned that PNS convection has previously been recognized as a GW source and as a possible g-mode excitation mechanism in 2D (Marek et al. 2009; Murphy et al. 2009; Müller et al. 2013) and in 3D (Müller et al. 2012), although not established as a dominant source pre-explosion. Furthermore, Torres-Forné et al. (2019b; in 2D) point toward the stable inner core as the region responsible for the dominant emission in their study, and they support the hypothesis that the dominant emission is always a g-mode (Torres-Forné et al. 2019a). For Morozova et al. (2018) and Sotani & Takiwaki (2020) in 2D, and for Radice et al. (2019) in 3D, their dominant source of emission after roughly 0.4 s stems from the fundamental quadrupolar mode, having transitioned from a g-mode. In 3D, Powell & Müller (2019) attribute their high-frequency emission to g-modes, and in their complimentary work using other progenitors in Powell & Müller (2020), to the f-mode.

It is possible that the dominant GW emission is strongly model dependent and that the scenery of strong GW emission in CCSNe is as rich as the literature suggests. For example, with regard to the GW emission from the PNS, the results of Mezzacappa et al. (2020) resemble those of Andresen et al. (2017) and Andresen et al. (2019), where PNS convection is the source or driving force, while the findings in Andresen et al. (2021) support Radice et al. (2019) and Powell & Müller (2019), where nonradial motion in the post-shock region is recognized as the main forcing mechanism. Using analytic arguments, Powell & Müller (2019) speculate that the relative importance between these two excitation mechanisms may be determined by the difference in convective energy inside and outside the PNS, where the latter is closely tied together with the explodability of the model (such as the neutrino heating and convective Mach number). It is further possible that the situation is more generic and a stronger consensus will emerge as PNS asteroseismology techniques and classification methods become even more sophisticated and widely used. However, only a couple (Andresen et al. 2017; Mezzacappa et al. 2020) spatially decompose their GW signals from the actual simulation. Thus, to shed light on the excitation mechanism, we encourage the combination of self-consistent perturbation analysis and spatially resolving the emission. In this work, we refrain from firmly diagnosing our dominant emission with a particular mode; we leave this to future work. Instead, we attempt to localize the source of the GWs in our simulations while simultaneously highlighting two analytical procedures, hoping to provide value to future studies in unveiling the mechanism and source of GWs:

(1) In Appendix A we remark that caution needs to be taken when calculating the spatial distribution of GWs since the often-used analytic expression for the first time derivative of the quadrupole moment (Equation (10) in Cartesian and Equation (A5) in spherical coordinates) neglects a surface term when considering only a portion of the PNS. This term, of course, vanishes when calculating the GW emission from the whole star. In that appendix we show how, when instead numerically taking two time derivatives of the quadrupole moment itself, the apparent emission region changes. In particular, emission that was previously seemingly from the convection zone, where the surface terms are largest, decreases.

(2) In Appendix B we highlight that the location where the energy density associated with a particular mode (Morozova et al. 2018; Torres-Forné et al. 2018, 2019b; Sotani & Takiwaki 2020) is most concentrated, does not necessarily overlap with the region emitting the largest contribution of GWs owing to the excitation of that mode. This is for two reasons: (a) because significantly more mass is located farther out at larger radii, even though the energy density is lower, and (b) the added lever arm of the radius (squared) in the quadrupole moment. As a consequence of this analysis, we show that the radial profiles of the GW emission at any particular frequency, computed either from a perturbation analysis or dynamically from the simulation, agree remarkably well in the PNS itself.

With these in mind, we attempt to localize the source of the GW emission in our 2D simulations and place the results in the context of the literature. In Figure 10, we show the spatial distribution of the radial velocity (top panels) and radial profiles of the GW spectrogram (bottom panels) for t = 400 ms (left), t = 800 ms (middle), and t = 1000 ms (right) post-bounce for our baseline simulation, m0.75. The GW signals in these spectra are calculated in the post-processing of high-cadence data by taking two numerical time derivatives of the contribution to the zz component of the quadrupole moment (i.e., Equation (9)) contained in each spherical shell (and multiplying by the constants outside the double derivative in Equation (11)). We note that this is different than using a radially decomposed Equation (10) and taking one additional time derivative; we refer the reader to Appendix A for the important details. We obtain the frequency dependence by multiplying this radially decomposed GW signal by a 40 ms Hann window centered at each time specified above, and finally Fourier transform it. To the left of each spectrogram profile we show the normalized projected power spectral density. We note that both the profile and projection are logarithmic. In each spectrogram profile we also plot the turbulent energy flux (white solid lines for positive values and dashed for negative values). Similar to Equations (15)–(17) in Andresen et al. (2017), the turbulent energy flux, f_m , is calculated via

$$\begin{eqnarray}&&{f}_{m}=\langle \rho ^{\prime} {v}_{r}^{\prime} \rangle ,\end{eqnarray} \tag{ 15 }$$

where angled brackets denote an angular average, ρ denotes mass density, and v_r is the radial velocity of the fluid. Here the primed quantities are the deviations from the angle average, e.g., $\rho ^{\prime} =\,\rho -\langle \rho \rangle$. To reduce the stochastic noise of f_m in our 2D simulations, the curves in Figure 10 are an average of all the turbulent energy flux profiles obtained over the 40 ms window. We also include the Brunt–Väisälä frequency (black lines) following Mezzacappa et al. (2020),

$$\begin{eqnarray}&&{f}_{\mathrm{BV}}^{2}=-\displaystyle \frac{\partial \phi /\partial r}{{\left(2\pi \right)}^{2}\rho }\left(\displaystyle \frac{\partial \rho }{\partial s}{\left|{}_{P,{Y}_{\mathrm{lep}}}\displaystyle \frac{{ds}}{{dr}}+\displaystyle \frac{\partial \rho }{\partial {Y}_{\mathrm{lep}}}\right|}_{P,s}\displaystyle \frac{{{dY}}_{\mathrm{lep}}}{{dr}}\right),\end{eqnarray} \tag{ 16 }$$

where Y_lep = Y_e + Y_ν , and both P and s include contributions from the neutrinos as well. For these neutrino contributions, we assume a thermal distribution of neutrinos at densities higher than ∼ ρ = 10¹² g cm⁻³ (Kaplan et al. 2014). While in principle the convective region should have negative ${f}_{\mathrm{BV}}^{2}$, we only see this for the outer part of the convection zone for the t = 400 ms profile, and therefore we rely on the turbulent energy flux to define the convection zone. For each of the three times, we mark the location of ρ = 10¹¹, 10¹³, and 2 × 10¹⁴ g cm⁻³ with dashed vertical lines in the spectrograms and as contours in the radial velocity plots.

Zoom In Zoom Out Reset image size

In Figure 10, the convective region is traced very well by the width of the negative part of the turbulent energy flux profile. The negative sign here is understood by the nature of fluid motion in a convection zone. The heavier (denser) fluid is advected down and the lighter (less dense) fluid is buoyantly lifted up, in both cases $\rho ^{\prime} {v}_{r}^{\prime} \lt 0$. The turbulent energy flux will thus always be negative in the convective layer. The reverse is true in the overshoot layer above the PNS convection zone, as the overshoot plumes are denser than the local average, explaining the small positive peak in the turbulent energy flux just above the convection zone in Figure 10. Hence, as seen in the figure, the radius corresponding to ρ = 2 × 10¹⁴, 10¹³, and 10¹¹ g cm⁻³ roughly divides the PNS core, the PNS convection zone, the convectively stable PNS surface layer, and the gain region (and beyond). This is further highlighted by the density contours in the radial velocity plots in the top panel of Figure 10, where funnels of fast-moving matter in opposite directions are seen between ρ = 10¹³ g cm⁻³ and ρ = 2 × 10¹⁴ g cm⁻³, namely, convection sustained and driven by continued neutrino diffusion out of the core and the maintenance of the lepton gradients. For completeness, we show in Figure 11 properties of the PNS convection zone over the course of the simulations for each nonrotating simulation. Following the PNS radius trend in Figure 4, the PNS convection zone is located at a larger radius for simulations with a lower effective mass. While the smallest (largest) in volume, the m0.95 (m0.55) simulation shows the largest (smallest) PNS convection zone mass and kinetic energy at early times (t ≲ 300 ms). At late times, the mass contained in the convection zone levels out and the ordering inverts, with the m0.95 simulation having less mass than the others. This may be due to the explosion in the m0.95 simulation. At the level explored here, K_sat does not impact any of the PNS convection properties.

Zoom In Zoom Out Reset image size

Regarding the source of the GWs and the radially decomposed spectrograms in Figure 10, resolving a precise location in the star where the emission occurs is difficult. The emission region is extended, which is supported by the results shown in Appendix B, where we show, for these select times, the radial profile of the GW emission compared to the predictions from a perturbation analysis. However, the figure provides a strong indication that the dominant emission in our simulations stems from the convectively stable surface layer, with a tendency toward the bottom end of this layer, the convective overshoot region. Characteristic of g-modes, the ρ = 10¹³ g cm⁻³ line, which closely traces the location of the overshoot region, crosses the Brunt–Väisälä frequency in close proximity to the dominant emission region, especially at t = 800 ms (middle) and t = 1000 ms (right). While our analysis suggests that the dominant source of GW emission does not stem from the convection zone itself, we cannot rule out PNS convection as a forcing mechanism that accounts for a substantial fraction of the GWs generated. In fact, we think that this is plausible and supports the findings of Andresen et al. (2017), who also find significant emission in the PNS convective overshoot region. On the other hand, considering the correlation we see between the GW amplitudes and the neutrino heating (Sections 3.1 and 3.2), as well as the noncorrelation with the total PNS convection layer kinetic energy seen in Figure 11, we do not rule out accretion plumes as a strong excitation mechanism. Furthermore, via entropy field plots we observe clear funnels of low-entropy matter (from the gain region) striking the PNS at times coinciding with some of the largest amplitude spikes in the GW signal. Perhaps these simulations serve as an example where both mechanisms are at work.

Several groups report on GW amplitudes decreasing by a factor of 10−20 when going from 2D to 3D (Andresen et al. 2017; O'Connor & Couch 2018b; Mezzacappa et al. 2020). This is generally expected since the axisymmetric assumption leads to an artificial enhancement of coherent motion. Andresen et al. (2017) gather the contributions to the GW signal from different regions, such as the convection zone plus its overshoot layer (region A in their paper) and the PNS surface layer (region B in their paper). They report on a suppression of the GW signal by a factor of ∼10 in region A when going from 2D to 3D, and a higher such factor in region B. Andresen et al. (2017) discuss why a higher suppression of GWs excited by accretion plumes, compared to overshoot from below, is expected when going from 2D to 3D. In 2D, turbulence structures cascade in reverse (Xiao et al. 2009) and the Kelvin–Helmholtz instability is suppressed at the edge of supersonic downflows (Müller 2015); this may result in artificially high accumulation of the turbulent energy on large scales and higher downflow velocities onto the PNS. Furthermore, Andresen et al. (2017) demonstrate that in 2D the frequency range of the chaotic turbulent downflows from the gain region may have a stronger overlap with the natural g-mode frequency, such that the accretion mechanism may produce more resonant excitation in 2D compared to 3D. Speculatively, then, if both the convective overshoot and the accretion mechanisms drive the GWs in our simulations, and with these effects taken into account, 3D counterparts to our models may have convective overshoot as the main forcing of PNS oscillations; however, based on Appendix B, we do not expect the radial distributions to change dramatically.

We end our results section by exploring the power gap present in all of our simulations. We show the evolution of the characteristic strain for the m0.75 run in the top panel of Figure 12, where we have sequentially windowed 200 ms of the GW signal centered at the epoch we specify on the right-hand side. For a direct comparison, each but the first strain is shifted vertically using even multiplication factors on a logarithmic scale. The strains are plotted in color, with smoothed versions in black. This figure provides yet another perspective of the drift of the peak frequency corresponding to the dominant PNS oscillation mode that goes from ∼550 Hz at 280 ms to ∼1450 Hz at 880 ms. The power gap, which we mark by crosses at the local minima of the black lines, is located around 1.2 kHz and is particularly easy to see between 0.6 and 0.9 s, when modes are excited both above and below the gap. In the bottom panel, we plot the location of the minimum over many time steps for the three models varying in effective mass and for the three models varying in incompressibility modulus. Recall that m0.75 and k230 are the same, and our baseline simulation.

Zoom In Zoom Out Reset image size

Addressing the models varying only in the effective mass first, before ∼0.3 s, the gap locations are ordered in frequency with the value of the effective mass. For m0.55 and m0.75, which share the same general gap evolution, the frequency of the gap increases until it flattens out around ∼0.4 s and stays roughly constant thereafter. The gap evolution of the m0.95 simulation differs from the others since its frequency decreases with time.

The evolution of the gap does not appear to change when varying the incompressibility modulus. However, the quantitative value of the gap frequency does. It is located at lower frequencies for progressively higher values of K_sat. Recall that the k200 and k260 simulations are only carried out until ∼600 ms after bounce.

Morozova et al. (2018) speculate that one may view the PNS as a coupled system of the inner core and the outer convectively stable shell, mediated by the PNS convection zone, such that the power gap may originate via a repelling interaction between a trapped PNS core mode and a surface layer mode. Their perturbation analysis roughly reproduces the morphology of the gap when plotting the frequency of modes after having moved the outer boundary to the core surface (their Figure 8).

The different gap frequencies (bottom panel of Figure 12) in our study further hint toward an involvement of the inner core in producing the power gap. The incompressibility modulus has only a significant effect in the core, where the density is sufficiently high and the thermal component of the EOS is less dominant (Schneider et al. 2019b). Recall that the PNS radius defined at ρ = 10¹¹ g cm⁻³ and the PNS convection zone properties, which reflect the structure between ∼ 10¹³ and 10¹⁴ g cm⁻³, remain equal between the models varying in incompressibility modulus. Thus, the outer structure of the PNS alone is likely not the cause of the gap, nor does it set the frequency. However, the inner core structure may. Varying the incompressibility (effective mass) above and below baseline changes the central density by ∼4% (∼12%). Furthermore, the central density and the gap frequency share the same ordering between all the models. That is, the higher the central density, the higher the frequency of the gap. The m0.95 run is the only deviation to this trend owing to its deviation in gap morphology.

Related to the gap origin suggested by Morozova et al. (2018), other asteroseismology studies see examples of avoided crossings between modes (Torres-Forné et al. 2019b; Sotani & Takiwaki 2020). Torres-Forné et al. (2019b) introduce mode grouping/classification via the shape of the mode functions, rather than by the number of nodes. With this method, some of their earlier-predicted avoided crossings (when counting nodes) are replaced by plain crossings, shedding light on the exchange of properties when modes cross. If what we see in our simulations is an avoided crossing, perhaps like the analytical eigenmode solutions in Sotani & Takiwaki (2020), the location of the avoided crossing coincides with the gap (Figure 3). While this is likely no coincidence and would naturally explain the absence of emission in a small region in the spectrogram (where the modes interact), we have not untangled the process responsible for the power gap that persists during the entire simulation and postpone further investigation to future work. However, our preliminary perturbation analysis yields eigenmodes predicting the two main modes seen in the spectrograms and further reveals that the m0.95 model differs from the others by predicting another avoided crossing at an earlier epoch between less excited modes, a possible explanation for its unique gap morphology. From this preliminary analysis, we also see an eigenmode close to the gap whose energy is mostly confined in the inner core.

For completeness, we remark that the power gap is also evident (and marked with a black dashed line in the projected power spectral density) in the spatially decomposed spectrograms, see Figure 10 and Appendix A. Other regions lack emission too, which very roughly correspond to nodes from p-modes above the gap and g-modes below (although these are not necessarily eigenmodes but perhaps sporadic wave-like transient perturbations). These "g-mode" nodes appear to converge at the gap, which in principle could produce it.

Zoom In Zoom Out Reset image size

To finalize our discussion of the power gap, we speculate on the reason why some investigations in the literature do not seem to show evidence for this gap (e.g., Torres-Forné et al. 2019b; Mezzacappa et al. 2020; Jardine et al. 2021). Our results and those of Morozova et al. (2018) show clear evidence that the gap is tied to the physics at the highest densities in the inner (≲10 km) core of the PNS. One potential cause for this disparity is the use of a spherical, 1D, inner core, which prevents aspherical motions, or, in this case, potential mode interactions. Morozova et al. (2018), O'Connor & Couch (2018b), Radice et al. (2019), and this work all see the power gap (or the power gap can be seen in the publicly available data) and do not have such an inner core. The one exception we can find is the work of Pan et al. (2018), who use a very similar setup to our FLASH simulations, but whose published GW spectrograms show no sign of a power gap. On the other hand, several studies see a much richer set of isolated excited modes (Torres-Forné et al. 2019b) compared to what is seen in our study. It very well could be that aspects of our numerical setup are enhancing interactions between modes (including the interaction responsible for the power gap); further work still is needed to assess the power gap.