# General Relativistic Neutrino-Driven Turbulence in One-Dimensional Core-
Collapse Supernovae

Luca Boccioli Center for Astrophysics, Department of Physics, University of
Notre Dame, 225 Nieuwland Science Hall, Notre Dame, IN 46556, USA  Grant J.
Mathews Center for Astrophysics, Department of Physics, University of Notre
Dame, 225 Nieuwland Science Hall, Notre Dame, IN 46556, USA  Evan P. O’Connor
The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-106 91 Stockholm, Sweden

###### Abstract

Convection and turbulence in core-collapse supernovae (CCSNe) are inherently
three-dimensional in nature. However, 3D simulations of CCSNe are
computationally demanding. Thus, it is valuable to modify simulations in
spherical symmetry to incorporate 3D effects using some parametric model. In
this paper, we report on the formulation and implementation of general
relativistic (GR) neutrino-driven turbulent convection in the spherically
symmetric core-collapse supernova code GR1D. This is based upon the recently
proposed method of Supernova Turbulence in Reduced-dimensionality (STIR) in
Newtonian simulations from Couch et al. (2020). When the parameters of this
model are calibrated to 3D simulations, we find that our GR formulation of
STIR requires larger turbulent eddies to achieve a shock evolution similar to
the original STIR model. We also find that general relativity may alter the
correspondence between progenitor mass and successful vs. failed explosions.

Supernovae, Supernova dynamics, General Relativity, Gravitational Collapse,
Relativistic Astrophysics

††journal: ApJ

##  1 Introduction

Despite the fact that core-collapse supernovae (CCSNe) have been the focus of
computational simulations for more than 50 years, many fundamental details of
the mechanism driving the explosion of massive stars remain unknown.
Nevertheless, our knowledge has significantly improved since the first
theoretical (Colgate & White, 1966; Arnett, 1966) and observational (Sonneborn
et al., 1987) efforts. This has been aided by an exponential growth of
computing power.

In the early ’80s, one-dimensional (1D) spherically symmetric simulations
revealed the crucial role of neutrinos in triggering the explosion through the
so-called delayed neutrino-heating mechanism (Bethe & Wilson, 1985; Bruenn,
1985). A few years later, the first two-dimensional (2D) axi-symmetric
numerical investigations of CCSNe were also performed (Miller et al., 1993;
Herant et al., 1994), while more recently the extension to three dimensions
(3D) has finally become feasible (Fryer & Warren, 2002). Since then, multi-D
simulations have continued to improve, and have shed light on the details of
the explosion mechanism (cf. review in Ref. Müller (2016)). For example, the
importance of neutrino-driven turbulent convection is now widely recognized
(Radice et al., 2016, 2018; Mabanta & Murphy, 2018). Moreover, new phenomena
such as the Standing Accretion Shock Instability (SASI) (Blondin et al.,
2003), the Lepton-number Emission Self-sustained Asymmetry (LESA) (Tamborra et
al., 2014), Collective Neutrino Oscillations (CNO) (Duan et al., 2010) and
others, have been revealed in multi-D simulations. However, the impact of
these phenomena on the explosion is still a topic of active research.

After many years of theoretical and observational work, it is now well
established that spherically symmetric 1D simulations only lead to self-
consistent explosions of low mass (typically around 9-11 M⊙) progenitors.
These stars develop an oxygen-neon-magnesium core, whose collapse is triggered
by electron captures on Ne and Mg. Contrary to these so-called Electron
Capture Supernovae (ECSNe), the more common Fe-core CCSNe do not self-
consistently explode in spherically symmetric simulations that employ a modern
treatment of the nuclear equation of state (EOS) and neutrino interactions,
with the possible exception of a zero-metallicity, 9.6 M⊙ progenitor (Melson
et al., 2015; Betranhandy & O’Connor, 2020). Nevertheless, there have been a
number of simulations in 2D and 3D that have led to successful explosions
(Müller et al., 2012; Lentz et al., 2015; Janka et al., 2016; Bruenn et al.,
2016; O’Connor & Couch, 2018a; Müller et al., 2019; Burrows et al., 2020).
However, axial symmetry has recently been shown (Couch & Ott, 2015) to favor
explosions by artificially enhancing neutrino-heating behind the shock through
an inverse turbulence cascade that is not present in 3D. Therefore, it is
generally believed that the final explanation as to what causes the explosion
must come from simulations in three spatial dimensions.

Unfortunately, despite the technological improvements of the last few decades,
three-dimensional core-collapse supernova simulations still pose a difficult
computational challenge, even for modern supercomputers. In addition to this,
3D simulations from different groups only agree for a few tens of milliseconds
after bounce (Cabezón et al., 2018). Afterward, they begin to significantly
deviate from each other, even when they start from the same initial
conditions. The more accessible 2D simulations, that require a relatively
affordable CPU time, are also still complicated enough to the point that
different groups often find different outcomes for the explosion (see Table 1
from O’Connor & Couch (2018a)), although some promising benchmarking work has
already been done (Pan et al., 2019).

In comparison, modern 1D simulations are much faster to run and more
consistent across different codes (O’Connor et al., 2018). Hence, when
starting from the same initial conditions, most groups obtain similar results
(while older 1D simulations, with outdated sets of neutrino interactions, lead
to some inconsistencies [cf. Mathews et al. (2020)]). This guarantees a
somewhat solid foundation and makes 1D simulations an ideal bench-marking tool
for studying how different input physics can affect supernova explosions.

Precisely due to the low computational cost associated with 1D simulations,
there have been several attempts to artificially drive explosions in 1D codes
that would otherwise yield failed supernovae. The advantage of having models
capable of artificially driving explosions is quite straightforward: if the
models are reliable, one can efficiently explore the parameter space of
different physical phenomena. For example, one could explore the connection
between the explodability of CCSNe and the progenitor mass (O’Connor & Ott,
2011; Ugliano et al., 2012; Pejcha & Thompson, 2015; Ertl et al., 2016;
Sukhbold et al., 2016; Ebinger et al., 2020; Couch et al., 2020), or the
impact of different input physics, such as the nuclear EOS (Schneider et al.,
2019; Olson et al., 2016) and the weak interactions (Fischer et al., 2020; Guo
et al., 2020; Betranhandy & O’Connor, 2020)) on the explosion itself.
Furthermore, one can quantify the dependence of various nucleosynthesis
processes on the progenitor mass, or calculate the diffuse relic neutrino
background [c.f. review in Mathews et al. (2020)]. Moreover, 1D simulations
can efficiently study any aspect of supernovae that requires a large sample of
progenitors and/or evolution to late times post bounce.

A number of artificially driven explosion mechanisms have been proposed over
the past few decades, based upon different parametric models (Blinnikov &
Bartunov, 1993; Woosley & Weaver, 1995; Ugliano et al., 2012; Perego et al.,
2015; Sukhbold et al., 2016; Couch et al., 2020). For most of them, the
parameters were chosen by comparing the explosion energies and other relevant
observables to known measured properties of SN 1987a (Sonneborn et al., 1987)
and other SNe.

The advantage of comparing with observations is that one can directly test the
results with respect to real data. On the other hand, since these parametric
models are not derived from a physical mechanism, one cannot meaningfully
compare these results to 3D simulations and gain insight on the physics of the
explosion. Therefore, a different (and complementary) way of approaching the
problem is to develop a parametric model based upon a known physical mechanism
that is motivated by the success of 3D simulations (Müller et al., 2016; Couch
et al., 2020).

An obvious candidate for such a mechanism is neutrino-driven turbulent
convection, which plays an important role in CCSNe (Couch & Ott, 2015; Radice
et al., 2016, 2018; Mabanta & Murphy, 2018). Recently, Couch et al. (2020)
(hereafter CWO20) developed a model called STIR (Simulated Turbulence In
Reduced-dimensionality), based upon a modified Mixing Length Theory (MLT),
that simulates features of 3D turbulence in spherically symmetric simulations.

