Exploring Fundamentally Three-dimensional Phenomena in High-fidelity Simulations of Core-collapse Supernovae

We make use of the FLASH hydrodynamics framework (Fryxell et al. 2000; Dubey et al. 2009) that we have outfitted for CCSNe in Couch (2013a, 2013b), Couch & O'Connor (2014), and more recently, have included both multidimensional, energy-dependent neutrino-radiation transport based on the moment formalism and effective general relativistic gravity in O'Connor & Couch (2018). For these simulations we use FLASH's unsplit hydrodynamic solver, which makes use of the piecewise parabolic method (Colella & Woodward 1984) and the hybrid HLLC Riemann solver, which reduces to HLLE in the presence of shocks. During each FLASH timestep, we solve for the new hydrodynamic state (the (n + 1) state) before the neutrino transport step. We refer the reader to previous works for a summary of the hydrodynamics and gravity methods (Couch 2013a, 2013b; Couch & O'Connor 2014; O'Connor & Couch 2018) and instead focus on a presentation of the multidimensional neutrino transport.

Neutrinos play a critical role in CCSNe. First and foremost, they provide a tremendous cooling channel for the newly formed PNS, they also are a source of heat in the shocked layers above the PNS, the so-called gain region. This heat source is thought to be critical for the development of the explosion. Following the evolution of the neutrinos from the core of the PNS out through the gain region requires a complex treatment of neutrino-radiation transport. Furthermore, since the opacity of the matter has a strong dependence on the neutrino energy, this transport must be done in an energy-dependent way. In FLASH, we have implemented a multidimensional, multispecies, energy-dependent two-moment (with an analytic closure) neutrino-radiation transport scheme. It is based on the work of O'Connor (2015), Shibata et al. (2011), and Cardall et al. (2013). We have outlined our FLASH implementation in great detail in O'Connor & Couch (2018). For the majority of our 3D simulations presented here, we ignore the velocity dependence of the neutrino transport equations (i.e., we set ${\boldsymbol{v}}$ and its derivatives to 0). We also ignore the gravitational redshift term that moves neutrinos between energy bins, although we keep the overall source term from the gravitational redshift. These approximations are due to not having fully implemented this physics prior to beginning our 3D simulations. For comparison purposes, we do perform two simulations with full velocity and gravitational redshift dependence, which were started after these improvements were made to our neutrino transport code. For completeness, we show the form of our moment evolution equations here and refer the reader to the Appendix of O'Connor & Couch (2018) for implementation details. The coordinate frame, energy-dependent, neutrino energy density (E; zeroth moment) evolution equation in Cartesian coordinates is given as

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}E+{\partial }_{i}[\alpha {F}^{i}]-{\partial }_{\nu }[\alpha \nu ({L}^{{ij}}{\partial }_{i}{v}_{j}+{F}^{i}{\partial }_{i}\phi )]\\ \quad =\,\alpha (W[\eta -{\kappa }_{a}J]-[{\kappa }_{a}+{\kappa }_{s}]{H}^{t}-{F}^{i}{\partial }_{i}\phi ),\end{array}\end{eqnarray} \tag{ 1 }$$

and the coordinate frame, energy-dependent, neutrino momentum density (Fⁱ; first moment) evolution equation is

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}{F}^{i}+{\partial }_{j}[\alpha {P}^{{ij}}]-{\partial }_{\nu }[\alpha \nu ({N}^{{ijk}}{\partial }_{j}{v}_{k}+{P}^{{ij}}{\partial }_{j}\phi )]\\ \quad =\,-\alpha ([{\kappa }_{s}+{\kappa }_{a}]{H}^{i}-W[\eta -{\kappa }_{a}J]{v}^{i}+E{\partial }_{i}\phi ),\end{array}\end{eqnarray} \tag{ 2 }$$

where in these equations, ${\boldsymbol{v}}$ is the matter velocity (with the associated Lorentz factor of W); $\alpha =\exp (\phi )$ is the lapse (with ϕ being the gravitational potential), ν and ∂_ν denote the neutrino energy and the energy-space derivative; κ_a, κ_s, and η are the matter absorption opacity, scattering opacity, and emissivity, respectively; P^ij is the second moment of the neutrino distribution function, that we use an analytic closure for; L^ij and N^ijk are higher moment tensors used in the energy-space flux determination; and J and Hⁱ are the fluid-frame zeroth and first moments that can be expressed in terms of E, Fⁱ, and P^ij,

$$\begin{eqnarray}&&J={W}^{2}[E-2{F}^{i}{v}_{i}+{v}^{i}{v}^{j}{P}_{{ij}}],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{H}^{t}={W}^{3}[-(E-{F}^{i}{v}_{i}){v}^{2}+{F}^{i}{v}_{i}-{v}_{i}{v}_{j}{P}^{{ij}}],\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{H}^{i}={W}^{3}[-(E-{F}^{j}{v}_{j}){v}^{i}+{F}^{i}-{v}_{j}{P}^{{ij}}].\end{eqnarray} \tag{ 5 }$$

These equations are energy and species dependent. We use 12 energy groups and three neutrino species (ν_e, ${\bar{\nu }}_{e}$, and ${\nu }_{x}\,=\{{\nu }_{\mu }+{\bar{\nu }}_{\mu }+{\nu }_{\tau }+{\bar{\nu }}_{\tau }\}$). This amounts to 144 evolution equations. Since we perform the spatial flux and energy flux calculations explicitly, these differential equations are only coupled within an energy group and spatial zone. We perform the resulting 4 × 4 matrix inversions analytically.