STIR is a Newtonian model for turbulent convection, and the simulations from
CWO20 only partially include general relativistic effects through a General
Relativistic Effective Potential (GREP) from Marek et al. (2006). However, we
know that General Relativity (GR) plays an important role in the explosion of
supernovae (Wilson & Mathews, 2003; Müller et al., 2012). Hence, simulations
in full GR are desirable.

In this paper we extend the STIR model to a general relativistic treatment. We
then analyze the differences between full GR and GREP. In Section 2 we
describe the essential physics of STIR and its generalization to GR. Results
and a comparison to 3D simulations are described in Section 3, followed by a
study of explodability vs. progenitor mass in Section 4. Conclusions are given
in Section 5. Throughout the manuscript, we chose to adopt natural units, i.e.
$G=c=M_{\odot}=1$.

##  2 Methods

###  2.1 The STIR model of Couch et al. (2020)

The importance of convection in CCSNe was recognized very early on, and the
first attempts at applying MLT to CCSN date back to Wilson & Mayle (1988) and
Bruenn et al. (1995). We follow a more modern version of those models, that
can be found in Mabanta & Murphy (2018), Mabanta et al. (2019) and CWO20. We
highlight that this formulation does consider both pressure and composition
gradients (see Appendix A of Müller et al. (2016) for an explicit proof of
this). However, it doesn’t capture lepton number-driven convection, relevant
at the surface of the PNS, where neutrinos are trapped below the
neutrinosphere and free streaming above it, creating a gradient which could
drive convection. However, given the complicated thermodynamics derivatives
involved in the modeling of this type of convection (Wilson & Mathews, 2003),
we leave its treatment to a future work.

As described in the aforementioned papers, the effects of turbulence can be
treated as a perturbation on a background fluid. This can be done by
decomposing a generic fluid variable $\psi$ into a mean component $\psi_{0}$
and a perturbed component $\psi^{\prime}$: $\psi=\psi_{0}+\psi^{\prime}$, with
$\langle\psi\rangle=\psi_{0}$, where $\langle...\rangle$ represents an average
over a proper spatial and time domain. Because of the chaotic nature of the
motions that characterize turbulence, the average should vanish (i.e.
$\langle\psi^{\prime}\rangle=0$). This decomposition, called Reynolds
averaging, can be applied to the compressible Euler equations to derive
conservation laws that contain explicit contributions from the turbulent
motion for total energy, momentum and mass. In the Newtonian limit, these are:

$\displaystyle\partial_{t}\langle\rho\rangle+\nabla\cdot(\rho_{0}\mathbf{v}_{0})=0~{}~{},$
(1)
$\displaystyle\begin{split}\partial_{t}\langle\rho\mathbf{v}\rangle+\nabla\cdot(\rho_{0}\mathbf{v}_{0}\times\mathbf{v}_{0}+P_{0}\mathbf{I})=&-\rho\mathbf{g}\\\
&-\nabla\cdot\langle\rho\mathbf{R}\rangle\end{split}~{}~{},$ (2)
$\displaystyle\begin{split}\partial_{t}\langle\rho\epsilon\rangle&+\nabla\cdot\mathbf{u_{0}}(\rho_{0}\epsilon_{0}+P_{0})\\\
=&-\rho_{0}\mathbf{v}_{0}\cdot\mathbf{g}-\nabla\cdot\mathbf{v}_{0}\langle\rho\mathbf{R}\rangle-\nabla\cdot\mathbf{F}_{e}+\rho_{0}\epsilon^{\prime}~{}~{}.\end{split}$
(3)

Here, $\rho$ is the mass density, $\mathbf{v}$ is the velocity, $P$ is the
pressure, $\mathbf{I}$ is the identity matrix, $\mathbf{g}$ is the
gravitational acceleration, $\epsilon$ is the internal energy per unit mass,
$\mathbf{R}$ is the Reynolds stress tensor and $\mathbf{F}_{e}$ is the energy
flux due to turbulence. In the above equations, we have neglected terms that
are proportional to the turbulent pressure. This is a good approximation in
regimes of low Mach numbers (for more details on the derivation of equations
(1)-(3) see Murphy & Meakin (2011), Mabanta & Murphy (2018), and CWO20).

In addition to Eqs. (1)-(3), one needs another equation to describe the
evolution of the specific turbulent kinetic energy
$K=\text{Tr}(\mathbf{R})/2$. The derivation of this equation can be found in
Murphy & Meakin (2011) and CWO20. Here, we simply highlight the important
steps that are relevant to our model.

For isotropic turbulence in spherical symmetry, the evolution equation for the
turbulent kinetic energy is:

$\begin{split}\frac{\partial\langle\rho K\rangle}{\partial
t}&+\frac{1}{r^{2}}\frac{\partial}{\partial r}[r^{2}(\langle\rho K\rangle
v_{r}+\langle\rho Kv_{r}^{\prime}\rangle)]\\\
=&-\langle\rho\mathbf{R}_{rr}\rangle\frac{\partial v_{r}}{\partial
r}+\langle\rho^{\prime}v_{r}^{\prime}\rangle
g-\rho\epsilon^{\prime}~{}~{}.\end{split}$ (4)

One can rewrite almost all of the unknown expressions in (4) in terms of
gradients and divergences of $v_{\text{turb}}\equiv v^{\prime}$ and
$\epsilon_{\text{turb}}\equiv\epsilon^{\prime}$, where:

$\epsilon_{\text{turb}}=\frac{v_{\text{turb}}^{3}}{\Lambda}$ (5)

with $\Lambda$ being the largest turbulent eddy size. The only quantities that
remain to be determined are the buoyancy force
$\langle\rho^{\prime}v_{r}^{\prime}\rangle g$ and $\Lambda$. Therefore, to
close the system of equations, one needs an expression that relates the
transport of turbulent energy to the typical speed of a buoyant convective
element. To do this, one refers to MLT, for which the buoyant force that makes
the convective element rise through the stratified fluid is:

$\langle\rho^{\prime}v^{\prime}\rangle g\approx\rho
v_{\text{turb}}\omega_{\textsc{\scriptsize bv}}^{2}\Lambda_{\text{mix}}$ (6)

with:

$\displaystyle\Lambda\equiv\Lambda_{\text{mix}}$
$\displaystyle=\alpha_{\textsc{\scriptsize mlt}}\frac{P}{\rho g},$ (7)
$\displaystyle\omega_{\textsc{\scriptsize bv}}^{2}$
$\displaystyle=g_{\text{eff}}\left(\frac{1}{\rho}\frac{\partial\rho}{\partial
r}-\frac{1}{\rho c_{s}^{2}}\frac{\partial P}{\partial r}\right).$ (8)

In the above equations, $\Lambda_{\text{mix}}$ is the mixing length,
$\omega_{\textsc{\scriptsize bv}}$ is the Brunt-Väisälä frequency, $c_{s}$ is
the sound speed, and $g_{\text{eff}}$ is the magnitude of the local effective
acceleration. For a fluid in hydrostatic equilibrium, $g_{\text{eff}}$ simply
reduces to the local gravitational acceleration $g$. More in general, however,
in the rest frame one should take the acceleration of the fluid into account.
Therefore, the total acceleration $g_{\text{eff}}$ can be expressed as:

$g_{\text{eff}}=g-v_{r}\frac{\partial v_{r}}{\partial r}~{},$ (9)

as described in CWO20, where $v_{r}$ is the radial velocity of the fluid.
During the pre-bounce and early post-bounce phases, both $v_{r}$ and $\partial
v_{r}/\partial r$ are small, so that $g_{\text{eff}}$ simply reduces to $g$.
Once the explosion sets in, however, the contribution from the acceleration of
the fluid is non-negligible, hence $g_{\text{eff}}$ becomes smaller than $g$.
This effectively shuts off turbulent convection in the exploding material by
driving $\omega_{\textsc{\scriptsize bv}}$ to zero.

The mixing length $\Lambda_{\text{mix}}$ is the average distance that a
convective element will travel before being mixed with (and increasing the
internal energy of) the surrounding material. The Brunt-Väisälä frequency
$\omega_{\textsc{\scriptsize bv}}$ is the rate at which the convective
elements are rising. As one can notice from Eq. (8),
$\omega_{\textsc{\scriptsize bv}}^{2}$ can be either positive or negative:
when $\omega_{\textsc{\scriptsize bv}}^{2}>0$ the fluid is convectively
unstable, i.e. convection is generated; when $\omega_{\textsc{\scriptsize
bv}}^{2}<0$ the fluid is convectively stable, i.e. convection is damped.
Ultimately, the main parameter of the model is $\alpha_{\textsc{\scriptsize
mlt}}$, which scales the mixing length to the pressure scale height in Eq.
(7). Typically $\alpha_{\textsc{\scriptsize mlt}}\sim O(1)$.

After closing the system of equations derived from the Reynolds decomposition,
one obtains the following evolution equation for the turbulent kinetic energy:

$\begin{split}\frac{\partial\rho v_{\text{turb}}^{2}}{\partial
t}&+\frac{1}{r^{2}}\frac{\partial}{\partial r}[r^{2}(\rho
v_{\text{turb}}^{2}v_{r}-\rho D_{K}\nabla v_{\text{turb}}^{2})]\\\ &=-\rho
v_{\text{turb}}^{2}\frac{\partial v_{r}}{\partial r}+\rho
v_{\text{turb}}\omega_{\textsc{\scriptsize
bv}}^{2}\Lambda_{\text{mix}}-\rho\frac{v_{\text{turb}}^{3}}{\Lambda_{\text{mix}}}~{}~{}.\end{split}$
(10)

In the above equation, $D_{K}$ is a diffusion coefficient due to turbulence,
defined as:

$D_{K}=\alpha_{K}v_{\text{turb}}\Lambda_{\text{mix}}~{}.$ (11)

Similar terms appear in the internal energy, electron fraction, and neutrino
energy density evolution equations (for the complete set of hydrodynamic
equations used in the model, see Eqs. 25-29 and 33 in CWO20). Therefore,
strictly speaking, STIR has 4 additional free parameters: $\alpha_{K}$,
$\alpha_{e}$, $\alpha_{\text{Ye}}$, $\alpha_{\nu}$. However, the convective
motions are not very sensitive to the value of these parameters, so we set
them to $1/6$ for simplicity, consistent with the choices of Müller et al.
(2016) and CWO20.

In the next section we describe a general relativistic version of the model
described above.

###  2.2 STIR in General relativity

The first attempts to create a general relativistic model for convection date
back to Thorne (1966). We will follow the same approach, but using a slightly
different formalism.

All the simulations described in this paper were run with the open-source,
spherically symmetric, general relativistic code GR1D (O’Connor & Ott, 2010).
The Boltzmann equation for neutrino transport is solved using an M1-scheme,
with opacity tables generated using the open-source code NuLib (O’Connor,
2015).

The metric evolved in GR1D is Schwarzchild-like:

$\begin{split}ds^{2}&=g_{\mu\nu}x^{\mu}x^{\nu}\\\
&=-\alpha(r,t)^{2}dt^{2}+X(r,t)^{2}dr^{2}+r^{2}d\Omega^{2},\end{split}$ (12)

where $\alpha$ and $X$ can be expressed as functions of a metric potential
$\phi$ (which reduces to the Newtonian potential in the Newtonian limit) and
the enclosed gravitational mass $M_{\text{grav}}$:

$\begin{split}\alpha(r,t)&=\exp[\phi(r,t)],\\\
X(r,t)&=\left(1-\frac{2M_{\text{grav}}(r,t)}{r}\right)^{-1/2}\end{split}$ (13)

For the present work, we first note that turbulence is mostly relevant far
from the proto-neutron star (PNS) where GR effects can be treated as a
perturbation. Therefore, one can simply make a few changes to the terms in Eq.
(10) without having to re-derive the entire Reynold’s decomposition. The
expression for $\omega_{\textsc{\scriptsize bv}}^{2}$, however, must be
carefully re-derived (a formal derivation of $\omega_{\textsc{\scriptsize
bv}}^{2}$ for a fluid in hydrostatic equilibrium can be found in Appendix A).
Far from the PNS, we invoke the following:

  1. (i)

replace the conserved variable $\rho$ with its GR counterpart, i.e.
$D=WX\rho$, where
$W=(1-v^{2})^{-1/2}$ and $v=Xv_{r}$;

  2. (ii)

multiply the RHS of Eq. (10) by $\alpha X$.

  3. (iii)