The neutrino-matter interaction coefficients, κ_a, κ_s, and η are also energy and species dependent. They also depend on the matter density, temperature, and electron fraction. We use NuLib (O'Connor 2015) to generate a three-dimensional table (in density, temperature, and electron fraction), which we then tri-linearly interpolate, for each energy group and species, on the fly using the $(n+1)$-state matter variables from the hydrodynamic step. For the underlying neutrino interactions, we use isotropic scattering on nucleons, alpha particles, and heavy nuclei following Bruenn (1985) and Burrows et al. (2006), with corrections for weak magnetism and recoil from Horowitz (2002). For emission and absorption processes involving electron-type neutrinos and antineutrinos, we employ charged current interactions on neutron, protons, and heavy nuclei, specifically,

$$\begin{eqnarray}&&n+{e}^{+}\leftrightarrow p+{\bar{\nu }}_{e},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&p+{e}^{-}\leftrightarrow n+{\nu }_{e},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&(A,Z)+{e}^{-}\leftrightarrow (A,Z-1)+{\nu }_{e},\end{eqnarray} \tag{ 8 }$$

from Bruenn (1985) using weak magnetism and recoil corrections from Horowitz (2002). Finally, for thermal pair processes we use both electron–positron annihilation and nucleon–nucleon bremsstrahlung following Burrows et al. (2006) as implemented in O'Connor (2015). We only include these processes for heavy-lepton neutrinos.

For our suite of 3D simulations, we adopt the 20 M_⊙, solar metallicity progenitor model from Farmer et al. (2016) produced using a non-equilibrium 204 isotope nuclear network with MESA (Paxton et al. 2011, 2013, 2015, 2018). We use the SFHo nuclear equation of state (EOS) (O'Connor & Ott 2010; Hempel et al. 2012; Steiner et al. 2013) everywhere. We simulate the collapse phase and early (first 15 ms) post-bounce phase using GR1D (O'Connor 2015) and then transition to FLASH. We use the same EOS table and neutrino interactions as well as the same description of gravity (i.e., using a general relativistic (GR) effective potential instead of full GR). We do interpolate from the 18 energy groups used in GR1D to the 12 energy groups used in the 3D FLASH simulations. We see a very smooth transition between the codes with little to no sign of transients that could impact the future evolution.

It is worth discussing the MESA 20 M_⊙ progenitor and comparing it to another model frequently used in the literature. The 20 M_⊙ MESA model (referred to as MESA20) is fairly similar to the 20 M_⊙ model from Woosley & Heger (2007) (referred to here as s20). The compactness (ξ_M; O'Connor & Ott 2011) of the two 20 M_⊙ models are also similar: ${\xi }_{1.75}^{\mathrm{MESA}20}\sim 0.69$ and ${\xi }_{1.75}^{{\rm{s}}20}\sim 0.75$. The silicon–oxygen interface, where a sharp drop in density occurs (∼50% is both models), is located at a radius of ∼2400 km and ∼2600 km for MESA20 and s20, respectively. The baryonic mass coordinates of these interfaces are ∼1.75 M_⊙ and ∼1.8 M_⊙, respectively. In summary, the similarity of these properties results in a qualitative and quantitatively similar evolution of the mass accretion rate onto the PNS throughout the post-bounce phase.

We simulate the post-bounce phase in 3D using a Cartesian grid with adaptive mesh refinement. Our smallest grid zones are ∼488 m. The adaptive mesh refinement would ensure that all of the post-shock region is fully refined, however, for our baseline simulations we limit the refinement so as to only maintain a minimum effective angular resolution of Δx/r ≲0.009 or Δx/r ≲ 053. This restriction leads to a reduction in resolution from Δx ∼ 488 m to Δx ∼ 1 km at r ∼ 108 km, and further reduction by a factor of 2 at r ∼ 216 km. We refer to this as our standard resolution. We also have lower resolution simulations (denoted by LR in the model name) where we still take the smallest Δx = 488 m, but only enforce Δx/r ≲ 0.015 or Δx/r ≲ 088. This leads to the first resolution decrement at r ∼ 65 km, and subsequent decrements at r ∼ 130 km, r ∼ 260 km, etc. Our baseline simulations are full 3D (4π), however we perform some simulations in octant symmetry (denoted by oct in the model name). For two of our simulations we use an improved neutrino transport that includes full velocity dependence. We denote these simulations with a v in the model name. Finally, in two of our simulations we introduce velocity perturbations into the silicon and oxygen shell at the time of mapping to FLASH (denoted by pert in the model name). We describe the details of these perturbations below. For the simulations without progenitor perturbations we do not impose any seed perturbations to the simulations, instead, we let the Cartesian grid seed deviations from spherical symmetry.

In total, we perform eight 3D simulations. Our flagship simulations, mesa20 and mesa20_pert, use our standard resolution. We also perform these simulations at lower resolution: mesa20_LR and mesa20_LR_pert. We include velocity dependence in mesa20_v_LR. Finally, we do three octant simulations mesa20_oct, mesa20_LR_oct, and mesa20_v_LR_oct. We also perform a collection of 1D and 2D simulations in order to address the dimensional dependence.

2.3.1. Progenitor Perturbations

Our suite of eight 3D simulations, which all use the same progenitor model, nuclear EOS, and neutrino opacities, span a number of interesting parameters, including resolution, geometry, inclusion of progenitor perturbations, and inclusion of neutrino-transport velocity dependence. We have done this in order to be able to test the sensitivity of our 3D simulations to each of these aspects. In this section we present an overview of all our of 3D results and discuss each of these conditions in turn.

We discuss each of 3D simulations below with the help of Figure 1 through Figure 5. In each of these figures we will show four basic quantities, which we describe here. (1) The mean shock radius r_sh; (2) the lateral kinetic energy in the gain region ${T}_{\mathrm{gain}}^{\mathrm{lat}}$, which is the sum of all non-radial kinetic energy in the gain region; (3) the ratio of the advection timescale through the gain region to the heating timescale in the gain region,

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{\mathrm{adv}}}{{\tau }_{\mathrm{heat}}}=\displaystyle \frac{{M}_{\mathrm{gain}}/\dot{M}}{| {E}_{\mathrm{gain}}| /{\dot{Q}}_{\mathrm{heat}}},\end{eqnarray} \tag{ 9 }$$

where M_gain is the mass of the gain region, $\dot{M}$ is the accretion rate outside the shock, E_gain is the energy of the matter in the gain region (gravitational + kinetic + internal [relative to free neutrons]), and ${\dot{Q}}_{\mathrm{heat}}$ is the rate of energy injection via neutrino heating; and (4) the heating rate in the gain region, ${\dot{Q}}_{\mathrm{heat}}$. The ratio, τ_adv/τ_heat, gives a quantitative measure of how close a simulation is to an explosion, especially when directly comparing models with the same underlying progenitor model and EOS. Empirically it has been found that when this ratio reaches 1 an explosion is expected (Marek et al. 2009; O'Connor & Couch 2018). Physically, τ_adv/τ_heat = 1 means that the time it takes to change the energy of the matter in the gain region by of order itself via neutrino heating is equal to the time it takes to accrete through the gain region. In this time the matter can be heated and unbound before settling onto the PNS.

Zoom In Zoom Out Reset image size

3.1.1. Impact of an Aspherical Progenitor

To begin, we discuss our two flagship simulations: mesa20 and mesa20_pert. These differ in the inclusion of progenitor perturbations in the silicon and oxygen shell as presented in Section 2.3. In Figure 1, we show the four basic quantities described above, r_sh (top left), ${T}_{\mathrm{gain}}^{\mathrm{lat}}$ (bottom left), τ_adv/τ_heat (top right), and ${\dot{Q}}_{\mathrm{heat}}$ (bottom right). Also in this figure we show the low-resolution variants of each of these simulations, mesa20_LR and mesa20_LR_pert, which will be primarily discussed in the following section. We show volume renderings at ∼150, ∼190, ∼220, and ∼260 ms (from left to right) in Figure 2 for both the mesa20 (top) and mesa20_pert (bottom) simulations to aid our discussion. Immediately we can see the impact of the aspherical progenitor in mesa20_pert and mesa20_LR_pert. Starting as early as ∼100 ms after bounce we see the first signs of the imposed asphericity in the silicon shell impacting the post-bounce dynamics. The amount of lateral kinetic energy in the aspherical models is diverging from and significantly exceeding the spherical models. The impact on the shock radius and the neutrino heating is not appreciable at this time. The mass accretion rate is high at this time, and the turbulent motions that do form are quickly accreted out of the gain region. It is not until later, ∼180 ms, that we see an impact on these other measures. Since this time is well before when the interface between the silicon and oxygen shells accretes through the shock (at ∼250 ms), the behavior seen at this time is still due to the perturbations imposed on the silicon shell. At this time, we see a dramatic rise in the lateral kinetic energy (∼7 × 10⁴⁸ erg s⁻¹ for mesa20_pert compared to ∼2 × 10⁴⁸ erg s⁻¹ for mesa20 at 220 ms), the shock recession temporarily ceases, and the neutrino heating rate is enhanced, 7.5 × 10⁵¹ erg s⁻¹ in mesa20_pert compared to $5\times {10}^{51}$ erg s⁻¹ in mesa20 at ∼220 ms. We note that the increase in the lateral kinetic energy is not because there is significantly more lateral kinetic energy accreting through the shock in the aspherical models (there is not), but rather that the non-radial velocities entering the post-shock region are much more effective at seeding the convective instabilities that drive turbulence.

Zoom In Zoom Out Reset image size

The quantitative impact of the progenitor asphericity is clear in the graph of the τ_adv/τ_heat. During the phase where the silicon shell perturbations are accreting through the gain region both the gain region itself is larger (giving a long advection time) and the neutrino heating rate is higher (giving a smaller heating timescale). We note that the longer advection time dominates the contribution to the larger τ_adv/τ_heat. After the silicon–oxygen shell interface accretes through the shock the perturbations from the silicon shell cease and are quickly accreted through the gain region and onto the PNS. We find that the perturbations imposed on the oxygen shell do not have the same qualitative impact as those from the silicon shell. We see no significant excess of lateral kinetic energy during this phase and the evolution of the mesa20_pert simulation begins to qualitatively match that of mesa20. The shock radius quickly recedes, and the lateral kinetic energy and the neutrino heating drop. τ_adv/τ_heat also drops from its peak value of ∼0.7 to ∼0.45.

3.1.2. Impact of Reduced Resolution

3.1.3. Impact of Octant Geometry

In Figure 3, we show the effect of imposing octant symmetry on our simulations. We show our main simulation, mesa20 (blue), as well as its low-resolution variant mesa20_LR (green). For each of these, we show variants, mesa20_oct (red) and mesa20_LR_oct (pink), where we restrict the motion to the first octant by imposing reflecting boundary conditions on each of the main axis. In the early evolution, prior to ∼100 ms, we see very little impact of the octant symmetry. Until ∼200 ms, difference that do arise are smaller than the difference observed between the standard and low-resolution simulations. For mesa20_oct, we see a small decrease in the lateral kinetic energy and well as the neutrino heating, shock radius, and τ_adv/τ_heat. This is similar to the effect seen in Roberts et al. (2016) who also perform octant simulations. Since the octant symmetry prevents the growth of the largest scale modes that help support the shock, after ∼100 ms both the octant simulations of Roberts et al. (2016) and the octant simulations presented here have marginally smaller average shock radii and neutrino heating. For both the standard and low-resolution simulations we note that the increased noise for the quantities from the octant simulations in Figure 3 is due to the smaller number of zones that we average over (by a factor of 8).

Zoom In Zoom Out Reset image size

After 200 ms, the difference between the full 3D simulations and the octant variants becomes more interesting. Due to the octant symmetry, global modes are suppressed. As a result, our octant simulations behave much like failed 1D models. While the full 3D simulations experience repeated episodes of SASI growth accompanied by excursions of the mean shock radius, the shock radii in the octant models are rapidly receding. The heating is also notably reduced, at times up to ∼25% lower, from the full 3D models, consistent with the smaller shock radii. From the ratio of τ_adv/τ_heat, the octant simulations are quantitatively the farthest from explosion.

3.1.4. Impact of Velocity-dependent Transport

We also perform simulations using a velocity-dependent variant of our neutrino transport. Unlike the other simulations presented here, these simulations fully account for velocity-dependent effects in the transport, gravitational redshift of the neutrinos as they stream out of the gravitational well, as well as advection of neutrinos with the fluid flow when they are trapped. We perform two low-resolution simulations with velocity dependence, one in full 3D, mesa20_v_LR and one in octant symmetry mesa20_v_LR_oct. We show the results of these simulations in Figure 4, where we also include the non-velocity-dependent simulations mesa20_LR and mesa20_LR_oct. We find that including velocity dependence has several effects on the evolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the velocity-dependent results, there appears to be slightly less overall heating (∼10%) in the gain region before ∼200 ms after bounce. This is due to, unfortunately, slightly different definitions of the energy exchange term in our two simulations that cannot be corrected with post-processing, but is purely a bookkeeping difference. In the velocity-dependent simulations, our matter exchange term is recorded as the rate of total energy exchange with the matter. While in the non-velocity-dependent simulations we record the rate of internal energy exchange with the matter. The difference between these two rates is the energy change due to momentum absorption of the neutrinos, which in the case of the gain region, acts to reduce the kinetic energy of the matter (hence why this rate is lower). We note that in both cases we properly take into account the energy exchange due to momentum absorption for the hydrodynamic source term. In 1D simulations we observe that the heating rate derived from the rate of total energy exchange in both the velocity-dependent and non-velocity-dependent simulations are very similar. This similarity is naively not what one would expect, but is a coincidental cancellation of two effects that seemingly equally impact the amount of neutrino heating with opposite signs. The first effect is from the gravitational redshift, which reduces the energy of the neutrinos, hence the effectiveness of neutrino capture in the gain region. A competing effect is that the background flow of the material in the gain region is against the neutrino flow so that in the frame of the fluid, the neutrinos are boosted to higher energies.

One of the more visible impacts of the added velocity dependence seen in Figure 4 is the evolution of the shock radius. With velocity dependence, the mean shock radii of the simulations mesa20_v_LR and mesa20_v_LR_oct reach ∼10 km further out at ∼100 ms after bounce and then recede more slowly such that after ∼250 ms they are ∼30 km (or ∼30%) further out compared to mesa20_LR and mesa20_LR_oct. This larger shock radius coincides with an increased PNS radius. Since such a dramatic difference in both the PNS radius and the shock radius is not seen in our 1D models of this progenitor, we attribute the difference, as least in part, to the presence of strong PNS convection, which is effectively supporting the PNS and allowing it to maintain a larger radius and therefore a larger shock radius. This effect, in relation to its impact on the heavy-lepton neutrino luminosities, has been noted in several recent multidimensional studies of CCSNe, including O'Connor & Couch (2018) and Radice et al. (2017), but note that other multidimensional simulations see a smaller impact, e.g., Buras et al. (2006). We discuss this further when we discuss the neutrino luminosities from our simulations in Section 3.3. Related to this is the observation that in the one full 3D simulation with velocity dependence we see no presence of the SASI. In part, we attribute this to the large shock radius discussed above. A consequence, as can be seen from Figure 4, after 200 ms after bounce there is a larger (≳30% more) turbulent kinetic energy present in the gain region of mesa20_v_LR when compared to mesa20_LR. This can lead to a suppression of the growth of the SASI (e.g., Foglizzo et al. 2007; Fernández et al. 2014).

3.1.5. Impact of Dimensionality

We see characteristic signs of the SASI in most of our full 3D models. Using the spherical harmonic decomposition analysis presented in Couch & O'Connor (2014), we compute the spherical harmonic components of the shock radius decomposition as a function of time for all of the full 3D simulations, mesa20, mesa20_LR, mesa20_pert, mesa20_LR_pert, and mesa20_v_LR. Based on the shock decomposition, the SASI is clearly present in models mesa20, mesa20_LR, and mesa20_pert, to a lesser extent in model mesa20_LR_pert, and not discernable in model mesa20_v_LR. In Figure 6, we show the individual components of the ℓ = 1 spherical harmonic, ${a}_{1-1}/{a}_{00}=\langle {z}_{\mathrm{sh}}\rangle /\langle {r}_{\mathrm{sh}}\rangle$, ${a}_{10}/{a}_{00}=\langle {y}_{\mathrm{sh}}\rangle /\langle {r}_{\mathrm{sh}}\rangle$, and ${a}_{11}/{a}_{00}=\langle {x}_{\mathrm{sh}}\rangle /\langle {r}_{\mathrm{sh}}\rangle$, where $\langle {x}_{i,\mathrm{sh}}\rangle$ denotes the mean value of the Cartesian coordinate x_i = {x, y, z} of the shock front and $\langle {r}_{\mathrm{sh}}\rangle$ is the mean radius of the shock front. We also show the total relative power in the ℓ = 1 mode, $| {\bar{a}}_{1}| \sim \sqrt{{\sum }_{m}{a}_{1\,m}^{2}}/{a}_{00}$ as the solid black line. The peak amplitudes reach ∼0.15–0.2 for the strongest three models, and ∼0.05 in mesa20_LR_pert. Most of the SASI-dominated phases are predominantly spiral modes where one expects a constant total power rather than for a sloshing mode where one would expect an oscillatory total power (at twice the frequency). At times there are clear, but small, oscillatory features on top of the total power (at twice the spiral mode frequency) that are due to a relatively smaller sloshing mode component. For example, from ∼250–325 ms in model mesa20_LR, as discussed in detail below, as a quintessential example.

Zoom In Zoom Out Reset image size

Spiral SASI modes, like the one seen from ∼250–360 ms in model mesa20_LR, redistribute angular momentum so that the PNS is counter-rotating with the spiral mode, even in zero net-angular momentum configurations (Blondin & Shaw 2007; Blondin & Tonry 2007; Hanke et al. 2013; Guilet & Fernandez 2014). We observe this redistribution in our models. In Figure 7, we show the enclosed angular momentum as a function radius for 18 times between 240 and 410 ms, in 5 ms increments. At all times, there is net-angular momentum inside the PNS (from ∼10–25 km) due to PNS convection. However, for the most part, this cancels itself out near the edge of the PNS convection zone. Outside of ∼25 km, at early times (yellows and light greens) there is little net enclosed angular momentum. As the SASI spiral mode grows, starting around ∼280 ms, we begin to see the redistribution of this angular momentum (dark greens toward black). By ∼360 ms (black dashed line), where there is a peak of the ℓ = 1 power, $| {\bar{a}}_{1}|$, the angular momentum redistribution is reaching a maximum. Inside ∼40 km, the PNS is counter-rotating against the spiral mode with an enclosed angular momentum of ∼1.5 ×10⁴⁵ g cm² s⁻¹. There is an equal, but opposite amount trapped in the SASI spiral wave outside of ∼40 km, so that the net-angular momentum of the system is zero. After ∼360 ms, we see the spiral SASI rapidly disrupts. The angular momentum in the gain region accretes onto the PNS rather than being sustained in the gain region. This process takes place over an accretion timescale ${M}_{\mathrm{gain}}/\dot{M}\sim 30$ ms. The amount of kinetic energy that is ultimately dissipated at the surface of the PNS during this time is ∼10⁴⁸ erg. Even if this was to be completely transformed into heat over ∼30 ms, the additional heating would be much less than the steady-state neutrino heating at this time.

Zoom In Zoom Out Reset image size

Complementary to Figure 7, we show slices along the SASI plane in Figure 8 (animation available online showing the spiral mode buildup and then the destruction; otherwise we show the frame at ∼360 ms) of the angular velocity of the material in this plane. Positive values (purples) are rotating counterclockwise while negative values (greens) are rotating clockwise. The spiral SASI wave is rotating counterclockwise, the triple point is located near the right side of the figure. Zones that are tagged as shocks are shown as transparent gray. In the core we see the buildup of clockwise rotating material (green) on the surface of the PNS, while the bulk of the gain region is rotating counterclockwise (purple).

Just prior to ∼360 ms when the spiral SASI mode peaks in amplitude, at around ∼350 ms, the shock starts a rapid expansion (see Figure 1). This is coincident with an increase in the lateral kinetic energy and an increase in the neutrino heating, likely as a result of the shock expansion. After this time the spiral SASI wave is disrupted. These dynamics suggest that the SASI mode becomes highly nonlinear, disrupts, and subsequently leads to a transient reenergization of the shock. This is quantitatively demonstrated via the increase in the ratio τ_adv/τ_heat. Clearly, this transient reenergization is not enough to launch as explosion, but it shows a potential way for the SASI to lead to an explosion. This phenomena is seen in other simulations present in this work (e.g., the peaks of the SASI mode in mesa20 at ∼400 and ∼470 ms, and mesa20_pert at ∼370 and ∼510 ms). It is important to note that not all transient shock expansions seen in Figure 1 are followed by a collapsed SASI wave, nor is it the case that all SASI wave collapses follow a rapid shock expansion, yet there does seem to be a high correlation between these phenomena.

Provided an explosion was launched during such an event, it may well be the case that the PNS retains the distributed angular momentum after the explosion as suggested in (Blondin & Shaw 2007; Blondin & Tonry 2007; Guilet & Fernandez 2014). Assuming the neutron star ends up with a total angular momentum of 1.5 × 10⁴⁶ g cm² s⁻¹, and a gravitational mass of ∼1.6 M_⊙ (based on the baryonic mass of the PNS at 360 ms of ∼1.8 M_⊙ and the SFHo EOS) and therefore a moment of inertia of ∼1.7 × 10⁴⁵ g cm², the final period would be ω = 8.8 rad s⁻¹, or ∼700 ms. This is roughly a factor of 5 faster than the prediction of Guilet & Fernandez (2014),

$$\begin{eqnarray}\begin{array}{rcl}P & \sim & 290{I}_{45}\left(\displaystyle \frac{10}{\kappa }\right)\left(\displaystyle \frac{{P}_{\mathrm{SASI}}}{50\,\mathrm{ms}}\right)\left(\displaystyle \frac{120\,\mathrm{km}}{{r}_{\mathrm{sh}}-{r}_{* }}\right)\left(\displaystyle \frac{{v}_{\mathrm{sh}}}{3000\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\\ & & \times \,\left(\displaystyle \frac{0.3{M}_{\odot }{{\rm{s}}}^{-1}}{\dot{M}}\right){\left(\displaystyle \frac{150\mathrm{km}}{{r}_{\mathrm{sh}}}\right)}^{2}{\left(\displaystyle \frac{{r}_{\mathrm{sh}}}{3{\rm{\Delta }}r}\right)}^{2}\,\mathrm{ms},\end{array}\end{eqnarray} \tag{ 10 }$$

where for our case, κ ∼ 9, P_SASI ∼ 13 ms, r_sh ∼ 110 km, r_* ∼ 40 km, v_sh ∼ 8000 km s⁻¹, M_⊙ ∼ 0.3, and ${\rm{\Delta }}r/{r}_{\mathrm{sh}}\,\sim \sqrt{2}| {\bar{a}}_{1}| \sim 0.2$, leading to a prediction of P ∼ 3500 ms. This difference could be due to any number of differences between the idealized models of Guilet & Fernandez (2014) and the work here. In particular, our SASI spiral wave is confined to a much smaller radius as the mean shock radius is only 110 km, compared to the typical value of 150 km in Guilet & Fernandez (2014). Also, the more massive PNS (1.8 M_⊙), and the smaller mean shock radius, gives a larger typical radial velocity behind the shock (8000 km s⁻¹ compared to 3000 km s⁻¹. There could be nonlinear effects associated with the large deviations from the assumed values and scalings in Guilet & Fernandez (2014).

We show the neutrino luminosities produced by each one of our five full 3D simulations in Figure 9. For clarity, we divide them into our standard resolution (left) and the low resolution (right). The solid lines show the electron neutrino luminosity, the dashed lines denote electron antineutrino luminosity, and the dashed-dotted lines show the luminosity of a single heavy-lepton neutrino species. The simulations begin at 15 ms after bounce. In the beginning of the non-velocity-dependent simulations there is a sharp spike in the luminosities due to the transition from the velocity-dependent transport in GR1D and the transport in FLASH. This spike is much smaller in the velocity-dependent FLASH simulations. Since the underlying progenitor is the same for each simulation we do not expect large variations in the predicted neutrino signals between simulations. That being said, we do see several noteworthy differences that we will explicitly mention. As the system evolves we see the rise of both the ${\bar{\nu }}_{e}$ and ν_x luminosities and the fall of the ν_e luminosity as the remaining iron core accretes onto the PNS. The electron-type luminosities plateau at ∼65 × 10⁵¹ erg s⁻¹, while the ν_x luminosity plateaus at ∼33 × 10⁵¹ erg s⁻¹. At a post-bounce time of ∼220 ms we begin to see the dramatic drop in luminosity that is a result of the mass accretion rate drop associated with the accretion of the silicon–oxygen interface through the shock. We do notice a difference between the simulations with and without the imposed progenitor perturbations. The mesa20_pert and mesa20_LR_pert simulations (shown in orange) have a shallower drop in luminosity. It starts earlier and ends later. This is not due to differences in the mass accretion rate onto the shock, but rather due to a longer advection time through the gain region and therefore a slower (for a short time) mass accretion rate into the cooling region. The longer advection time is due to the additional support provided to the shock by the turbulence motions seeded by the progenitor perturbations. We see a similar phenomena at later times, when collapsing spiral SASI waves trigger increased turbulence activity (see Section 3.2), increased advection time, and small drops in the neutrino luminosity. For example, at ∼370 ms in mesa20_LR, ∼400 ms in mesa20_pert, and ∼410 and ∼480 ms in mesa20.

Zoom In Zoom Out Reset image size

The luminosities presented in Figure 9 are spherically averaged. However, since our transport code is multidimensional, we can easily examine the variation of the neutrino signal over the emission direction. We perform the same spherical harmonic decomposition on the individual neutrino signals (extracted at 500 km) as we applied on the shock surface. However, for neutrinos it is convenient to consider the decomposition in the form of (Tamborra et al. 2014; Melson 2016)

$$\begin{eqnarray}&&{L}_{\nu }\sim {L}_{\nu }^{\mathrm{monopole}}+{L}_{\nu }^{\mathrm{dipole}}\cos \theta ,\end{eqnarray} \tag{ 11 }$$

where the monopole component is related to the ℓ = 0 decomposition, ${L}_{\nu }^{\mathrm{monopole}}={({a}_{00}^{2})}^{1/2}$, the dipole component is related to the ℓ = 1 decomposition, ${L}_{\nu }^{\mathrm{dipole}}=3\times {({\sum }_{i=-1}^{1}{a}_{1i}^{2})}^{1/2}$, and $\cos (\theta )$ is the cosine of the angle between the observer and the direction of the dipole at any given time. We do find a nonspherical (i.e., dipole) component that is highly correlated with the SASI motions. For ν_e and ${\bar{\nu }}_{e}$ we find ${L}_{{\nu }_{e}/{\bar{\nu }}_{e}}^{\mathrm{dipole}}=3| {\bar{a}}_{1}{| }_{{\nu }_{e}/{\bar{\nu }}_{e}}\sim 0.03{L}_{{\nu }_{e}/{\bar{\nu }}_{e}}^{\mathrm{monopole}}$. The ν_x signal also has a directional variation. However, it is roughly four times smaller, i.e., ${L}_{{\nu }_{x}}^{\mathrm{dipole}}=3| {\bar{a}}_{1}{| }_{{\nu }_{x}}\sim 0.008{L}_{{\nu }_{x}}^{\mathrm{monopole}}$. As the spiral SASI wave rotates around, so to does the dipole component of the neutrino luminosities. Over the course of one SASI period, the neutrino luminosities as seen by an observer viewing the spiral SASI wave edge on would have a trough-peak variation of ∼6% (∼1.5%) for ν_e and ${\bar{\nu }}_{e}$ (ν_x). This is smaller than seen in Tamborra et al. (2013), where peak-trough variations of order ∼25% were seen for some spiral SASI waves. The phase of the neutrino modulation is interesting. After taking into account the light travel time from the PNS where the neutrinos are emitted and at 500 km where they are measured, we find that the ν_e and ${\bar{\nu }}_{e}$ luminosities peak in the direction opposite the average shock position, i.e., they lag the direction determined by the $\langle {x}_{i,\mathrm{sh}}\rangle$ by 180^◦. The ν_x luminosities generally peak at a different time than the ν_e and ${\bar{\nu }}_{e}$ neutrinos. We attribute this to a variation in the location of the peak values of the density and temperature within the spiral SASI plane as a function of radius. The electron-type neutrinos and antineutrinos are emitted at a much larger radius than the heavy-lepton neutrinos. The accretion waves take longer to reach the deeper radii where the heavy-lepton neutrinos are emitted. The exact phase difference between the neutrino types smoothly varies with time.

Finally, we remark on the largest difference seen in Figure 9. As discussed in Section 3.1.4, when we include velocity dependence into our simulations, i.e., in the mesa20_v_LR and mesa20_v_LR_oct simulations, we see an increased level of PNS convection. This dramatically increases (by upwards of ∼30% compared to the non-velocity-dependent simulations) the heavy-lepton neutrino luminosity starting from ∼100 ms after bounce. The convection brings hotter material from deeper down to larger radii where emitted neutrinos can more easily escape. This also increases the radii of the heavy-lepton emission region, giving a larger effective area and therefore a larger luminosity. As mentioned in Section 3.1.4, this effect has been noted in several recent multidimensional studies of CCSNe, including O'Connor & Couch (2018) and Radice et al. (2017). There is tension regarding the total expected enhancement due to multidimensional effects as we note that other multidimensional simulations see a smaller impact (Buras et al. 2006).

The LESA phenomena was first reported in Tamborra et al. (2014). The LESA is characterized by a spatial asymmetry in the total electron number emission via a strong dipole component. The LESA can be described with a monopole and dipole component,

$$\begin{eqnarray}&&{N}_{{\nu }_{e}}-{N}_{{\bar{\nu }}_{e}}={A}_{\mathrm{monopole}}+{A}_{\mathrm{dipole}}\cos \theta ,\end{eqnarray} \tag{ 12 }$$

where $\cos (\theta )=\hat{r}\cdot {\hat{r}}_{\mathrm{dipole}}^{\nu }$. In some cases reported in Tamborra et al. (2014), the dipole component was comparable to or larger than the monopole component. The consequence is that in the direction of $-{\hat{r}}_{\mathrm{dipole}}^{\nu }$, the net lepton-number emission was negative, i.e., more electron antineutrinos were emitted then electron neutrinos. Such a situation may have dramatic consequences for neutrino oscillations, neutrino detection (parameter inferences), and nucleosynthesis, to name a few.

Since Tamborra et al. (2014), the LESA has been reported in all 3D models from the Garching Collaboration using the PROMETHEUS-VERTEX code (Janka et al. 2016) but has never been conclusively demonstrated by an independent source. In addition to being only observed by one simulation code, another large criticism of the LESA instability is that the neutrino transport scheme employed made use of the ray-by-ray approximation. Our neutrino transport scheme is completely independent of PROMETHEUS-VERTEX. Importantly, it does not make the ray-by-ray approximation but rather solves the multidimensional transport of the neutrinos directly. For these reasons, a search for LESA in our simulations is warranted.

We search our simulations for evidence of LESA. We extract from our simulations the net lepton-number flux through a sphere with radius 500 km. In Figure 10, we show the relative net lepton flux $({N}_{{\nu }_{e}}-{N}_{{\bar{\nu }}_{e}})/\langle {N}_{{\nu }_{e}}-{N}_{{\bar{\nu }}_{e}}{\rangle }_{{\rm{\Omega }}}$ from the mesa20_v_LR simulation for two times (330 and 440 ms after bounce). The presence of a dipole is clear, located at ϕ ∼ 90° (∼15°) and θ ∼ −20° (∼−15°) for t_pb = 330 ms (440 ms). At 1 ms time intervals, we decompose the net lepton-number flux (${N}_{{\nu }_{e}}-{N}_{{\bar{\nu }}_{e}}$) sky map into spherical harmonics. Recall from Section 3.3, this decomposition gives both ${A}_{\mathrm{monopole}}={({a}_{00}^{2})}^{1/2}$ and ${A}_{\mathrm{dipole}}=3\times {({\sum }_{i=-1}^{1}{a}_{1i}^{2})}^{1/2}$, where a_lm are the spherical harmonic decomposition components (Couch & O'Connor 2014). This notation is such that if the monopole and dipole term are equal in amplitude (and the dipole term is positive), then an observer located along the dipole direction sees twice the average net lepton-number flux and an observer along the anti-dipole direction sees zero net lepton-number flux. In model mesa20_v_LR we see a strong presence of LESA, as expected based on the sky maps seen in Figure 10. In the top panel of Figure 11, we show A_monopole and A_dipole for model mesa20_v_LR. The dipole component begins to grow at ∼200 ms after bounce and reaches an appreciable fraction of the monopole component (at times ∼50%) by ∼300–400 ms after bounce. In addition to having a large ratio of the dipole to monopole component, we also see that the LESA is stable (but migratory) in its position. In the bottom panel of Figure 11, we show, in blue, the θ (solid) and ϕ (dashed) values of the dipole direction, ${\hat{r}}_{\mathrm{dipole}}^{\nu }$, determined from the individual spherical harmonic components a_1m of the net lepton-number flux. Prior to 200 ms the dipole strength is very low and therefore the dipole direction is not well constrained. We exclude this time from the figure for clarity. We also show the simultaneous position of an asymmetry in the electron fraction near the PNS convection zone, ${\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}$. We evaluate the direction of this dipole between the radii of 25 and 30 km (where the PNS convection is strongest) as

$$\begin{eqnarray}&&{\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}=\langle {x}_{i}{Y}_{e}\rangle /{\left(\displaystyle \sum _{j=\mathrm{1,2,3}}\langle {x}_{j}{Y}_{e}{\rangle }^{2}\right)}^{1/2}.\end{eqnarray} \tag{ 13 }$$

The direction of this dipole is remarkably similar to ${\hat{r}}_{\mathrm{dipole}}^{\nu }$. In fact, as shown in the inset of the bottom panel in Figure 11, the dipole direction in Y_e slightly precedes the dipole direction in the net lepton-number flux. This lag time is roughly ∼2 ms, and is due to the light travel time between the location of the Y_e dipole and the sphere at which the neutrino number fluxes are measured. To make this more clear, in the inset we show both the true (solid orange line) and a shifted version (dotted orange line; shifted by 500 km/c = 1.67 ms) of the ${\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}$. The alignment of these dipoles is not unexpected. The asymmetry in Y_e is directly responsible for the asymmetry in the net lepton-number flux. The matter in the direction of ${\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}$ is on average more electron rich and therefore emits relatively more electron neutrinos compared to electron antineutrinos. Unlike the asymmetry in the net lepton number, the asymmetry in the electron fraction can be seen earlier than 200 ms. However, its position is less stable at this time. Although not shown, we also see a small dipole component to the total ν_x flux that is anti-aligned with the dipole direction (i.e., anti-aligned with the excess flux of the ν_es and aligned with the excess flux of the ${\bar{\nu }}_{e}$ s). This amplitude of this ν_x dipole component is ∼25% of the amplitude of the ν_e and ${\bar{\nu }}_{e}$ dipole components, which is in agreement with Tamborra et al. (2014).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The cause of the asymmetry in Y_e is not yet clear. Tamborra et al. (2014) propose a self-sustaining mechanism where the increase of electron antineutrinos in the hemisphere of $-{\hat{r}}_{\mathrm{dipole}}^{\nu }$ (and the decrease in the hemisphere of ${\hat{r}}_{\mathrm{dipole}}^{\nu }$) leads to relatively stronger neutrino heating, a larger shock radius, and consequently funneled accretion onto the opposite hemisphere (i.e., near the direction of ${\hat{r}}_{\mathrm{dipole}}^{\nu }$). This accretion of predominately electron-rich material sources the imbalance of the lepton emission, and according to Tamborra et al. (2014) excites the enhanced PNS convection. Our simulations lend support to this idea. We show in Figure 12 the relative difference in the shock radius, neutrino heating, and mass accretion rate (measured between 55 and 65 km) as measured in the two hemispheres defined by the direction determined via the Y_e asymmetry. Specifically, we show

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{{X}^{+}-{X}^{-}}{\tfrac{1}{2}({X}^{+}+{X}^{-})},\end{eqnarray} \tag{ 14 }$$

where the + denotes the hemisphere defined with ${\boldsymbol{r}}\cdot {\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}\gt 0$ and − denotes the hemisphere where ${\boldsymbol{r}}\cdot {\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}\lt 0$. We see a significant (upward of ∼30% at times) higher mass accretion rate in the hemisphere aligned with ${\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}$, and therefore with ${\hat{r}}_{\mathrm{dipole}}^{\nu }$. Near the time when the LESA is strongest (between ∼300 and ∼450 ms), we also see lower average shock radii and neutrino heating rates (upward of ∼2%–3%) in the hemisphere aligned with ${\hat{r}}_{\mathrm{dipole}}^{{{\rm{Y}}}_{{\rm{e}}}}$. The sign of each of these relative ratios, agree with the observations of Tamborra et al. (2014). Furthermore, for the most similar progenitor used in Tamborra et al. (2014), s20, the magnitudes of these differences roughly correspond to each other. Interestingly, the mass accretion rate asymmetry begins to form early (before 200 ms) before significant asymmetry forms in the shock radius, neutrino heating, or even net lepton-number emission. This suggest that although asymmetric heating via a LESA dipole may funnel mass accretion into the opposite hemisphere and thereby closing the self-sustaining cycle, there may be other mechanisms that cause the mass accretion asymmetry in the first place in order to initiate the growth of the LESA dipole.

Zoom In Zoom Out Reset image size

In the full 3D models without velocity dependence, (mesa20, mesa20_pert, mesa20_LR, and mesa20_LR_pert) we do not see any conclusive evidence for LESA. As discussed in Tamborra et al. (2014) and explored above, PNS convection is clearly critical for the LESA to develop as the neutrino lepton-number asymmetry stems from the asymmetry in the electron fraction distribution near and below the neutrinospheres as shown above. We do not attribute the presence of LESA in the mesa20_v_LR simulations to the explicit inclusion of the velocity dependence in the transport equations themselves, rather, we suggest that it is due to the effect of the improved neutrino transport on the PNS convection. With the included velocity dependence, neutrinos are able to advect with the flow. We see an earlier onset of PNS convection in the mesa20_v_LR simulation (at ∼100 ms after bounce) compared to, for example, the mesa20_LR simulation (which does not have equivalent PNS convection until ∼200 ms after bounce). Furthermore, at a given time, the PNS convection is stronger (and the region is wider) in the velocity-dependent simulation compared to the simulations without velocity dependence. Evidence for this conclusion is also provided by an examination of the 3D models from Couch & O'Connor (2014). In these simulations we see the presence of a strong asymmetry in the electron fraction distribution in the PNS core and the expected effect on the neutrino emission. Due to the low fidelity of the neutrino transport (i.e., neutrino leakage), we find this inconclusive evidence for the LESA. However, and important for this discussion, these models did show signs of strong PNS convection (even stronger than seen in mesa20_v_LR).

Using the quadrupole formula (Reisswig et al. 2011), we estimate the gravitational wave (GW) signal from our 2D and 3D simulations by computing the first time derivative of the reduced mass-quadrupole tensor via

$$\begin{eqnarray}&&\displaystyle \frac{{{dI}}_{{jk}}}{{dt}}=\int ({x}^{j}{v}^{k}+{x}^{k}{v}^{j}-2/3{\delta }^{{jk}}{x}^{i}{v}^{i}){dm},\end{eqnarray} \tag{ 15 }$$

and numerically taking the second time derivative during a post-processing step. In Figure 13 we show the plus polarization of the GW as seen from an observer on the equator (${h}_{+}^{\mathrm{eq}};$ this is the only nonzero component of the gravitational waves from axisymmetric simulations) for mesa20, mesa20_pert, mesa20_2D, and mesa20_2D_pert along with the GWs from our low-resolution simulations, mesa20_LR, mesa20_LR_pert, and our simulation with velocity dependence in the neutrino transport, mesa20_v_LR.

Zoom In Zoom Out Reset image size

The first observation is that due to the symmetries imposed in 2D, the gravitational wave signal is between 10 and 20 times the strength of the 3D signal. Note the scale difference, of a factor of 10, for the top two panels of Figure 13 compared to the bottom three panels. Gravitational waves are generated from large coherent matter motions in the turbulent gain region and near the PNS surface. In 2D, the axisymmetric assumption leads to an artificial enhancement of this coherent motion. This is in addition to the consequences of the reverse turbulent cascade in 2D, which also generates large-scale structures not seen in 3D simulations. This observation is in agreement with other analyses (Mueller & Janka 1997; Andresen et al. 2017), although other comparisons between the gravitational wave signature in 2D and 3D find much less of a suppression in 3D (Yakunin et al. 2017).

In the third and fourth panels of Figure 13 we compare the standard (third panel) and low (fourth panel) resolution simulations with and without the imposed perturbations. The most striking impact of the perturbations is the strong period of GW emission soon after 200 ms. At this time the perturbations at the top of the silicon shell are being accreted through the shock and amplifying the turbulent motions in the gain region. In addition to boosting the pressure support behind the shock and, at least temporarily, increasing the shock radius (see also Section 3.1.1) these increased turbulent motions give rise to strong GW emission. This effect is also visible in the 2D simulations that include perturbations (second panel versus first panel of Figure 13). To locate the source of this enhanced emission we compute spectrograms of the GW signal. We show these in Figure 14 for the mesa20 (left panel) simulation and mesa20_pert (right panel) simulation. The spectrograms generally behave similarly. There is no appreciable GW signal until ∼100 ms after bounce when the flow becomes turbulent. There are two distinct sources of GWs seen. First, there is a component that begins at ∼300 Hz and grows with time. This component has been associated with the contracting (but mass-gaining) PNS (Marek et al. 2009; Müller et al. 2013; Kuroda et al. 2017). The frequency we observe is slightly higher (by ∼20% at 450 ms) than expected based on the analytic predictions in Müller et al. (2013) because our gravitational wave predictions arise from Newtonian hydrodynamics and therefore are not either time dilated or redshifted. The other visible component is the low frequency (∼50–200 Hz) GW emission, which is arising from the convective and turbulent region behind the shock. The main difference in the GW spectrograms between the mesa20 and the mesa20_pert simulations is the excess emission near 200 ms. These spectrograms imply that the primary frequencies of this enhanced emission in mesa20_pert are close to ∼600 Hz, which at this time is the characteristic frequency of the PNS. After the perturbations accrete through the shock and stimulate the growth of turbulence, they then accrete onto the PNS, excite it, and lead to the production of GWs. Since we use only monopole gravity, we have been cautious in interpreting our GW signals. To ensure that the presence of this excess emission associated with the imposed presupernova perturbations is not a result of the monopole gravity assumption, we have also performed 2D simulations using ℓ_max = 16 in our monopole gravity solver. All of the gravity moments higher than ℓ = 0 are based on the Newtonian gravitational potential. We also see the excess gravitational wave emission due to progenitor perturbations in these simulations.

Zoom In Zoom Out Reset image size

In addition to the bursts of GWs around 200 ms due to the imposed progenitor perturbations, we also see bursts of GW emission associated with the collapse of spiral SASI waves. In particular, the increased emission at ∼360 ms in mesa20 and at ∼400 ms in mesa20_pert. While not shown, we see a similar feature in mesa20_LR at ∼370 ms. The emission is concentrated at the peak PNS frequency, but does contribute power to other frequencies as well. In terms of the dynamics, both the accretion of progenitor perturbations and the collapsing spiral SASI waves are very similar. They both excite turbulence in the gain region marked by an increase in the lateral kinetic energy (see Figure 1).

Finally, we compare mesa20_LR and mesa20_v_LR. During the first ∼300 ms we see broad agreement in the GW signal between these two simulations. Starting at ∼350 ms, coincident with the increased SASI activity, the power in the GW signal of mesa20_LR grows while mesa20_v_LR stays roughly constant. This is in agreement with our previous observations, which show reduced SASI motions in the mesa20_v_LR simulation due to the larger shock radius and slower recession of the PNS core (see Figure 4).