multiply the spatial flux in Eq. (10) by $\alpha/X$ (see O’Connor & Ott (2010)
for more details on the derivation of the hydrodynamic equations in GR1D.

We note, however, that GR effects are not small within the PNS, where interior
lepton number gradients can drive convection and affect the neutrino emission
(Wilson & Mathews, 2003). We leave a proper treatment of lepton number-driven
convection to a future work.

The expression of $\omega_{\textsc{\scriptsize bv}}^{2}$, with the assumption
of the background fluid being in hydrostatic equilibrium, is:

$\omega_{\textsc{\scriptsize bv}}^{2}=\frac{\alpha^{2}}{\rho
hX^{2}}\frac{\text{d}\phi}{\text{d}r}\left(\frac{\partial\rho(1+\epsilon)}{\partial
r}-\frac{1}{\rho c_{s}^{2}}\frac{\partial P}{\partial r}\right).$ (14)

For an accelerating fluid, the GR derivation becomes much less obvious, since
non-inertial frames are involved. Fortunately, GR effects are very small in
the outer layers, where the fluid has a non-negligible acceleration. Hence, we
can simply use the Newtonian expression (9) with ad hoc corrections for
relativistic space-time curvature, so that $\omega_{\textsc{\scriptsize
bv}}^{2}$ becomes:

$\omega_{\textsc{\scriptsize
bv}}^{2}=\frac{\alpha^{2}}{X^{2}}\left(\frac{\text{d}\phi}{\text{d}r}-v\frac{\partial
v}{\partial r}\right)\left(\frac{1}{\rho
h}\frac{\partial\rho(1+\epsilon)}{\partial r}-\frac{1}{\rho
hc_{s}^{2}}\frac{\partial P}{\partial r}\right),$ (15)

where now $v=Xv_{r}$.

The main difference between Eqs. (8) and (15) is the inclusion of
$\partial\rho\epsilon/\partial r$ in the latter. In the gain region the
internal energy decreases with radius, i.e. $\partial\rho\epsilon/\partial
r<0$. This decreases the magnitude of $\omega_{\textsc{\scriptsize bv}}$ and
therefore the amount of turbulence that is generated. We will come back to
this in Section 3.2

###  2.3 A word of caution: conservation of energy

As briefly noted in CWO20, the equations of the STIR model do not conserve
energy. As pointed out by Müller (2019), the first and second term on the RHS
of Eq. (4) break total energy conservation. Since turbulence is driven by
neutrinos interacting with the material in the gain region, the rate of
turbulent energy generation by buoyant driving is of the same order as the
volume-integrated neutrino heating rate $\dot{\mathcal{Q}}_{\nu}$. This is
also what (Mabanta et al., 2019) assume in their turbulence model in the gain
region, where the buoyant driving is given by $\beta\dot{\mathcal{Q}}_{\nu}$.

The extra energy coming from buoyant driving breaks the energy-conserving
structure of the hydrodynamical equations. As pointed out in Mabanta & Murphy
(2018) this can be thought of as energy extracted from the free energy stored
in unstable compositional and thermodynamic gradients. Therefore, although the
energy is coming from a physically motivated phenomenon, it cannot be
consistently modeled in spherical symmetry, as recognized by Müller (2019).

If one looks at the problem of energy conservation in other artificially-
driven explosions in one dimension, it is clear that energy is not conserved
in any of them. This reflects the main issue with 1D models of CCSNe: they do
not self-consistently explode. Hence, generally speaking, to drive an
explosion in one dimension one needs more efficient neutrino heating. This
extra-heating is injected in the model following different prescriptions, that
can be based on physical phenomena (Mabanta & Murphy, 2018; Couch et al.,
2020) and/or observational constraints (Ugliano et al., 2012; Sukhbold et al.,
2016), or simply by well-reasoned tweaking of some parameters (O’Connor & Ott,
2011; Perego et al., 2015). Until multi-dimensional models become
computationally more affordable, the above methods can give different and
complementary insights into the physics behind the explosion mechanism of
CCSNe.

##  3 Results: Comparison with 3D simulations

###  3.1 Results using an effective potential

Following the approach of CWO20, we compare our GREP model to the mesa20_LR_v
3D simulation by O’Connor & Couch (2018b). In the original STIR paper, the
advantage of carrying out such a comparison is that both O’Connor & Couch
(2018b) and CWO20 use the code FLASH, so they can directly compare their 1D
and 3D simulations with everything else unchanged. Despite using a different
code, we still expect our GREP model to show very little differences to the
original STIR, since GR1D and FLASH have been shown to yield similar results
in 1D (O’Connor et al., 2018).

In all the simulations presented in this section we used the same setup chosen
by CWO20. That is, we simulated the collapse of a 20 M⊙ progenitor from Farmer
et al. (2016), adopting the SFHo EOS (Steiner et al., 2013) and assuming
Nuclear Statistical Equilibrium (NSE) everywhere. The neutrino radiation
transport solver used in FLASH closely resembles the one used in GR1D
(O’Connor, 2015), and we used the same set of NuLib opacities adopted in
CWO20. For the neutrinos, we used 12 energy groups geometrically spaced up to
$\sim$250 MeV.

In The upper panels of Figure 1 we plot the shock radius versus time and the
turbulent velocity profile at $\sim$135 ms post-bounce for our GREP model (to
be compared to Figures 1 and 2 of CWO20). GR1D consistently yields larger
values than FLASH for the turbulent velocity at a given
$\alpha_{\textsc{\scriptsize mlt}}$. This translates into larger values of the
shock radius at a given time. Except for these small differences, the two
models agree very well, yielding explosions for $\alpha_{\textsc{\scriptsize
mlt}}>1.2$. It is worth noting, however, that the 3D results show features
that are not present in our MLT-like model. Specifically, the convective speed
in 3D has a slower-decaying tail at 50-80 km. As pointed out by CWO20, this is
due to angular variations in the 3D model, rather than presence of convection
in the region below the gain layer. As one would expect, in fact, convection
in our model shuts off at 80 km, which is approximately the location of the
gain layer. More interesting is the lack of PNS convection at 25 km, not
captured by our model. A possible explanation for this is that we are not
taking lepton number-driven convection into account, and with a more careful
treatment of this type of convection, this feature might be correctly
reproduced, but this goes beyond the scope of this paper.

Figure 1: The plots on the upper row were generated using our GREP model,
while the plots on the bottom row were generated using full GR. Panels (a,c)
show the time evolution of the shock radius for different values of the
parameter $\alpha_{\textsc{\scriptsize mlt}}$, and can be compared to Figure 2
from Couch et al. (2020). Panels (b,d) show a snapshot at $\sim 135$ ms post
bounce of $v_{\text{turb}}$, and can be compared to Figure 1 from Couch et al.
(2020). The dashed lines are from the 3D simulation of O’Connor & Couch
(2018b).

###  3.2 Results using GR

We approached the analysis of the simulations in full GR from two different
points of view: on the one hand, we can compare the GR and GREP models with
everything else being unchanged; on the other hand, we can compare the GR
model to a 3D simulation in full GR.

For the first part of the analysis, we show the results of our GREP model in
the upper panels of Figure 1, while the bottom panels refer to the runs with
full GR. The first and most important difference to notice is the range of
$\alpha_{\textsc{\scriptsize mlt}}$ used in the GREP and GR simulations. To
achieve turbulent velocities and shock radii that are similar to the GREP
results, the model in full GR requires an $\alpha_{\textsc{\scriptsize mlt}}$
that is $\sim 20\%$ larger. To understand the origin of this increase in
$\alpha_{\textsc{\scriptsize mlt}}$, it’s useful to look at the form of
$\omega_{\textsc{\scriptsize bv}}$. As already pointed out in Section 2.2, the
inclusion of the internal energy gradient in eq. (15), compared to eq. (8), is
in fact the main difference between the GREP and GR models.

In the gain region, which is where turbulence is most relevant, the internal
energy gradient is negative, and therefore it decreases
$\omega_{\textsc{\scriptsize bv}}^{2}$, making the fluid more stable against
convection. The inclusion of $\rho\epsilon$ in the definition of
$\omega_{\textsc{\scriptsize bv}}$ is therefore needed for a realistic
characterization of turbulent convection. Indeed, if one takes the form of
$\omega_{\textsc{\scriptsize bv}}$ from GR and implements it in the GREP
model, the value of $\alpha_{\textsc{\scriptsize mlt}}$ needed to achieve an
explosion increases, becoming comparable to the one we use for the GR model.

The second part of the analysis was carried out by comparing the full GR model
to the 3D octant simulations of Schneider et al. (2019). For this part of the
analysis, we simulated a 20 M⊙ progenitor from Woosley & Heger (2007) with two
of the EOS’s used in Schneider et al. (2019): SLy41.0 (Chabanat et al., 1997)
and LS220 (a version of the original EOS with incompressibility
$K_{\infty}=220$ MeV from Lattimer & Swesty (1991) recomputed using the SRO
code from da Silva Schneider et al. (2017)). For the neutrino opacities, we
adopted the same NuLib tables used by Schneider et al. (2019), with 24 energy
groups logarithmically spaced up to 280 MeV.

Like Schneider et al. (2019), at 20 ms after bounce we switch from inelastic
electron-neutrino scattering to a simpler elastic treatment, and we also shut
off velocity-dependent terms in the transport equation. Unlike them, we
consistently keep 24 energy groups throughout the whole simulation, while they
switch to 16 energy groups at 20 ms post-bounce. Our choice ensures a smooth
transition in our code from inelastic to elastic electron-neutrino scattering.

A comparison between the 3D data and our 1D model is presented in Figure 2. As
one can immediately see, for the LS220 case, both the shock radius and the
turbulent velocity profile show very good agreement between the 3D and our 1D
results. As already pointed out in section 3.1, minor differences in the
turbulent velocity profile at small radii can be explained as the consequence
of angular variations in the 3D simulation. The SLy4 case is a bit more
complicated, and the agreements is not as good as in the LS220 case. After
$\sim$300 ms the shock trajectory in 1D starts deviating from the 3D case,
suggesting that turbulence in our MLT-like model might be slightly more
efficient than in 3D, hence generating a successful explosion. A similar
feature can be seen in the LS220 case, in which after $\sim$300ms the 3D shock
trajectory drops below the one in 1D. Finally, it’s worth mentioning that
since the 3D simulations reported in Figure 2 were run imposing octant
symmetry, SASI cannot develop. This represents a further confirmation that our
GR model accurately predicts the development of turbulence in the pre-
explosion phase of CCSNe.

It is also worth noting that, for larger values of
$\alpha_{\textsc{\scriptsize mlt}}$, the turbulence inside the PNS is weaker.
This is consistent with the findings of CWO20. In their paper, they (semi-
quantitatively) assert that the convection inside the PNS is too fast to allow
a build-up of turbulence energy. This is because the gradients are quickly
flattened so there isn’t enough time for turbulence to develop. This is the
reason why for small values of $\alpha_{\textsc{\scriptsize mlt}}$ our model
is able to generate some turbulence outside the PNS, but as convection becomes
more efficient (i.e. for larger values of $\alpha_{\textsc{\scriptsize mlt}}$)
the gradients are flattened and $\omega_{\textsc{\scriptsize bv}}^{2}$ becomes
smaller, resulting in less turbulent energy generated above the PNS.

Figure 2: Comparison of our GR turbulent convection model (solid lines) as a
function of $\alpha_{\textsc{\scriptsize mlt}}$ with the 3D octant simulations
from Schneider et al. (2019) (dashed lines). The upper (lower) panels refer to
the simulations using the LS220 (SLy4) EOS. The panels on the left show the
shock evolution, while the panels on the right show the turbulent velocity
profile at $\sim 220$ ms after bounce.

###  3.3 Another word of caution: Resolution

It is well known that, in 3D, low resolution simulations can overestimate the
amount of turbulence produced in the gain region. The turbulence cascade is in
fact artificially truncated at larger scales than what would otherwise happen
in higher resolution simulations (Radice et al., 2016). Low resolution leads
to an accumulation of turbulent energy at the largest scales, resulting in an
additional pressure support in the gain region, which facilitates an
explosion. In an MLT approach, this can typically be avoided, but other
problems may arise (for a more detailed discussion see Appendix B).

For the remainder of the paper, we chose to perform simulations at a higher
spatial resolution and with more energy groups with respect to the results
presented in Section 3. As explained in Appendix B, a higher resolution in
space and in neutrino energy groups can change how much turbulent energy is
injected in the model. However, it is quite unlikely that this will affect the
qualitative results (i.e. explosion or fallback) of our simulation, because we
are using a parametric model where one can choose how much turbulent energy
will be injected in the model simply by changing the value of
$\alpha_{\textsc{\scriptsize mlt}}$. Therefore, if a simulation run with
higher spatial resolution develops more (or less) turbulence, one can obtain
the same qualitative result in a lower resolution simulation by increasing (or
decreasing) the value of $\alpha_{\textsc{\scriptsize mlt}}$.

For the remainder of the paper, however, we chose to perform all of the
simulations with a higher spatial resolution and with more neutrino energy
groups with respect to Section 3.

##  4 Progenitor Study

After validating our turbulent convection model by comparing it to the 3D
results of Couch & O’Connor (2014) and Schneider et al. (2019), we simulated
the collapse and subsequent shock revival of 20 progenitors from Sukhbold et
al. (2016) for different values of $\alpha_{\textsc{\scriptsize mlt}}$, chosen
such that the fraction of successful explosions is approximately between 25%
and 80%. All of the simulations were run using a high spatial resolution,
shown in the upper panel of Figure 4 as a dashed line, with 18 neutrino energy
groups geometrically spaced from 0 to 280 MeV.

From the analysis carried out in section 3, it is clear that once an explosion
sets in for a certain value $\alpha_{\textsc{\scriptsize mlt}}^{*}$, that
progenitor will explode for all $\alpha_{\textsc{\scriptsize
mlt}}>\alpha_{\textsc{\scriptsize mlt}}^{*}$ (see Figures 1 and 2). Therefore,
if a given progenitor explodes for a given value of
$\alpha_{\textsc{\scriptsize mlt}}$, it will also explode for higher values of
$\alpha_{\textsc{\scriptsize mlt}}$. Analogously, if a given progenitor ends
up in a failed explosion for a given value of $\alpha_{\textsc{\scriptsize
mlt}}$, it will also fail for lower values of $\alpha_{\textsc{\scriptsize
mlt}}$. Hence, we did not run simulations for those values of
$\alpha_{\textsc{\scriptsize mlt}}$ that could be already predicted by the
outcome of the simulations at other values of $\alpha_{\textsc{\scriptsize
mlt}}$.

The results of our GREP model agree with the ones obtained by CWO20, although
shifted by $\Delta\alpha_{\textsc{\scriptsize mlt}}\simeq 0.05$, as can be
seen by comparing the upper panel of Figure 3 with Figure 6 from their paper.
The shift is due to the fact that we use a finer resolution in space and
energy, and also to the fact that the numerical methods used to solve the
hydrodynamics and the neutrino transport are slightly different.

More importantly, however, Figure 3 shows that the GR treatment of turbulent
convection does not simply reproduce the STIR results for a higher value of
$\alpha_{\textsc{\scriptsize mlt}}$, as Figure 1 might suggest, but it
modifies the explosion pattern of CCSNe. To accurately characterize the
differences between the explodability patterns of the GR and GREP models, a
study with more progenitors and a finer resolution in
$\alpha_{\textsc{\scriptsize mlt}}$ would be desired. That, however, is beyond
the scope of the present work.

We can conclude from Figure 3 that general relativity changes how likelihood
for certain progenitors to explode. Focusing on the
$\alpha_{\textsc{\scriptsize mlt}}=1.27$ and $\alpha_{\textsc{\scriptsize
mlt}}=1.48$ panels, we find that the GREP model gives successful explosions
for the 24, 25 and 30 M⊙ progenitors, and failed explosions for the 18 M⊙. In
the GR model, the explodability is the exact opposite.

The explodability pattern of the GR model with $\alpha_{\textsc{\scriptsize
mlt}}=1.48$ is intermediate between the results of CWO20 and Sukhbold et al.
(2016): (i) like the former (and unlike the latter), we find that low mass
progenitors with M = 13 -15 M⊙ result in failed explosions; (ii) like the
latter (and unlike the former) we find that higher mass progenitors with M =
24 -25 M⊙ result in failed explosions.

It is worth mentioning that the explodability pattern found in the GR model
with $\alpha_{\textsc{\scriptsize mlt}}=1.48$ cannot be reproduced using the
GREP model. Focusing on the upper panel of Figure 3, one can notice how, in
the $\alpha_{\textsc{\scriptsize mlt}}=1.25$ case, the 25 M⊙ progenitor fails,
but at the same time all progenitors below 22 M⊙ fail as well. This shows that
the GREP model cannot produce successful explosions for low mass progenitors
and at the same time failed explosions of the 24-25 M⊙ progenitors. Notably,
the GR model with $\alpha_{\textsc{\scriptsize mlt}}=1.5$ reproduces the same
explodability pattern found in the GREP model with
$\alpha_{\textsc{\scriptsize mlt}}=1.27$ (with the only exception of the 18 M⊙
progenitor). This tells us that: (i) the threshold between successful and
failed explosions is a very steep function of $\alpha_{\textsc{\scriptsize
mlt}}$; and (ii) GR can indeed reproduce the GREP results. However, GR also
shows another possible pattern not yet observed in the GREP model of CWO20.

Multi-D simulations cannot yet produce results analogous to Figure 3, but the
number of 3D simulations in the last 5 years has dramatically increased. In
particular, a (relatively) large set of progenitors has been recently
simulated by Burrows et al. (2020). In their work, they simulate 14
progenitors from Sukhbold et al. (2016), 12 of which are in the range $M=9-20$
M⊙, and the remaining two have $M=25$ M⊙ and $M=60$ M⊙. Moreover, they use the
same SFHo EOS that we employed for our progenitor study. However, like the
vast majority of 3D simulations, they only include general relativistic
effects in the hydrodynamics through the GREP model of Marek et al. (2006),
which they also use to account for redshift effects in the neutrino transport.
Hence, we can qualitatively compare the outcome of our simulations with
theirs.

In Burrows et al. (2020) all of the progenitors explode except for the 13, 14
and 15 M⊙, which is exactly what we find in the GREP model with
$\alpha_{\textsc{\scriptsize mlt}}=1.27$. The GR model with
$\alpha_{\textsc{\scriptsize mlt}}=1.48$ finds that the 25 M⊙ should also
fail, and it is interesting to note that in Burrows et al. (2020) that
progenitor is the one with the explosion occurring at the latest time, i.e. it
is triggered at $\sim 400$ ms post bounce. Therefore, it can be considered a
weak explosion, and a slightly less effective neutrino heating would cause the
supernova to fail. This is consistent with what the GR model is showing. That
is, among all progenitors in the range 9-25 M⊙, the 25 M⊙ progenitor is the
one that requires the largest value of $\alpha_{\textsc{\scriptsize mlt}}$ to
explode (i.e. $\sim 1.5$).

This shows that including general relativistic effects can change how
convective turbulence behaves in this MLT model. We cannot predict if this
effect translates to multi-dimensional simulations, but these results are
suggesting that GR (rather than GREP) might have a significant impact on the
explosion of CCSNe. Only a detailed comparison, done with 2D and 3D codes,
between full-GR and Newtonian simulations, and across multiple progenitors,
can clarify whether this effect translates to higher dimensions, where
turbulent convection is generated self-consistently.

Figure 3: Explosion pattern of CCSN for the GREP (upper panel) and GR (lower
panel) models as a function of the Zero Age Main Sequence mass. Green bands
represent successful explosions (i.e. the shock has reached 500 km), while
black bands represent failed explosions.

##  5 Conclusions

In this article, we have shown how to extend STIR, the MLT-like model of
CWO20, to a full general relativistic formalism. After verifying that we can
reproduce the results of CWO20 using the same GREP model that they developed,
we showed how a GR version of this model needs a larger value of
$\alpha_{\textsc{\scriptsize mlt}}$ to reproduce similar shock dynamics. The
reason behind this is the inclusion of the internal energy gradient in eq.
(15), which reduces the magnitude of $\omega_{\textsc{\scriptsize bv}}$, hence
inhibiting the generation of turbulence, so that larger values of
$\Lambda_{\rm mix}$ are needed to develop convective mixing as strong as in a
Newtonian formalism.

The GR model reproduces very well the shock dynamics and turbulent energy
generation in the octant 3D simulations of Schneider et al. (2019), performed
with two different equations of state.

After comparing our results to 3D results, we ran simulations of 20 different
progenitors from (Sukhbold et al., 2016) for different values of
$\alpha_{\textsc{\scriptsize mlt}}$, for both the GREP and GR models. We find
that GR changes the explodability pattern of CCSNe. Specifically, the 24-25 M⊙
progenitors seem to be hared to explode with the GR model. This results in an
explodability as a function of progenitor mass that is intermediate between
the results of CWO20 and (Sukhbold et al., 2016). However, the GR model also
shows, for $\alpha_{\textsc{\scriptsize mlt}}=1.5$, an explodability that is
compatible with the results of CWO20.

The authors would like to thank Sean Couch, Andre da Silva Schneider and
Mackenzie Warren for fruitful discussions. Work at the University of Notre
Dame supported by the U.S. Department of Energy under Nuclear Theory Grant DE-
FG02-95-ER40934. EOC would like to acknowledge Vetenskapsrådet (the Swedish
Research Council) for supporting this work under award numbers 2018-04575 and
2020-00452.

##  Appendix A Derivation of the Brunt-Väisälä frequency

A derivation of $\omega_{\textsc{\scriptsize bv}}$ in general relativistic
hydrostatic equilibrium can be found in Appendix C of Müller et al. (2013).
For clarity, here we go over a similar derivation, but for a different gauge,
trying not to skip any non-obvious steps. In the fluid reference frame, the
velocity is simply the convective velocity $v_{c}$. Hence, the 4-velocity is:

$u^{\mu}=(W/\alpha,Wv_{c},0,0)$ (A1)

Conservation of momentum requires that $\nabla_{\mu}T^{\mu r}=0$, where:

$\displaystyle T^{\mu\nu}$ $\displaystyle=(\rho h+\delta\rho
h)u^{\mu}u^{\nu}+g^{\mu\nu}P$ (A2) $\displaystyle T^{tt}$
$\displaystyle=\frac{(\rho h+\delta\rho h)W^{2}-P}{\alpha^{2}},\qquad
T^{tr}=\frac{(\rho h+\delta\rho h)W^{2}v_{c}}{\alpha X},\qquad
T^{rr}=\frac{(\rho h+\delta\rho h)W^{2}v_{c}^{2}+P}{X^{2}}$ (A3)
$\displaystyle T^{\theta\theta}$ $\displaystyle=\frac{P}{r^{2}},\qquad
T^{\phi\phi}=\frac{P}{r^{2}\sin^{2}\theta}$ (A4)

Note that the convective eddy is in pressure equilibrium with the
surroundings. Therefore, there is a single pressure term in Eq. (A2). Writing
out the conservation of momentum one gets:

$\begin{split}\nabla_{\mu}T^{\mu
r}&=0=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}T^{\mu
r})+\Gamma^{r}_{\rho\sigma}T^{\rho\sigma}\\\
&=\frac{1}{\sqrt{-g}}\partial_{t}(\sqrt{-g}T^{tr})+\frac{1}{\sqrt{-g}}\partial_{r}(\sqrt{-g}T^{rr})+\Gamma^{r}_{tt}T^{tt}+\Gamma^{r}_{rr}T^{rr}+\Gamma^{r}_{\theta\theta}T^{\theta\theta}+\Gamma^{r}_{\phi\phi}T^{\phi\phi}+2\Gamma^{r}_{rt}T^{rt}\\\
&=\frac{1}{\alpha X}\frac{\partial}{\partial t}\left[(\rho h+\delta\rho
h)W^{2}v_{c}\right]+\frac{1}{\alpha Xr^{2}}\frac{\partial}{\partial
r}\left(\frac{\alpha r^{2}}{X}\left[P+(\rho h+\delta\rho
h)W^{2}v_{c}^{2}\right]\right)\\\ &+\frac{1}{X^{2}}[(\rho h+\delta\rho
h)W^{2}-P]\partial_{r}\phi+\frac{P}{X^{3}}\partial_{r}X-\frac{r}{X^{2}}\frac{P}{r^{2}}-\frac{r\sin^{2}\theta}{X^{2}}\frac{P}{r^{2}\sin^{2}\theta}+\frac{2(\rho
h+\delta\rho h)W^{2}v_{c}}{\alpha X^{2}}\partial_{t}X\end{split}$ (A5)

The Christoffel symbols for the metric (12) can be found in Table A.1 from
O’Connor & Ott (2010).

By explicitly evaluating the derivative in the second term of Eq. (A5), it can
be easily shown that:

$\frac{1}{\alpha Xr^{2}}\frac{\partial}{\partial r}\left(\frac{\alpha
r^{2}}{X}P\right)=\frac{1}{X^{2}}\frac{\partial P}{\partial
r}+\frac{P}{X^{2}}\frac{\partial\phi}{\partial
r}+\frac{2P}{rX}-\frac{P}{X^{3}}\partial_{r}X$ (A6)

In addition to this, using Eqs. A.5-7 from O’Connor & Ott (2010) one has:

$v_{c}\partial_{t}X=v_{c}\frac{X^{3}}{r}\partial_{t}m\propto\cancelto{0}{v_{c}^{2}}$
(A7)

By combining Eqs. (A5) and (A6), and using $\rho h_{t}=\rho h+\delta\rho h$ we
arrive at:

$\frac{1}{X^{2}}\frac{\partial P}{\partial r}+\frac{1}{\alpha
X}\frac{\partial}{\partial t}[\rho h_{t}W^{2}v_{c}]+\frac{1}{\alpha
Xr^{2}}\frac{\partial}{\partial r}\left(\frac{\alpha r^{2}}{X}\left[\rho
h_{t}W^{2}\cancelto{0}{v_{c}^{2}}\right]\right)=-\frac{\rho
h_{t}W^{2}}{X^{2}}\frac{\partial\phi}{\partial r}~{}~{}.$ (A8)

Here we have ignored terms of order $v_{c}^{2}$. Now one can use the condition
of hydrostatic equilibrium:

$\frac{\partial P}{\partial r}=-\rho h\frac{\partial\phi}{\partial r}.$ (A9)

By considering subsonic motions in the gain region (hence
$v_{c}<<c_{s}<<c\Rightarrow W^{2}=1$), one obtains:

$\frac{1}{\alpha X}\frac{\partial}{\partial t}\left(\rho
h_{t}v_{c}\right)=-\frac{\delta\rho h}{X^{2}}\frac{\partial\phi}{\partial
r}=-\frac{1}{X^{2}}\frac{\partial\phi}{\partial r}C_{L}\delta r~{}~{},$ (A10)

where for $C_{L}$ we have used the Schwarzschild discriminant from Thorne
(1966):

$C_{L}=\left(\frac{\partial\rho(1+\epsilon)}{\partial
r}-\frac{1}{c_{s}^{2}}\frac{\partial P}{\partial r}\right)~{}~{}.$ (A11)

By using conservation of energy and rest mass (i.e. equations A.13 and A.17
from O’Connor & Ott (2010)), and assuming pressure equilibrium, it can be
shown that:

$v_{c}\frac{\partial\rho h_{t}}{\partial t}\propto\cancelto{0}{v_{c}^{2}}$
(A12)

Therefore Eq. (A10) becomes:

$\frac{\partial v_{c}^{r}}{\partial t}=-\frac{\alpha}{\rho
h_{t}X^{2}}\frac{\partial\phi}{\partial r}C_{L}\delta r~{}~{}.$ (A13)

Since $\ddot{r}=\alpha\dfrac{\partial v_{c}^{r}}{\partial t}$ (where the
factor of $\alpha$ converts $t$ to proper time), we finally have:

$\ddot{r}=-\frac{\alpha^{2}}{X^{2}\rho h_{t}}\frac{\partial\phi}{\partial
r}C_{L}\delta
r\;\Longrightarrow\;\omega_{\text{BV}}^{2}=\frac{\alpha^{2}}{\rho
h_{t}X^{2}}\frac{\partial\phi}{\partial
r}\left(\frac{\partial\rho(1+\epsilon)}{\partial
r}-\frac{1}{c_{s}^{2}}\frac{\partial P}{\partial r}\right)$ (A14)

Figure 4: Evolution of the shock radius for the mesa20_LR_v progenitor by
O’Connor & Couch (2018b) for the GR (blue lines) and GREP models (red lines).
The GR models were all run for $\alpha_{\textsc{\scriptsize mlt}}=1.5$, while
the GREP models were all run with $\alpha_{\textsc{\scriptsize mlt}}=1.3$. (a)
The upper panel compares low (solid lines) and high (dashed lines) spatial
resolution simulations. All cases were run with 18 energy groups for
neutrinos. (b) The lower panel compares simulations that use 18 (dashed lines)
and 12 (solid lines) energy groups for neutrinos. These were run using the low
resolution spatial grid shown as a solid line in panel (a)

##  Appendix B Dependence of MLT-like turbulence on spatial and energy
resolution

Because turbulence is most effective behind the shock, it is crucial to have
high resolution in the proximity of the shock, especially in the first few
hundred milliseconds after bounce. In fact, in this time window, the shock is
stalling at $\sim 150$ km, and turbulent energy can provide the extra-pressure
necessary to trigger the explosion. In our numerical treatment of equation
(10) we use a 5-point interpolation formula to evaluate the gradient of
$v_{r}$, and to avoid numerical instabilities we shut off turbulence within
two zones from the shock. This ensures that the turbulence is not advected
through the shock, which, being a discontinuity in pressure and radial
velocity, would make the evaluation of gradients problematic. However, this
means that in our model there is no turbulence in the two zones right behind
the shock. Therefore, if the resolution behind the shock is poor (typically
$\Delta r\sim 1/1.5$ km at $r\sim 100/150$ km), one ends up shutting off
turbulence in the most important regions (i.e. close to the shock).

We ran some tests by changing the spatial resolution for the GREP and GR model
(see upper panel of Fig. 4) for the same mesa20_LR_v progenitor used in
Section 3. From these tests it appears that, with the increased resolution in
the GR model, the shock is pushed to larger radii at early times. This allows
a widening of the gain region, generating a stronger explosion. This feature
isn’t present in the GREP model, which on the contrary shows weaker explosions
at higher spatial resolutions. We also ran some tests by changing the
resolution in energy. For this case the GR and GREP models showed the same
behaviour, i.e. higher energy resolution yields weaker explosions (see lower
panel of Fig. 4).

The dependence of the shock dynamics on resolution is not as straightforward
as what might be expected based on simple arguments. However, the effects are
relatively small, limited to a few percents in $\alpha_{\textsc{\scriptsize
mlt}}$, and therefore do not significantly affect our analysis.

## References

  * Arnett (1966) Arnett, W. D. 1966, Canadian Journal of Physics, 44, 2553, doi: 10.1139/p66-210
  * Bethe & Wilson (1985) Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14, doi: 10.1086/163343
  * Betranhandy & O’Connor (2020) Betranhandy, A., & O’Connor, E. 2020, Phys. Rev. D, 102, 123015, doi: 10.1103/PhysRevD.102.123015
  * Blinnikov & Bartunov (1993) Blinnikov, S. I., & Bartunov, O. S. 1993, A&A, 273, 106
  * Blondin et al. (2003) Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003, ApJ, 584, 971, doi: 10.1086/345812
  * Bruenn (1985) Bruenn, S. W. 1985, ApJS, 58, 771, doi: 10.1086/191056
  * Bruenn et al. (1995) Bruenn, S. W., Mezzacappa, A., & Dineva, T. 1995, Phys. Rep., 256, 69, doi: 10.1016/0370-1573(94)00102-9
  * Bruenn et al. (2016) Bruenn, S. W., Lentz, E. J., Hix, W. R., et al. 2016, ApJ, 818, 123, doi: 10.3847/0004-637X/818/2/123
  * Burrows et al. (2020) Burrows, A., Radice, D., Vartanyan, D., et al. 2020, MNRAS, 491, 2715, doi: 10.1093/mnras/stz3223
  * Cabezón et al. (2018) Cabezón, R. M., Pan, K.-C., Liebendörfer, M., et al. 2018, A&A, 619, A118, doi: 10.1051/0004-6361/201833705
  * Chabanat et al. (1997) Chabanat, E., Bonche, P., Haensel, P., Meyer, J., & Schaeffer, R. 1997, Nucl. Phys. A, 627, 710, doi: 10.1016/S0375-9474(97)00596-4
  * Colgate & White (1966) Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626, doi: 10.1086/148549
  * Couch & O’Connor (2014) Couch, S. M., & O’Connor, E. P. 2014, ApJ, 785, 123, doi: 10.1088/0004-637X/785/2/123
  * Couch & Ott (2015) Couch, S. M., & Ott, C. D. 2015, ApJ, 799, 5, doi: 10.1088/0004-637X/799/1/5
  * Couch et al. (2020) Couch, S. M., Warren, M. L., & O’Connor, E. P. 2020, The Astrophysical Journal, 890, 127, doi: 10.3847/1538-4357/ab609e
  * da Silva Schneider et al. (2017) da Silva Schneider, A., Roberts, L., & Ott, C. 2017, Phys. Rev. C, 96, doi: 10.1103/PhysRevC.96.065802
  * Duan et al. (2010) Duan, H., Fuller, G. M., & Qian, Y.-Z. 2010, Annual Review of Nuclear and Particle Science, 60, 569, doi: 10.1146/annurev.nucl.012809.104524
  * Ebinger et al. (2020) Ebinger, K., Curtis, S., Ghosh, S., et al. 2020, ApJ, 888, 91, doi: 10.3847/1538-4357/ab5dcb
  * Ertl et al. (2016) Ertl, T., Janka, H. T., Woosley, S. E., Sukhbold, T., & Ugliano, M. 2016, ApJ, 818, 124, doi: 10.3847/0004-637X/818/2/124
  * Farmer et al. (2016) Farmer, R., Fields, C. E., Petermann, I., et al. 2016, The Astrophysical Journal Supplement Series, 227, 22, doi: 10.3847/1538-4365/227/2/22
  * Fischer et al. (2020) Fischer, T., Typel, S., Röpke, G., Bastian, N.-U. F., & Martínez-Pinedo, G. 2020, Phys. Rev. C, 102, 055807, doi: 10.1103/PhysRevC.102.055807
  * Fryer & Warren (2002) Fryer, C. L., & Warren, M. S. 2002, ApJ, 574, L65, doi: 10.1086/342258
  * Guo et al. (2020) Guo, G., Martínez-Pinedo, G., Lohs, A., & Fischer, T. 2020, Phys. Rev. D, 102, 023037, doi: 10.1103/PhysRevD.102.023037
  * Herant et al. (1994) Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, ApJ, 435, 339, doi: 10.1086/174817
  * Janka et al. (2016) Janka, H.-T., Melson, T., & Summa, A. 2016, Annual Review of Nuclear and Particle Science, 66, 341, doi: 10.1146/annurev-nucl-102115-044747
  * Lattimer & Swesty (1991) Lattimer, J. M., & Swesty, D. F. 1991, Nucl. Phys. A, 535, 331, doi: 10.1016/0375-9474(91)90452-C
  * Lentz et al. (2015) Lentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, ApJ, 807, L31, doi: 10.1088/2041-8205/807/2/L31
  * Mabanta & Murphy (2018) Mabanta, Q. A., & Murphy, J. W. 2018, ApJ, 856, 22, doi: 10.3847/1538-4357/aaaec7
  * Mabanta et al. (2019) Mabanta, Q. A., Murphy, J. W., & Dolence, J. C. 2019, ApJ, 887, 43, doi: 10.3847/1538-4357/ab4bcc
  * Marek et al. (2006) Marek, A., Dimmelmeier, H., Janka, H. T., Müller, E., & Buras, R. 2006, A&A, 445, 273, doi: 10.1051/0004-6361:20052840
  * Mathews et al. (2020) Mathews, G. J., Boccioli, L., Hidaka, J., & Kajino, T. 2020, Modern Physics Letters A, 35, 2030011, doi: 10.1142/S0217732320300116
  * Melson et al. (2015) Melson, T., Janka, H.-T., & Marek, A. 2015, ApJ, 801, L24, doi: 10.1088/2041-8205/801/2/L24
  * Miller et al. (1993) Miller, D. S., Wilson, J. R., & Mayle, R. W. 1993, ApJ, 415, 278, doi: 10.1086/173163
  * Müller et al. (2013) Müller, B., Janka, H.-T., & Marek, A. 2013, The Astrophysical Journal, 766, 43, doi: 10.1088/0004-637x/766/1/43
  * Müller (2016) Müller, B. 2016, PASA, 33, e048, doi: 10.1017/pasa.2016.40
  * Müller (2019) —. 2019, MNRAS, 487, 5304, doi: 10.1093/mnras/stz1594
  * Müller et al. (2012) Müller, B., Janka, H.-T., & Marek, A. 2012, ApJ, 756, 84, doi: 10.1088/0004-637X/756/1/84
  * Müller et al. (2016) Müller, B., Viallet, M., Heger, A., & Janka, H.-T. 2016, ApJ, 833, 124, doi: 10.3847/1538-4357/833/1/124
  * Müller et al. (2019) Müller, B., Tauris, T. M., Heger, A., et al. 2019, MNRAS, 484, 3307, doi: 10.1093/mnras/stz216
  * Murphy & Meakin (2011) Murphy, J. W., & Meakin, C. 2011, ApJ, 742, 74, doi: 10.1088/0004-637X/742/2/74
  * O’Connor (2015) O’Connor, E. 2015, ApJS, 219, 24, doi: 10.1088/0067-0049/219/2/24
  * O’Connor & Ott (2010) O’Connor, E., & Ott, C. D. 2010, Classical and Quantum Gravity, 27, 114103, doi: 10.1088/0264-9381/27/11/114103
  * O’Connor & Ott (2011) —. 2011, ApJ, 730, 70, doi: 10.1088/0004-637X/730/2/70
  * O’Connor et al. (2018) O’Connor, E., Bollig, R., Burrows, A., et al. 2018, Journal of Physics G: Nuclear and Particle Physics, 45, 104001, doi: 10.1088/1361-6471/aadeae
  * O’Connor & Couch (2018a) O’Connor, E. P., & Couch, S. M. 2018a, ApJ, 854, 63, doi: 10.3847/1538-4357/aaa893
  * O’Connor & Couch (2018b) —. 2018b, ApJ, 865, 81, doi: 10.3847/1538-4357/aadcf7
  * Olson et al. (2016) Olson, P. J., Warren, M., Meixner, M., et al. 2016, arXiv e-prints, arXiv:1612.08992.  https://arxiv.org/abs/1612.08992
  * Pan et al. (2019) Pan, K.-C., Mattes, C., O’Connor, E. P., et al. 2019, Journal of Physics G Nuclear Physics, 46, 014001, doi: 10.1088/1361-6471/aaed51
  * Pejcha & Thompson (2015) Pejcha, O., & Thompson, T. A. 2015, The Astrophysical Journal, 801, 90, doi: 10.1088/0004-637x/801/2/90
  * Perego et al. (2015) Perego, A., Hempel, M., Fröhlich, C., et al. 2015, ApJ, 806, 275, doi: 10.1088/0004-637X/806/2/275
  * Radice et al. (2018) Radice, D., Abdikamalov, E., Ott, C. D., et al. 2018, Journal of Physics G Nuclear Physics, 45, 053003, doi: 10.1088/1361-6471/aab872
  * Radice et al. (2016) Radice, D., Ott, C. D., Abdikamalov, E., et al. 2016, ApJ, 820, 76, doi: 10.3847/0004-637X/820/1/76
  * Schneider et al. (2019) Schneider, A. S., Roberts, L. F., Ott, C. D., & O’Connor, E. 2019, Phys. Rev. C, 100, 055802, doi: 10.1103/PhysRevC.100.055802
  * Sonneborn et al. (1987) Sonneborn, G., Altner, B., & Kirshner, R. P. 1987, ApJ, 323, L35, doi: 10.1086/185052
  * Steiner et al. (2013) Steiner, A. W., Hempel, M., & Fischer, T. 2013, ApJ, 774, 17, doi: 10.1088/0004-637X/774/1/17
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H.-T. 2016, The Astrophysical Journal, 821, 38, doi: 10.3847/0004-637x/821/1/38
  * Tamborra et al. (2014) Tamborra, I., Hanke, F., Janka, H.-T., et al. 2014, The Astrophysical Journal, 792, 96, doi: 10.1088/0004-637x/792/2/96
  * Thorne (1966) Thorne, K. S. 1966, ApJ, 144, 201, doi: 10.1086/148595
  * Ugliano et al. (2012) Ugliano, M., Janka, H.-T., Marek, A., & Arcones, A. 2012, ApJ, 757, 69, doi: 10.1088/0004-637X/757/1/69
  * Wilson & Mathews (2003) Wilson, J. R., & Mathews, G. J. 2003, Relativistic Numerical Hydrodynamics (Cambridge University Press)
  * Wilson & Mayle (1988) Wilson, J. R., & Mayle, R. W. 1988, Phys. Rep., 163, 63, doi: 10.1016/0370-1573(88)90036-1
  * Woosley & Heger (2007) Woosley, S. E., & Heger, A. 2007, Phys. Rep., 442, 269, doi: 10.1016/j.physrep.2007.02.009
  * Woosley & Weaver (1995) Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181, doi: 10.1086/192237