# Effect of the Nuclear Equation of State on Relativistic-Turbulence Induced
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  In-Saeng Suh
Center for Astrophysics, Department of Physics, University of Notre Dame, 225
Nieuwland Science Hall, Notre Dame, IN 46556, USA  Center for Research
Computing, University of Notre Dame, 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

The nuclear equation of state is an important component in the evolution of
core collapse supernovae. In this paper we make a survey of various equations
of state in the literature and analyze their effect on spherical core-collapse
models in which the effects of three-dimensional turbulence is modeled by a
general relativistic formulation of Supernova Turbulence in Reduced
dimensionality (STIR). We show that the viability of the explosion is quite
EOS dependent and that it best correlates with the early-time interior entropy
density of the proto-neutron star. We check that this result is not progenitor
dependent, although low-mass progenitors show different explosion properties,
due to the different pre-collapse nuclear composition. Larger central
entropies also induce more vigorous proto-neutron-star convection in our one-
dimensional turbulence model, as well as a wider convective layer.

##  1 Introduction

The detailed explosion mechanism for core-collapse supernovae (CCSNe) has been
one of the key problems in nuclear astrophysics for more than 50 years, and
resolving this problem still poses huge physical and computational challenges.
The thermodynamic conditions typical of the collapse and subsequent explosion
of massive stars are characterized by wide ranges of densities
($10^{3}-10^{15}$ g cm-3), temperatures ($1-100$ MeV) and electron fractions
($0.01-0.6$). Therefore, a detailed knowledge of the properties of matter
under all these thermodynamic conditions is required. However, the Equation of
State (EOS) for nuclear matter, especially at high densities, is still poorly
constrained (e.g. Klähn et al. (2006); Hebeler et al. (2013); Zhang et al.
(2018); Burgio et al. (2021)). In addition to this, an extremely large number
of neutrinos is produced during the collapse, some of which can be trapped
inside the newly formed Proto-Neutron Star (PNS). Therefore, an accurate
description of the interactions between neutrinos and regular matter is
needed, as well as advanced algorithms describing how they are transported
from the center to the outer layers of the star (Liebendörfer et al., 2005;
Kato et al., 2020; Mezzacappa et al., 2020).

In this article, we focus on the impact that the EOS can have on the explosion
mechanism of CCSNe by using the spherically symmetric, general relativistic
model described in Boccioli et al. (2021). In that paper, the Supernova
Turbulence In Reduced-dimensionality model (STIR), developed by Couch et al.
(2020), was extended to incorporate a general relativistic treatment of
turbulent convection. Then, it was studied how turbulent convection behaves
with and without taking General Relativity (GR) into account. Lastly, the
effect that this has on the explodability as a function of progenitor mass was
also analyzed. However, in that work only a single EOS was considered, and the
properties of the proto-neutron star (PNS) were not discussed. In the present
study, we extend that analysis by investigating the effects of the EOS on the
explosion and on the properties of the PNS.

There have been a number of studies in the past regarding the impact of the
nuclear equation of state on CCSNe (e.g. Lattimer & Swesty (1991); O’Connor &
Ott (2011); Hempel et al. (2012); Janka (2012); Steiner et al. (2013); Couch
(2013); Suwa et al. (2013); Fischer et al. (2014); Char et al. (2015); Olson
et al. (2016); Furusawa et al. (2017); Richers et al. (2017); Nagakura et al.
(2018); Morozova et al. (2018); Burrows et al. (2018); Schneider et al.
(2019); Harada et al. (2020); Yasin et al. (2020)), but in most of those
studies only one progenitor and a few EOSs were considered (although see the
recent paper by Ghosh et al. (2021), where they study the impact of the EOS on
the nucleosynthesis of several progenitors). Here, we study four progenitors
modeled with a broad range of EOSs, some of which were constructed using
Skyrme density functionals and others using Relativistic Mean Field (RMF)
theory.

There exist other recent studies that have analyzed the impact of the EOS on
the explosion properties and on the PNS structure (Schneider et al., 2019;
Yasin et al., 2020). However, they only focused on EOSs generated using Skyrme
energy-density functionals, since they were interested in identifying how
different nuclear parameters could impact the physics of CCSNe. In those
studies, they found that the effective nucleon mass $m_{n}^{*}$ is correlated
with the strength of the explosion. In this paper we take a different
approach, and compare EOSs calculated using different theoretical frameworks,
i.e. both Skyrme density functionals and relativistic mean field theory. In
this way we are able to analyze how differences in the underlying theoretical
framework translate into differences in explosion properties and PNS
structure.

This is the first extensive EOS study based upon a relativistic explosion
model driven by turbulent convection, and therefore we are also able to
analyze the impact that the EOS has on the convection in the gain region and
inside the PNS. In Section 2 we describe the EOSs adopted, the numerical setup
chosen and the essential physics of STIR. In Section 3 we describe how we
calibrate our parametric model by comparing to 3D simulations. Results
regarding the explosion properties, the structure of the PNS and neutrinos-
related quantities are reported in Section 4. Conclusions are given in Section
5. Throughout the manuscript, we adopt natural units, i.e. $G=c={\rm
M}_{\odot}=1$.

##  2 Methods

###  2.1 Equations of State

To study the effect of the EOS on the explosion properties, we simulated the
explosion of four representative progenitors with initial main-sequence masses
of 9, 15, 20 and 25 M⊙ from Sukhbold et al. (2016) using 7 different EOSs
summarized in Table 1. We employ some of the most common EOSs used for CCSNe
simulations to be consistent with previous results presented in the literature
(e.g. Lattimer & Swesty (1991); O’Connor & Ott (2011); Hempel et al. (2012);
Janka (2012); Steiner et al. (2013); Couch (2013); Suwa et al. (2013); Fischer
et al. (2014); Char et al. (2015); Olson et al. (2016); Furusawa et al.
(2017); Richers et al. (2017); Nagakura et al. (2018); Morozova et al. (2018);
Burrows et al. (2018); Schneider et al. (2019); Harada et al. (2020); Yasin et
al. (2020)), despite the fact that some of these equations of state are
inconsistent with recent observational and/or experimental constraints.

The nuclear properties that can most significantly impact the dynamics of
CCSNe and determine the properties of the proto-neutron star formed in the
collapse (Schneider et al., 2019; Yasin et al., 2020) are: (i) the zero
temperature saturation density of symmetric nuclear matter $n_{0}$; (ii) the
effective nucleon mass in symmetric matter at saturation density $m^{*}_{n}$;
(iii) the nuclear incompressibility $K_{0}$; (iv) the density-dependent
symmetry energy $\epsilon_{\rm sym}(n)$. The latter is defined as the iso-
vector component of the energy per baryon, obtained by expanding the energy
per baryon around its value for symmetric matter:

$\epsilon_{B}(n,y)=\epsilon_{B}(n,y=1/2)+\epsilon_{\rm sym}\delta(y)^{2}~{},$
(1)

where $y$ is the proton fraction and $\delta(y)=1-2y$. Hence, the symmetry
energy is defined as:

$\epsilon_{\rm sym}=\epsilon_{B}(n,y=1/2)-\epsilon_{B}(n,y=0).$ (2)

It is also common to encounter an alternative definition of the symmetry
energy, which can be derived by expanding $\epsilon_{B}(n,y)$ around its value
for symmetric matter (i.e. $\epsilon_{B}(n,y=1/2)$), yielding:

$\epsilon_{\rm sym}=\epsilon_{B}(n,y=1/2)+\epsilon_{\rm
sym}\delta(y)^{2}+\mathcal{O}[\delta(y)^{4}].$ (3)

The two definitions agree up to quartic terms.

A brief summary of these equations of state is as follows:

  * •

The APR EOS (Akmal et al., 1998; Schneider et al., 2019) is based upon a
thorough variational calculation, including realistic two- and three-body
nuclear forces  to derive the ground state energies of pure neutron matter and
symmetric nuclear matter. It utilizes a Hamiltonian that includes the two-body
Argonne V18 potential extracted by fitting two-nucleon experimental scattering
data, as well as the more empirical three-body Urbana IX potential, which is
fit in part to light nuclei. This is the only EOS considered here that
includes a phase transition to a neutral pion condensate. It is a
generalization of the original EOS at zero temperature from Akmal et al.
(1998), where the nuclear potentials are extended to the finite temperature
part. At low densities it is connected to a network of $\sim$3300 nuclei in
Nuclear Statistical Equilibrium (NSE) using the SRO code (Schneider et al.,
2017).

  * •

The LS220 (Lattimer & Swesty, 1991), KDE0v1 (Agrawal et al., 2005) and SLy4
(Chabanat et al., 1997) parameterization are taken from Schneider et al.
(2017). For these EoSs nuclear matter is described using Skyrme models for
nuclear interactions, and in the transition to nuclear saturation density
nuclei are described as a compressible liquid drop using the SNA. The version
of LS220 utilized here is in excellent agreement with the original work by
Lattimer & Swesty (1991). At low densities these equations of state are
coupled to an NSE network of $\sim$3300 nuclei using the SRO code (Schneider
et al., 2017).

  * •

The DD2 and SFHo EOSs are calculated using the Relativistic Mean-Field (RMF)
approach of Hempel & Schaffner-Bielich (2010) with the parametrization of the
interactions between nucleons from Typel et al. (2010) and Steiner et al.
(2013), respectively. Nuclear matter is treated using an NSE approach for the
transition from non-uniform to uniform nuclear matter. These two EOSs differ
from the EOSs described above in that they do not employ the SNA in the
vicinity of the saturation density. Rather, several thousand nuclei are taken
into account, whose masses and binding energies, when available, are taken
from experiment.

  * •

The HShen EOS from Shen et al. (2011) is also based on RMF theory. This EOS
differs from the SFHo and DD2 in that it uses a Single Nucleus Approximation
(SNA) based upon a Thomas-Fermi model to describe matter in the approach to
nuclear matter density.

The first four EOSs were generated using the SROEOS code from Schneider et al.
(2017). Therefore, they share the same treatment of inhomogeneous phases at
sub-nuclear densities. This consists of a finite temperature compressible
liquid-drop model in which heavy nuclei are described using a single-nucleus
approximation. All of the adopted equations of state in Table 1 can thus be
grouped in two categories: SRO-type (the first four), and the RMF-type (the
last three). The gravitational mass versus radius relationships for cold
neutron are shown in Figure 1, alongside some observational constraints on the
mass and radius of Neutron Stars. As can be seen, not all these EOSs are
consistent with the current observational constraints. Nevertheless, these
formulations should represent a reasonable gamut of possible equations of
state.

EOS | Type | $n_{0}$ | $J$ | $L$ | $K_{0}$ | $m^{*}_{n}/m_{\rm n}$ | $m^{*}_{p}/m_{\rm p}$ | Mmax
---|---|---|---|---|---|---|---|---
|  | (fm-3) | (MeV baryon-1) | (MeV baryon-1) | (MeV baryon-1) | - | - | M⊙
LS220 | SRO | 0.1549 | 28.61 | 73.81 | 219.85 | 1.0000 | 1.0000 | 2.04
APR | SRO | 0.1600 | 32.59 | 58.47 | 266.00 | 0.6980 | 0.6987 | 2.19
KDE0v1 | SRO | 0.1646 | 34.58 | 54.70 | 227.53 | 0.7440 | 0.7443 | 1.97
SLy4 | SRO | 0.1595 | 32.00 | 45.96 | 229.90 | 0.6950 | 0.6953 | 2.05
SFHo | RMF | 0.1583 | 31.57 | 47.10 | 245.40 | 0.7609 | 0.7606 | 2.05
DD2 | RMF | 0.1491 | 32.73 | 57.94 | 242.70 | 0.5628 | 0.5622 | 2.14
HShen | RMF | 0.1455 | 36.95 | 110.99 | 281.00 | 0.6340 | - | 2.21
Table 1: Parameters characterizing the equations of state utilized in this
study: $n_{0}$ is the nuclear saturation density, $J$ is the symmetry energy
at saturation density, $L$ is the slope of the symmetry energy, $K_{0}$ is the
nuclear incompressibility, $m_{n}^{*}$ is the nucleon effective mass and
$m_{n}$ is the nucleon mass. The first four equations of state are calculated
using the SRO code from Schneider et al. (2017), while the last three are
calculated using RMF theory.  Figure 1: Gravitational mass vs radius, from the
various equations of state adopted in the present work. The green shaded
region shows the NS mass-radius constraints from model A of Nättilä et al.
(2016). The grey contour represents the mass and radius of the millisecond
pulsar PSR J0030+0451 from Miller et al. (2019), while the blue contour
represents the mass and radius of the millisecond pulsar PSR J0740+6620 from
Miller et al. (2021), both measured by the NICER collaboration.

###  2.2 Numerical model

All of the simulations presented in this paper were run using a modified
version of the open-source code GR1D (O’Connor & Ott, 2010; O’Connor, 2015).
The original code is based on general relativistic, spherically symmetric
hydrodynamics and radiation transport. The Boltzmann equation that describes
the propagation of neutrinos is solved with a two-moment scheme, the so-called
M1 transport (Shibata et al., 2011; Cardall et al., 2013), while neutrino
opacities were generated using the open-source code NuLib (O’Connor, 2015).

We use the same spatial grid for all the simulations, i.e. 700 zones with a
constant resolution of 0.3 km in the inner 20 km and then logarithmically
spaced outer zones extending to a radius of $15\,000$ km. The only exception
are the progenitors with very steep density gradients (i.e. the ones with
masses smaller than 12 M⊙), for which we place the outer boundary at the
radius where the density reaches $6\times 10^{4}\ {\rm g\ cm}^{-3}$, which can
vary between $5\,000$ km and $10\,000$ km.

For the neutrinos, we use 18 energy groups logarithmically spaced from 1 to
280 MeV. We adopted the opacities described in Bruenn (1985), with weak-
magnetism and recoil corrections from Horowitz (2002) and a virial correction
to neutrino-nucleon scattering from Horowitz et al. (2017). We include
inelastic neutrino-electron scattering and velocity-dependent terms in the
transport equation according to O’Connor (2015).

###  2.3 Spherical relativistic turbulence model

Neutrino heated turbulence plays a crucial role in the shock revival and
explodability, as many previous studies have noted (Couch & Ott, 2015; Radice
et al., 2016, 2018; Mabanta & Murphy, 2018). In essence, material behind the
expanding shock is heated by neutrinos escaping the proto-neutron star. This
increases the entropy and leads to a negative entropy gradient which triggers
turbulent convection. The rising turbulent fluid elements can then transfer
energy to the shock. This makes spherical explosions triggered by neutrino
heated turbulent convection physically well-motivated.

To analyze the impact of each EOS on the explosion, we used a modified version
of GR1D described in Boccioli et al. (2021) based upon the Supernova
Turbulence In Reduced-dimensionality (STIR) model of Couch et al. (2020). With
that, one can include the effects of turbulent convection in our spherically
symmetric simulations by using a time-dependent Mixing-Length Theory (MLT)
approach, and therefore trigger an explosion.

The difference between our model and the original STIR is that we explicitly
take General Relativity (GR) into account. Despite yielding small differences
in the shock dynamics, GR can have a significant effect in the explodability
of SNe as a function of progenitor mass (Boccioli et al., 2021). Here we
summarize the main features of this model, but more details can be found in
Mabanta & Murphy (2018), Couch et al. (2020) and Boccioli et al. (2021). The
metric adopted in GR1D is the following:

$\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}$ (4)

where the $t-t$ component $\alpha$ is the lapse function and the $r-r$
component $X$ is a function of the gravitational mass of the system at radius
$r$.

In addition to the standard hydrodynamic equations solved in GR1D (O’Connor &
Ott, 2010), we have added an equation to account for the time evolution of the
turbulent energy:

$\begin{split}&\frac{\partial Dv_{\rm turb}^{2}}{\partial
t}+\frac{1}{r^{2}}\frac{\partial}{\partial r}[\frac{\alpha r^{2}}{X}D(v_{\rm
turb}^{2}v-D_{\rm K}\nabla v_{\rm turb}^{2})]\\\ &=-\alpha X\left(\rho v_{\rm
turb}^{2}\frac{\partial v}{\partial r}+\rho v_{\rm turb}\omega_{\rm
BV}^{2}\Lambda_{\text{mix}}-\rho\frac{v_{\rm
turb}^{3}}{\Lambda_{\text{mix}}}\right)~{},\end{split}$ (5)

where $v_{\rm turb}$ is the turbulent velocity and $D=X\rho W$ is the
conserved density, with $W=(1-v^{2})^{-1/2}$ being the Lorentz factor. The
proper physical velocity is defined as $v=Xv^{r}$, where $v_{r}$ is the
coordinate velocity, and $D_{\rm K}=\alpha_{\rm K}v_{\rm
turb}\Lambda_{\text{mix}}$ is a diffusion coefficient. The two main quantities
characterizing the model are the mixing length

$\Lambda_{\text{mix}}=\alpha_{\rm
MLT}\frac{P}{\rho~{}\text{d}\phi/\text{d}r}~{},$ (6)

and the Brunt-Väisälä frequency

$\omega_{\rm BV}^{2}=\frac{\alpha^{2}}{\rho
hX^{2}}\left(\frac{\text{d}\phi}{\text{d}r}-v\frac{\partial v}{\partial
r}\right)\left(\frac{\partial\rho(1+\epsilon)}{\partial
r}-\frac{1}{c_{s}^{2}}\frac{\partial P}{\partial r}\right)~{},$ (7)

where $\phi=\ln\alpha$ reduces to the Newtonian potential in the non-
relativistic limit. The quantity $h=1+\epsilon+P/\rho$ is the relativistic
enthalpy, where $\epsilon$ is the internal energy, $P$ is the pressure, and
$c_{s}$ is the sound speed.

The above equations depend on two parameters: $\alpha_{\rm MLT}$ and
$\alpha_{\rm K}$. Additional diffusion coefficients analogous to $D_{\rm K}$
are added to the evolution equations for internal energy, electron fraction
and neutrino energy density, increasing the number of parameters in the model
to 5. However, we fix $\alpha_{\rm K}$ and the other mixing length parameters
to $1/6$, consistent with the choice of Müller et al. (2016) and Couch et al.
(2020), since the model is not very sensitive to their value. The only
parameter left is $\alpha_{\rm MLT}$, which determines the magnitude of the
mixing length and therefore the scale of convection.

##  3 Calibration of the turbulence model

The size of the convective cells in our model depends upon the mixing length
$\Lambda_{\rm MLT}$, which is proportional to the mixing length parameter
$\alpha_{\rm MLT}$. The larger $\alpha_{\rm MLT}$, the larger the convective
eddies and the associated turbulent energy, which can therefore lead to a
stronger shock revival. Following the approach of Couch et al. (2020), we
calibrate the value of $\alpha_{\rm MLT}$ by comparing to realistic convection
in 3D simulations. As pointed out in Boccioli et al. (2021), the calibrated
value of $\alpha_{\rm MLT}$ can vary by up to a few percent depending upon
spatial resolution and the number of neutrino energy groups adopted.
Therefore, we maintain the same spatial resolution and number of neutrino
energy groups in all our simulations, as described in Section 2.2.

For the calibration of $\alpha_{\rm MLT}$ we utilize the 20 M⊙ solar-
metallicity progenitor of Farmer et al. (2016), and we employ the SFHo
equation of state. We then compare our simulations using 6 different values of
$\alpha_{\rm MLT}$ with the 3D model of Couch & O’Connor (2014), which
utilized the same progenitor, the same EOS and a very similar M1 scheme for
the neutrino transport.

A comparison of the turbulent velocities at $\sim$ 135 ms after bounce is
shown in Figure 2. As already pointed out in Couch et al. (2020), the 3D
profile exhibits a longer tail of convective motion between 50 and 75 km. This
results from non-radial motions. Indeed, the width of the gain region in 3D
simulations varies as a function of angle. As a consequence, the turbulent
velocity can extend down to smaller radii. Lastly, the 3D profile shows a
significant amount of convection in the PNS. This convection is present but
the convective velocity is much reduced in the 1D results. We will come back
to this in Section 4.2.

To complete the calibration, we have analyzed two more quantities shown in
Figure 3. The first is the shock radius as a function of time and the second
is the integrated turbulent energy at $\sim$ 135 ms after bounce.

Since convection in the PNS appears to be different in the 3D and 1D cases, we
do not include it in our calibration analysis. Hence, we define the integrated
turbulent energy as:

$E_{\rm turb}(r)=\int_{R^{*}}^{r}\rho v_{\rm turb}^{2}r^{\prime 2}{\rm
d}r^{\prime 2}~{},$ (8)

where we set $R^{*}=30$ km to exclude the convection in the PNS.

From the comparison shown in Figure 3, the value of $\alpha_{\rm MLT}$ that
most closely matches the 3D results appears to be in the range between 1.45
and 1.5. Specifically, for $\alpha_{\rm MLT}=1.48$, the integrated turbulent
energy at radii $r\gtrsim 150$ km is in excellent agreement with the 3D case
(see the right panel of Figure 3). In other words, at $\sim 135$ ms after
bounce, our simulation with $\alpha_{\rm MLT}=1.48$ exhibits the same total
turbulent energy in the gain region as the 3D case.

The excellent agreement with the total turbulent energy comes at the cost of a
slightly larger shock radius in our 1D simulation compared to the 3D at times
later than $\sim 250$ ms. The same discrepancy in the shock radius versus time
can be found in the original STIR simulations by Couch et al. (2020). However,
we choose to normalize the convection to the total turbulent energy at earlier
times ($\sim 135$ ms) as this is when convection sets in and the outcome of
the explosion is determined. One should also keep in mind that the shock
radius in our 1D model tends to be a little larger than in the 3D case in the
vicinity of the shock (see Figure 2). This is because the turbulent velocity
is only in the radial directions and therefore delivers energy to the outgoing
shock more efficiently. Nevertheless, we adopt this as the best representation
of the contribution of neutrino-heated turbulence to the shock even though
this can result in slightly larger shock radii for our 1D simulation.

Given the discrepancy between 1D and 3D results seen in some of the
hydrodynamical variables analyzed above, we also analyze how the value of
$\alpha_{\rm MLT}$ impacts the explodability as a function of progenitor mass.
We have simulated 58 progenitors from Sukhbold et al. (2016), with masses
ranging from 9 to 120 M⊙, using the SFHo EOS for three values of $\alpha_{\rm
MLT}$: 1.45, 1.48, and 1.5. The outcome of these simulations is shown in
Figure 4. The explodability with $\alpha_{\rm MLT}=1.48$ is consistent with
the previous findings reported in Boccioli et al. (2021).

With these results, one can calculate the total explosion fraction as a
function of $\alpha_{\rm MLT}$, weighed by a Salpeter initial mass function.
This is shown in Figure 5. For comparison, we also show the observationally
deduced range from Adams et al. (2017) (grey area) and the more recent results
of Neustadt et al. (2021) (green area).

The explosion fraction for $\alpha_{\rm MLT}=1.48$ lies just outside the 90%
confidence interval of the most recent results, while for $\alpha_{\rm
MLT}=1.5$ it is fully consistent with observations. However, since our
hydrodynamical quantities are better reproduced by $\alpha_{\rm MLT}=1.48$ and
the discrepancy with observational results is quite small, we adopt this value
for the rest of the paper, although the same discussion and results are
equally valid for simulations done with $\alpha_{\rm MLT}=1.5$.

After calibrating the value of $\alpha_{\rm MLT}$ we simulate four progenitors
with masses of 9, 15, 20 and 25 M⊙ for all of the EOSs listed in Section 2.1.

Figure 2: Turbulent velocity $v_{\rm turb}$ as a function of radius for
different values of $\alpha_{\rm MLT}$, taken at $\sim 135$ ms after bounce.
The dashed line represents the results from the 3D simulation of Couch &
O’Connor (2014), performed with FLASH. Figure 3: Panel (a) shows the shock
radius as a function of time post bounce for different values of $\alpha_{\rm
MLT}$. Panel (b) shows the integrated turbulent energy as a function of radius
for a snapshot at $\sim$135 ms post bounce. Notice that, given its definition
in Eq. 8, $E_{\rm turb}(\rm r)$ is zero below 30 km and then it’s constant for
radii larger than $\sim 150$ km, which is approximately the location of the
shock, above which no convection is developed. Figure 4: Explodability of
progenitors as a function of mass for different values of $\alpha_{\rm MLT}$.
Green bands represent successful explosions (defined as those simulations for
which the shock crosses a radius of 500 km), while black bands represent
failed explosions. We color in grey the progenitors we did not simulate.
Figure 5: The explosion fraction is calculated by weighing the explodability
as a function of progenitor mass shown in Figure 4 with a Salpeter initial
mass function. The black dots are the results from our simulations for
different values of $\alpha_{\rm MLT}$. The black and green horizontal lines
show the median from Adams et al. (2017) and Neustadt et al. (2021),
respectively. The shaded regions represent the 90% confidence interval. Notice
that in those papers they actually report the fraction of failed SN $f_{\rm
FSN}$, therefore here we assume that the fraction of explosion is $1-f_{\rm
FSN}$

##  4 Results

###  4.1 Explosion properties

Two of the main quantities of interest for an exploding supernova are the
trajectory of the shock as a function of time and the explosion energy. We
display these in Figure 6 for the four progenitors described in Section 1.
Since the simulations we performed extend to less than one second of the post
bounce phase, a calculation of the final total explosion energy is not
possible. Instead, one can define a diagnostic explosion energy (Buras et al.,
2006; Marek & Janka, 2009; Müller et al., 2012; Bruenn et al., 2016) as the
integral of the binding energy $E_{\rm bind}$ in regions where $E_{\rm
bind}>0$.

In the GR case, following Müller et al. (2012), one can define $E_{\rm bind}$
in terms of the lapse function $\alpha$, the density $\rho$, the specific
internal thermal energy $\epsilon_{\rm th}$, the pressure $P$ and the Lorentz
factor $W$:

$E_{\rm bind}=\alpha\left[(\rho+\rho\epsilon_{\rm th}+P)W^{2}-P\right]-\rho
W~{}.$ (9)

Here, it is important to distinguish between the traditional internal energy
adopted in classical thermodynamics and the internal energy used in the
context of CCSNe (see Appendix A of Bruenn et al. (2016)). The latter
definition incorporates the binding energy of nuclei in the expression of the
internal energy. This allows negative values of $\epsilon$, hence making the
definition of the total energy of the fluid problematic. A way around this is
to define $\epsilon_{\rm th}$ as the difference between the internal energy
$\epsilon$ — taken from the EOS at a given density, temperature and electron
fraction — and the internal energy $\epsilon_{0}$, defined as the internal
energy at zero temperature and at the same density and electron
fraction111From a practical point of view, what we do is to take the internal
energy corresponding to the lowest temperature of the table (typically $\sim
0.01$ MeV), which for the thermodynamic conditions of the supernova can be
safely treated as zero..

From this definition, it is clear that $\epsilon_{0}$ represents the binding
energy of the nuclei, which at low temperature is the only significant
contribution to the internal energy. Finally, the diagnostic explosion energy
can be defined as:

$E_{\rm diag}=\int_{E_{\rm bind}>0}E_{\rm bind}{\rm d}V~{},$ (10)

where d$V$ is the proper volume.

As can be seen from the shock radius and the diagnostic energy, the LS220 EOS
gives the earliest and strongest explosion, followed by the APR, KDE0v1, SLy4,
SFHo, DD2 and HShen. This hierarchy is maintained across all progenitors, with
the only exception being the 9 M⊙ progenitor, for which the hierarchy of the
explosion energies is different. Specifically, the SFHo and the DD2 generate
comparatively stronger explosions in the 9 M⊙ progenitor than in the more
massive progenitors, as can be seen from the diagnostic energy for this
progenitor.

In addition to this, it should be noted that the 9 M⊙ progenitor leads to very
early and faint explosions (the diagnostic energy in Figure 6 has been
multiplied by 8 so that it could more clearly be represented on the plot), and
the decrease around 450 ms post bounce is a numerical artifact due to the
shock leaving the outer boundary of the simulation. This trend for the 9 M⊙
progenitor relates to the fact that these stars explode as electron capture
supernovae (Isern et al., 1991). A more detailed analysis of the differences
between this progenitor and the more massive ones is given in Appendix A. For
now we will focus on the 15, 20 and 25 M⊙ progenitors.

Since our convection model is designed to reproduce turbulence in the gain
region, one can analyze how the explosion is affected by this mechanism. To
make the comparison among different EOSs, which have different gain and shock
radii, we define a dimensionless quantity:

$\xi_{\rm gain}=\frac{r-r_{\rm gain}}{r_{\rm shk}-r_{\rm gain}}~{},$ (11)

such that $\xi_{\rm gain}=0$ when $r=r_{\rm gain}$ and $\xi_{\rm gain}=1$ when
$r=r_{\rm shock}$. The profiles of the turbulent velocity in the gain region
are plotted in Figure 7 for two different time snapshots.

At $\sim 50$ ms post bounce the gain region forms and the turbulent velocity
starts building up. After that, as can be seen in the snapshot at $\sim 150$
ms, the profiles become smoother and the peak velocity steadily grows.
Overall, the turbulent velocity profiles do not significantly change with
progenitors and EOSs.

Nevertheless, some global trends can be identified, since for example, the
LS220 and the SFHo consistently yield the largest peak turbulent velocities,
while the smallest corresponds to the SLy4. By looking at Figure 6, however,
it is clear that the strongest explosions are given by the LS220 and the APR,
which yields a smaller peak turbulent velocity than the SFHo. Hence, there
doesn’t seem to be any direct correlation between the strength of convection
in the gain region and the strength of the explosion as a function of EOS.

It is interesting to note that if one only considers SRO-type EOSs or only
RMF-type EOSs, the hierarchy of turbulent velocities follows the hierarchy of
explosion strengths. Namely, for SRO-type EOSs, from largest to lowest peak
turbulent velocity, we find: LS220, APR, KDE0v1, SLy4, which is the same as
the hierarchy for explosion strengths. For RMF-type EOSs the hierarchy of peak
turbulent velocities reads: SFHo, DD2, HShen, which is also the same as the
hierarchy of explosion strength.

Putting SRO and RMF types together, however, disrupts the correlation between
explosion strength and convective velocity, as can be seen from the zoomed-in
panels of Figure 7. There, one can see that SFHo, DD2 and HShen yield a quite
large turbulent velocity compared to the SRO-type EOSs, despite having low
explosion energies, or not exploding at all. We will come back to this when we
discuss neutrino luminosities in Section 4.4. We suggest that this difference
in the turbulent velocity profiles (and neutrino luminosities, as discussed
later) between SRO-type and RMF-type EOSs lies most likely in the different
treatment of nuclear matter as it approaches saturation density.

Figure 6: The upper panels display shock radius vs time for different
progenitors and EOSs, while the lower panels display the diagnostic energy vs
time. Notice that for the 9 M⊙ progenitor the diagnostic energies have been
multiplied by a factor of 5 so that they could be more clearly shown on the
same scale used for the other progenitors. Figure 7: Turbulent velocity
$v_{\rm turb}$ in the gain region for different progenitors (columns) and EOSs
(colors). The horizontal axis is the “normalized radius” defined in Eq. 11
which is 0 at $r_{\rm gain}$ and 1 at $r_{\rm shk}$. The upper panels show
profiles of turbulent velocity taken at $~{}50$ ms after bounce, when
convection sths developing, while the lower panels are taken at $~{}150$ ms
after bounce, when convection is fully developed.

###  4.2 PNS interior and the role of central entropy

The PNS is a very high density environment, and therefore its interior can be
greatly impacted by the EOS. In Figure 8 we show the central density,
temperature and entropy per baryon as a function of time for different
progenitors and EOSs. Note that the central entropy in the first 50 ms after
bounce correlates extremely well with the strength of the explosion. Indeed,
by comparing it to the explosion energy and shock trajectory, one can clearly
see that the EOS with the highest central entropy (LS220) yields the earliest
and strongest explosion, while the EOS with the lowest central entropy (HShen)
does not explode at all, and after the initial stalling of the shock at $\sim
150$ km is the EOS with the quickest re-collapse. In general, the hierarchy of
central entropies is the same as the hierarchy of explosion strength for all
progenitors, with the possible exception of the SLy4 and KDE0v1 EOSs, which we
will discuss at the end of this Section.

As a historical note, this correlation between explodability and early-time
central entropy was suggested by Hans Bethe (Bethe, 1990). The explanation
being the fact that entropy in the core is generated by the onset of
nonequilibrium processes due to deleptonization during collapse and energy
exchange between matter and the non-thermal neutrinos. The amount of entropy
generated will then affect the shock formation radius, and hence, the
explodability. For example, the lower central $Y_{e}$ and temperature for the
LS220 EoS are consistent with stronger deleptonization and larger entropy.

Another way to illustrate the same concept, from the point of view of MLT, is
that the minimum value of $\alpha_{\rm MLT}$ needed to trigger an explosion is
very low (typically around 1.2/1.3, depending on the progenitor) for the EOS
with the highest central entropy (LS220), while it is very large (typically
around 1.7/1.8, depending on the progenitor) for the EOS with the lowest
central entropy (HShen). On the other hand, central temperature and density do
not correlate with the strength of explosion.

Figure 8: The first three rows show the main central thermodynamic variables
for different progenitors (columns) and EOSs (colors) as a function of time
post bounce. The last row shows the ratio of the nucleon effective mass to the
nucleon mass calculated for the central values of densities and temperatures
as a function of time post bounce. The central entropy at early times tpb
$\lesssim 5$ ms is larger for EOSs that yield stronger explosions. See the
text for a more detailed discussion.

A similar relationship was found in the studies of Schneider et al. (2019) and
Yasin et al. (2020), who focused only on one type of EOS, i.e. a liquid drop
model which uses Skyrme interactions between nucleons at high densities. With
their setup they were able to vary different parameters (such as the effective
mass of nucleons at saturation density, the symmetry energy, the
incompressibility etc…) one by one, while keeping all the other parameters
fixed. Therefore, they were able to pinpoint which property of the EOS has the
largest impact on the explosion.

Their analysis revealed that the effective mass of nuclei at saturation
density $m^{*}_{n}$ is the quantity that most impacts the structure of the
PNS, and therefore the explosion, since it determines how fast its contraction
is (Yasin et al., 2020) and changes the temperature of the neutrinospheres for
all flavors (Schneider et al., 2019). They also pointed out that the effective
mass affects the central entropy of the PNS, in that a larger effective mass
would lead to a larger central entropy.

Our approach is somewhat different compared to these previous studies. We
include not only EOSs based on Skyrme interactions, like the LS220, KDE0v1,
and SLy4, but also the APR EOS and EOSs based upon RMF theory. The APR EOS
shares the same properties as the Skyrme-type EOSs near saturation density,
therefore we group them together under the label “SRO-type”. However, the
nuclear potentials for the APR EOS are derived from nucleon-nucleon
scattering, rather than from the energy density of nuclear matter. Alongside
the EOS developed by Togashi et al. (2017), it is one of the few EOSs
available for supernovae simulations developed starting from realistic nuclear
potentials. The EOSs calculated using RMF theory use either an NSE (SFHo and
DD2) or a Thomas-Fermi (HShen) treatment of nuclei near saturation density. By
considering this variety of EOS formulations we can investigate the subtle
differences among the most common frameworks used to calculate the nuclear
EOS.

When we compare EOSs calculated using Skyrme interactions (i.e. LS220, KDE0v1
and SLy4), we find that the effective mass indeed correlates with the strength
of explosion. However, this is not true when we compare all of the RMF-type
and SRO-type EOSs, since for example the SFHo has a larger effective mass than
the SLy4, but generates weaker and somewhat delayed explosions. Also, despite
the APR sharing the same properties as the Skyrme-type EOSs near saturation
density, it does not fit into the correlation between effective mass and
strength of explosion. For example, it has a smaller effective mass than
KDE0v1 but yields stronger explosions. This is not completely unexpected since
the density dependence of the effective mass is more complicated in the APR
model (Schneider et al., 2019). It also shows that the effective mass is not
the only parameter that can predict the strength of explosion. Instead, the
only quantity that seems to correlate with the strength of the explosion is
the central entropy, regardless of the framework used to calculate the nuclear
EOS.

A slight deviation from the correlation between explosion and central entropy
is represented by the case of KDE0v1 and SLy4. These two EOSs yield shock
trajectories and explosion energies that are very similar across progenitors,
with KDE0v1 always giving a slightly stronger explosion. If one looks at the
central entropy right after bounce, however, KDE0v1 gives a slightly smaller
value than SLy4. This shows that, despite $m^{*}_{n}$ being the parameter that
impacts the PNS structure the most, not surprisingly the other nuclear
parameters can also change the structure of the PNS. In particular, assuming a
Fermi liquid theory, one expects the central entropy to be roughly given by
$S_{c}\sim m^{*}_{n}T_{c}/\rho_{c}^{2/3}$ (Baym & Pethick, 1991), where
$S_{c}$, $T_{c}$ and $\rho_{c}$ are the central entropy, temperature and
density, respectively. KDE0v1 has a larger nucleon effective mass $m^{*}_{n}$,
which increases the central entropy, but it also has a larger saturation
density, which has the effect of increasing the central density, hence
lowering the central entropy. These competing effects can slightly disrupt the
hierarchy of central entropy, but since in our analysis the EOSs differ for
several nuclear parameters, we don’t expect a perfect correlation.
Nonetheless, the correlation between central entropy right after bounce and
the strength of the explosion remains quite robust.

It is not completely unexpected that the correlation between strength of
explosion and nucleon effective mass breaks down when comparing SRO-type and
RMF-type EOSs. The definition of effective mass and entropy is different in
the two frameworks (and within SRO-type, it varies between the APR and the
Skyrme-type EOSs). This is especially true since relativistic effects, which
are taken into account by RMF models but not by the SRO models, can be
significant. As can be seen in Figure 9, however, for densities between
$10^{14}$ and $10^{15}$ g cm-3 the hierarchy of effective masses is the same
as the hierarchy of entropies. At first, this seems to be in contrast with the
fact there is no clear correlation between effective nucleon mass at
saturation density and central entropy in our simulations. However, from
Figure 8 one can see that different EOSs yield very different thermodynamic
conditions inside the PNS. Temperatures, densities and electron fractions span
ranges of 10-20 MeV, $3-4\times 10^{14}$ g cm-3, and 0.27-0.3, respectively.
Therefore, comparing different EOSs at a fixed temperature and electron
fraction can be misleading if one wants to analyze the impact of the EOS on
the explosion of CCSNe. The bottom panel of Figure 8 confirms this. Here, one
can see that the hierarchy of the nucleon effective masses at the center of
the newly formed PNS is different than the hierarchy of entropies, and
therefore does not correlate with the strength of explosion.

Figure 9: Panel (a) and (b) show the entropy and the effective nucleon mass,
respectively, as a function of density for the EOSs considered in this paper.
Temperature and electron fraction are 15 MeV and 0.25, respectively, and
represent typical thermodynamic conditions reached in core-collapse supernovae
environments. The nucleon effective mass is the average of the proton and
neutron effective masses, which for most EOSs are different, while for HShen
and LS220 are equal.

###  4.3 PNS interior and the role of convection

Convection within the PNS has been shown to be very important since it can
change the neutrino signal (Roberts et al., 2012; Mirizzi et al., 2016;
Müller, 2020), the gravitational wave signal (Marek et al., 2009; Yakunin et
al., 2010; Müller et al., 2013; Andresen et al., 2017; Morozova et al., 2018)
and trigger the so-called Lepton-number Emission Self-sustained Asymmetry
(LESA) (Tamborra et al., 2014). However, its impact on the explosion is still
unclear, and more high-resolution 3D simulations are needed to clarify its
role.

In a recent paper, Nagakura et al. (2020) (hereafter NBRV20) analyzed the
convection inside the PNS for 14 progenitors from Sukhbold et al. (2016), with
masses ranging from 9 to 60 M⊙, using the SFHo EOS. Their conclusion was that
there doesn’t seem to be any direct influence of the shock revival on the PNS
convection at early times (and vice-versa). The most significant correlation
they found is between more massive PNSs and more vigorous convection.

The quantity chosen to represent the “vigor” of convection is the total
turbulent energy in the PNS, defined as:

$E^{\rm pns}_{\rm turb}=\int_{0}^{R_{\rm pns}}\rho v_{\rm turb}^{2}dV~{}.$
(12)

The top panels in Figure 10 show the evolution of $E^{\rm pns}_{\rm turb}$ for
the various progenitors and equations of state.

In Figure 11 we show the time evolution of baryonic mass of the PNS for the
various equations of state and progenitors considered in this paper. This
figure, together with Figure 10, confirms the trend found by NBRV20, i.e. that
more massive PNSs show more vigorous convection. This is true for all EOSs,
hence we confirm that the trend noted in the simulations of NBRV20 holds true
regardless of the EOS used.

Figure 10: The upper panels show the total turbulent energy in the PNS as a
function of time post bounce, as defined in Eq. 12; the central panels show
the width of the convective layer inside the PNS as a function of time post
bounce; The bottom panels show the baryonic mass contained in the convective
layer. Different columns represent different progenitors, while different
colors represent different EOSs.

It is also interesting to analyze each progenitor separately, and study how
$E^{\rm pns}_{\rm turb}$ correlates with the strength of the explosion. For
each progenitor, as can be seen in Figure 10, the hierarchy of $E^{\rm
pns}_{\rm turb}$ is the same as the hierarchy of the explosion energy, i.e.
earlier explosions with larger explosion energies have a more vigorous PNS
convection. From the analysis carried out in Section 4.2, one could then
conclude that a larger central entropy leads to stronger convection, since
both are correlated with a stronger explosion.

Figure 11: The baryonic mass of the PNS as a function of time post bounce.
Here, different colors refer to different EOSs, while different line styles
refer to different progenitors. The apparent correlation between PNS mass and
progenitor mass is only a coincidence. If more progenitors were included one
would not see this correlation. The real correlation is the one discussed in
the text between the PNS mass and turbulent energy.

Another detail that emerges from our analysis is that the SRO-type EOSs seem
to generate overall stronger PNS convection and form larger convective layers.
By looking at the top two panels of Figure 10, it is evident how there is a
clear separation between the curves which represent SRO-type EOSs and the
curves that represent RMF-type EOSs.

We postpone a more thorough analysis of PNS convection and this effect in
particular to future work. The reason for this is that the turbulent
convection model that we used is designed to reproduce convection in the gain
region, and therefore it is not accurate inside the PNS. In fact, by looking
at the turbulent velocity near the PNS in Figure 2, it is clear that there is
a discrepancy with the 3D results. Overall, the turbulent velocities we find
are roughly a factor of 10-20 smaller than in 3D. Hence, the turbulent energy
in the PNS $E^{\rm pns}_{\rm turb}$ is also about two orders of magnitude
smaller than the expectation from 3D simulations.

As pointed out by Couch et al. (2020), the convection inside the PNS is
extremely efficient in STIR, and therefore it quickly erases the unstable
thermodynamic gradients. This shuts off further growth of the turbulent
velocity, which explains why in STIR it is so small. This is confirmed by the
fact that, counterintuitively, smaller values of $\alpha_{\rm MLT}$ exhibit
larger turbulent velocities inside the PNS. This indicates that, for large
values of $\alpha_{\rm MLT}$, convection is too efficient to allow the
turbulent velocity to build-up.

Another reason why STIR is inadequate to predict the PNS convection is that it
does not take lepton gradients into account. The neutrinos trapped inside the
PNS alter the full lepton gradient. This is not merely equal to the electron
fraction gradient as it is in the gain region, and as it is assumed in our
STIR model.

Calculating the full lepton gradient, however, is not trivial (Roberts et al.,
2012) since it involves complicated thermodynamic derivatives. Hence, we leave
a more detailed treatment of PNS convection to a future work.

Despite the inadequacy of STIR in predicting the magnitude of the turbulent
velocity in the PNS, we find that it correctly predicts the width and the mass
of the PNS convective layer, when compared to the results of NBRV20, and we
find the same correlation between the mass of the PNS and the strength of
convection.

###  4.4 EOS dependence of the neutrino signal

Figure 12: The top, middle and bottom panels show the radius, density and
temperature of the neutrinosphere as a function of time post bounce for the 20
M⊙ progenitor and the EOSs listed in Table 1, with $\alpha_{\rm MLT}=1.48$.
The curves have been smoothed out, since the sharp thermodynamic gradients
lead to large jumps in density and temperature between adjacent zones, hence
resulting in a saw-tooth pattern. Different columns refer to different
neutrino flavors, i.e. electron neutrino $\nu_{e}$, electron antineutrino
$\bar{\nu}_{e}$ and heavy lepton neutrino $\nu_{x}$. The densities of the
electron and heavy-lepton neutrinos have been multiplied by 2 and 0.5,
respectively. Figure 13: The top and bottom panels show the average energy and
the luminosity of neutrinos, respectively, for the 20 M⊙ progenitor, with
$\alpha_{\rm MLT}=1.48$, and the EOSs listed in Table 1. Both quantities are
taken at 500 km from the center of the star, where neutrinos no longer
interact with matter, i.e. are free streaming. The luminosity of the heavy-
lepton neutrinos $\nu_{X}$ has been multiplied by a factor of 2 to present
them on the same scale as the other neutrinos.

Since the pioneering work of Bethe & Wilson (1985), it has been clear that
delayed neutrino heating is a primary mechanism enabling the explosion. The
EOS changes the thermodynamic conditions of the PNS, and therefore it
indirectly changes the neutrino-matter interaction. This in turn changes the
neutrino properties, such as the luminosity evolution and the energy spectrum.
The bulk of emission of neutrinos happens inside the PNS, where neutrinos
decouple from matter. This is usually identified as the region in the vicinity
of the neutrinosphere, defined as the layer inside the PNS at which the
neutrino opacity is equal to $2/3$. Therefore, changes in temperature and
density inside the PNS will affect the temperature and density profiles where
the neutrinos decouple, thereby modifying the emission properties of
neutrinos.

As discussed in Schneider et al. (2019), a larger value of the effective mass
at saturation changes the properties of neutrinos, since it causes a faster
contraction of the PNS. This shifts the location of the different
neutrinospheres to smaller radii, larger temperatures and lower densities. The
only exception are the heavy-lepton neutrinos. In this case the density of the
neutrinosphere increases as the effective mass increases. It should also be
noted that neutrinos of different energy decouple at different radii, since
the opacity depends upon the energy-dependent neutrino interaction cross
section, but for simplicity one usually refers to a single neutrino-averaged
neutrinosphere, defined using an energy-averaged opacity.

To analyze the impact of the EOS on the neutrino properties, we focus on the
20 M⊙ progenitor, although a similar behavior can be seen in the other
progenitors as well. The location and thermodynamic properties of the
neutrinosphere for different flavors are shown in Figure 12, and the
luminosities and average energies are shown in Figure 13.

If we limit the analysis to the SRO-type EOSs, the trends found by Schneider
et al. (2019) are confirmed in our simulations. That is, the SRO-type EOSs
with a larger effective mass yield a hotter neutrinosphere located at smaller
radii. Similarly, larger effective masses lead to larger neutrino luminosities
and energies for all flavors. If we include the RMF-type EOSs in this
analysis, however, these trends are disrupted, and for example the SFHo, with
a smaller effective mass at saturation than the LS220, has larger
neutrinosphere temperatures, neutrino energies and luminosities for all
flavors.

In their paper, Schneider et al. (2019) show the neutrino properties for non-
exploding models. In the present work, however, some of the EOSs lead to
explosions (due to the turbulent convection), and therefore some neutrino
properties will be different. In particular, as the explosion sets in, the
accretion onto the PNS stops, and the neutrino emission leads to fast cooling
of the inner regions. This accelerates the contraction of the PNS and moves
the neutrinospheres to regions of higher density. In addition to this, the
temperature stops increasing and becomes roughly constant. As a consequence,
the average neutrino energy also tends to a constant value, while the
luminosity sharply decreases, since without accretion the neutrino production
rapidly diminishes. Hence, in these simulations, the EOS is responsible for
differences in the neutrino signal only in the first 150-200 ms, when
turbulent convection hasn’t yet triggered an explosion. After that, the
different explosion dynamics will dramatically change the neutrino signal.

Figure 14: Electron neutrino luminosity for the 20 M⊙ progenitor and the EOSs
listed in Table 1, with $\alpha_{\rm MLT}=0$ (i.e. no turbulent convection
present). Highlighted is the accretion phase when the luminosity becomes
constant. The hierarchy of luminosities is the same for electron neutrinos and
antineutrinos, while the EOS effects on the heavy lepton neutrinos are
smaller.

As pointed out in previous studies (Hempel et al., 2012; Steiner et al., 2013)
the EOS can impact the neutrino signal in several ways. For example, the
presence of light nuclei near the neutrinospheres could in principle modify
the neutrino emission and absorption rates, and therefore affect the neutrino
signal. However, since we use the standard Bruenn (1985) neutrino opacities,
which don’t include interactions with light nuclei, we do not see this effect
in our simulations. The main quantities that can affect the neutrino
luminosities and average energy are mass accretion and the location of the
neutrinosphere. Larger mass accretions result in larger luminosities, whereas
neutrinospheres located at smaller radii (and therefore higher temperatures)
yield larger neutrino energies. A more detailed analysis of these effects can
be found in Hempel et al. (2012).

It is also worth mentioning that, given the different treatment of nuclei near
saturation density (namely, SNA or NSE), the fraction of heavy nuclei in
inhomogeneous phases will be different in each EOS. This can also alter the
neutrino signal, since a different treatment of heavy nuclei at intermediate
densities will change the electron fraction (and therefore the
deleptonization) of the core. Smaller electron fractions lead to more compact
configurations, and therefore larger luminosities and energies (for a more
detailed description, see Hempel et al. (2012)). The contraction of the PNS
can also be influenced by the presence of heavy nuclei below saturation
density, as well as by the inclusion of general relativistic effects (taken
into account only by RMF-type EOSs). Fast contractions will quickly increase
the temperatures of the neutrinosphere, and therefore increase the average
energies of neutrinos. A more detailed discussion can be found in Hempel et
al. (2012) and Steiner et al. (2013).

Overall, our simulations don’t show a clear correlation between neutrino
properties and strength of explosion. However, if one analyzes the neutrino
luminosity during the accretion phase in models with no turbulence (i.e.
simple spherically symmetric simulations without the inclusion of STIR) one
can see the correlation between effective mass and luminosity for SRO-type
models (see Figure 14). On the other hand, the RMF-type models show
comparatively larger luminosities, although the hierarchy (SFHo, DD2, HShen)
is the same as the hierarchy of explosion strength. Hence, we see that within
SRO-type or RMF-type frameworks the correlation between luminosity and
explosion strength holds. However, this is disrupted when analyzing EOSs
across both frameworks, as already pointed out in Section 4.1 for the case of
the turbulent velocity in the gain region.

One can conclude that, compared to SRO-type EOSs, RMF-type EOSs have larger
luminosities and turbulent velocities in the gain region but generate weaker
explosions. The reason behind this most likely lies in the different treatment
of nuclei near saturation density. A similar difference can be found in the 9
M⊙ progenitor, discussed in Appendix A, where a different treatment of nuclear
matter near saturation density leads to larger central densities which then
modify the explosion energies. Hence, not only the nuclear properties, but
also the treatment of nuclear matter near saturation density can affect the
explosion of core-collapse supernovae.

##  5 Conclusions

In this article we have shown how the EOS can affect the explosion of CCSNe
using one-dimensional, fully general relativistic simulations that employ a
parametric treatment of relativistic turbulent convection. We calibrated the
main free parameter of our model, $\alpha_{\rm MLT}$, by comparing our results
to 3D simulations. We then simulated the collapse of four stellar progenitors
using 7 different EOSs.

We find a remarkable correlation between the strength of the explosion and the
central entropy immediately after bounce. This is in agreement with previous
results obtained by Schneider et al. (2019) and Yasin et al. (2020). In these
works, EOSs calculated using Skyrme interaction models were used to analyze
the impact of different nuclear parameters — such as the effective nucleon
mass, the symmetry energy and the incompressibility — on the explosion. They
concluded that the effective nucleon mass is the parameter that most strongly
correlates with the explosion, as well as correlating with central entropy.
Our results using the three Skyrme EOSs considered in this work confirm this.

In addition to Skyrme EOSs, however, we have added the APR EOS and three RMF-
type EOSs to the comparison. This alters the correlation between strength of
explosion and effective nucleon mass. If one considers entropy and effective
nucleon mass for densities between $10^{13}$ and $10^{15}$ g cm-3 and for a
fixed temperature and electron fraction, larger nucleon effective masses are
correlated with larger entropies. However, in CCSNe simulations, different
EOSs yield different central temperatures and densities, and therefore the
correlation between effective nucleon mass and entropy breaks down. Similarly,
the correlation between effective nucleon mass and strength of explosion
reported in other studies (Schneider et al., 2019; Yasin et al., 2020) is no
longer present in our simulations. However, we still obtain a remarkable
correlation between strength of the explosion and central entropy. Hence, we
conclude that the central entropy immediately after bounce is the best
indicator of explodability, and must therefore play a key role in determining
the strength of the explosion.

We also analyzed how the central entropy is correlated with convection in the
PNS. We found that more vigorous convection induces larger convective zones
and is correlated with a larger central entropy. Moreover, the SRO-type EOSs
generate much larger convective zones and integrated turbulent energies than
RMF-type EOSs. The role of the PNS turbulent convection should be more
thoroughly investigated, possibly using non-parametric 2D and 3D simulations.
Our current PNS convection, however, shows some discrepancies with 3D results,
mainly in the radial profiles of the convective velocity. We therefore leave
an improved treatment of the PNS convection to a future work.

Finally, we analyzed the neutrino signal for various EOSs. We did not find any
clear correlation between the strength of the explosion and neutrino
luminosities or energies. However, we noticed that the RMF-type EOSs tend to
produce larger neutrino luminosities and turbulent velocities in the gain
region compared to the SRO-type EOSs. This is true even despite them yielding
weaker explosions. Therefore, we conclude that the approach used to calculate
the EOS, and more specifically the treatment of nuclear matter near saturation
density, has a significant impact on the explosion of CCSNe. ¡

The authors would like to thank Sean Couch and Mackenzie Warren for fruitful
discussions, Andre da Silva Schneider for his help in calculating the symmetry
energy of the EOSs adopted, and Erika Holmbeck for helping us plotting the
NICER constraints. 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 Discussion on 9 M⊙ progenitor

A separate discussion for the 9 M⊙ progenitor is warranted since it does not
follow the same correlations discussed in the text regarding explosion energy,
central entropy and strength of the PNS convection. In particular, the DD2 and
SFHo EOSs yield smaller central entropies and larger shock radii than what
would be expected from the trends found in the other progenitors (see Figure
15). To understand the origin of this discrepancy, it’s useful to define the
quantity.

$\Delta\rho_{\rm bounce}=\frac{\rho_{\rm bounce}^{9M_{\odot}}-\rho_{\rm
bounce}^{20M_{\odot}}}{\rho_{\rm bounce}^{20M_{\odot}}}~{}.$ (A1)

This represents how much larger (or smaller) the density is at bounce than
that in the heavier, e.g. 20 M⊙, progenitor. We arbitrarily chose to compare
to the 20 M⊙ progenitor, although the same conclusions apply if one compares
to either of the other progenitors studied in this work.

Figure 15: The vertical axis shows the radius of the shock at 300 ms post
bounce. We use this as a proxy for the strength of the explosion, since the
diagnostic energy is non-zero only for exploding progenitors. The larger the
shock radius, the closer the star is to an explosion. On the horizontal axis
we show the value of the central entropy at 5 ms after bounce, and show that
these two quantities are very well correlated (see the text for a more
detailed discussion). Different panels represent different progenitors, while
different colors represent different EOSs. Note that the shock radius for the
9 M⊙ progenitor is multiplied by a factor of (1/3) to fit on the plot. Figure
16: Relative difference between the density at bounce for the 9 M⊙ and the 20
M⊙ progenitors, as defined in Eq. (A1), as a function of radius. In the 9 M⊙
progenitor, the DD2 and SFHo EOSs yield a much larger central density than the
other EOSs, for example compared to the 20 M⊙ progenitor. This is due to the
different treatment of nuclear matter near the saturation density, i.e. using
an NSE approach versus the SNA model.

As can be seen in Figure 16, the increase of the central density at bounce in
the 9 M⊙ progenitor with the DD2 and SFHo EOSs is much larger than with the
other EOSs, which show an increase (or decrease, for the HShen) of up to a few
percent. Both the DD2 and SFHo EOSs were calculated using an NSE-approach near
saturation density. Since low-mass stars are characterized by a different core
composition (mostly 54Fe, 40Ca and 44Ti, versus heavier iron-peak elements in
higher-mass stars) the NSE approach changes the thermodynamic properties
significantly compared to the other EOS models.

Overall, the explosion energy for the 9 M⊙ progenitor is much smaller than
that of the higher mass progenitors. This is because the density drops off
much more rapidly for low-mass stars. Hence, when the material in the outer
layers is ejected, its density is very small. As a consequence the internal
and kinetic energies are decreased, making the total explosion energy smaller
than that of the higher mass stars.

## References

  * Adams et al. (2017) Adams, S. M., Kochanek, C. S., Gerke, J. R., & Stanek, K. Z. 2017, MNRAS, 469, 1445, doi: 10.1093/mnras/stx898
  * Agrawal et al. (2005) Agrawal, B. K., Shlomo, S., & Au, V. K. 2005, Phys. Rev. C, 72, 014310, doi: 10.1103/PhysRevC.72.014310
  * Akmal et al. (1998) Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804, doi: 10.1103/PhysRevC.58.1804
  * Andresen et al. (2017) Andresen, H., Müller, B., Müller, E., & Janka, H. T. 2017, MNRAS, 468, 2032, doi: 10.1093/mnras/stx618
  * Baym & Pethick (1991) Baym, G., & Pethick, C. 1991, doi: https://doi.org/10.1002/9783527617159
  * Bethe (1990) Bethe, H. A. 1990, Rev. Mod. Phys., 62, 801, doi: 10.1103/RevModPhys.62.801
  * Bethe & Wilson (1985) Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14, doi: 10.1086/163343
  * Boccioli et al. (2021) Boccioli, L., Mathews, G. J., & O’Connor, E. P. 2021, ApJ, 912, 29, doi: 10.3847/1538-4357/abe767
  * Bruenn (1985) Bruenn, S. W. 1985, ApJS, 58, 771, doi: 10.1086/191056
  * 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
  * Buras et al. (2006) Buras, R., Rampp, M., Janka, H. T., & Kifonidis, K. 2006, A&A, 447, 1049, doi: 10.1051/0004-6361:20053783
  * Burgio et al. (2021) Burgio, G. F., Schulze, H. J., Vidaña, I., & Wei, J. B. 2021, Progress in Particle and Nuclear Physics, 120, 103879, doi: 10.1016/j.ppnp.2021.103879
  * Burrows et al. (2018) Burrows, A., Vartanyan, D., Dolence, J. C., Skinner, M. A., & Radice, D. 2018, Space Sci. Rev., 214, 33, doi: 10.1007/s11214-017-0450-9
  * Cardall et al. (2013) Cardall, C. Y., Endeve, E., & Mezzacappa, A. 2013, Phys. Rev. D, 87, 103004, doi: 10.1103/PhysRevD.87.103004
  * 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
  * Char et al. (2015) Char, P., Banik, S., & Bandyopadhyay, D. 2015, ApJ, 809, 116, doi: 10.1088/0004-637X/809/2/116
  * Couch (2013) Couch, S. M. 2013, ApJ, 765, 29, doi: 10.1088/0004-637X/765/1/29
  * 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
  * 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. (2014) Fischer, T., Hempel, M., Sagert, I., Suwa, Y., & Schaffner-Bielich, J. 2014, European Physical Journal A, 50, 46, doi: 10.1140/epja/i2014-14046-5
  * Furusawa et al. (2017) Furusawa, S., Togashi, H., Nagakura, H., et al. 2017, Journal of Physics G Nuclear Physics, 44, 094001, doi: 10.1088/1361-6471/aa7f35
  * Ghosh et al. (2021) Ghosh, S., Wolfe, N., & Fröhlich, C. 2021, arXiv e-prints, arXiv:2107.13016.  https://arxiv.org/abs/2107.13016
  * Harada et al. (2020) Harada, A., Nagakura, H., Iwakami, W., et al. 2020, ApJ, 902, 150, doi: 10.3847/1538-4357/abb5a9
  * Hebeler et al. (2013) Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11, doi: 10.1088/0004-637X/773/1/11
  * Hempel et al. (2012) Hempel, M., Fischer, T., Schaffner-Bielich, J., & Liebendörfer, M. 2012, ApJ, 748, 70, doi: 10.1088/0004-637X/748/1/70
  * Hempel & Schaffner-Bielich (2010) Hempel, M., & Schaffner-Bielich, J. 2010, Nucl. Phys. A, 837, 210, doi: 10.1016/j.nuclphysa.2010.02.010
  * Horowitz (2002) Horowitz, C. J. 2002, Phys. Rev. D, 65, 043001, doi: 10.1103/PhysRevD.65.043001
  * Horowitz et al. (2017) Horowitz, C. J., Caballero, O. L., Lin, Z., O’Connor, E., & Schwenk, A. 2017, Phys. Rev. C, 95, 025801, doi: 10.1103/PhysRevC.95.025801
  * Isern et al. (1991) Isern, J., Canal, R., & Labay, J. 1991, ApJ, 372, L83, doi: 10.1086/186029
  * Janka (2012) Janka, H.-T. 2012, Annual Review of Nuclear and Particle Science, 62, 407, doi: 10.1146/annurev-nucl-102711-094901
  * Kato et al. (2020) Kato, C., Nagakura, H., Hori, Y., & Yamada, S. 2020, ApJ, 897, 43, doi: 10.3847/1538-4357/ab97b2
  * Klähn et al. (2006) Klähn, T., Blaschke, D., Typel, S., et al. 2006, Phys. Rev. C, 74, 035802, doi: 10.1103/PhysRevC.74.035802
  * Lattimer & Swesty (1991) Lattimer, J. M., & Swesty, D. F. 1991, Nucl. Phys. A, 535, 331, doi: 10.1016/0375-9474(91)90452-C
  * Liebendörfer et al. (2005) Liebendörfer, M., Rampp, M., Janka, H. T., & Mezzacappa, A. 2005, ApJ, 620, 840, doi: 10.1086/427203
  * Mabanta & Murphy (2018) Mabanta, Q. A., & Murphy, J. W. 2018, ApJ, 856, 22, doi: 10.3847/1538-4357/aaaec7
  * Marek & Janka (2009) Marek, A., & Janka, H. T. 2009, ApJ, 694, 664, doi: 10.1088/0004-637X/694/1/664
  * Marek et al. (2009) Marek, A., Janka, H. T., & Müller, E. 2009, A&A, 496, 475, doi: 10.1051/0004-6361/200810883
  * Mezzacappa et al. (2020) Mezzacappa, A., Endeve, E., Messer, O. E. B., & Bruenn, S. W. 2020, Living Reviews in Computational Astrophysics, 6, 4, doi: 10.1007/s41115-020-00010-8
  * Miller et al. (2021) Miller, M., Lamb, F. K., Dittmann, A. J., et al. 2021, NICER PSR J0740+6620 Illinois-Maryland MCMC Samples, Zenodo, doi: 10.5281/zenodo.4670689
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, NICER PSR J0030+0451 Illinois-Maryland MCMC Samples, 1.0.0, Zenodo, doi: 10.5281/zenodo.3473466
  * Mirizzi et al. (2016) Mirizzi, A., Tamborra, I., Janka, H. T., et al. 2016, Nuovo Cimento Rivista Serie, 39, 1, doi: 10.1393/ncr/i2016-10120-8
  * Morozova et al. (2018) Morozova, V., Radice, D., Burrows, A., & Vartanyan, D. 2018, ApJ, 861, 10, doi: 10.3847/1538-4357/aac5f1
  * Müller (2020) Müller, B. 2020, Living Reviews in Computational Astrophysics, 6, 3, doi: 10.1007/s41115-020-0008-5
  * 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. (2013) Müller, B., Janka, H.-T., & Marek, A. 2013, The Astrophysical Journal, 766, 43, doi: 10.1088/0004-637x/766/1/43
  * Nagakura et al. (2020) Nagakura, H., Burrows, A., Radice, D., & Vartanyan, D. 2020, MNRAS, 492, 5764, doi: 10.1093/mnras/staa261
  * Nagakura et al. (2018) Nagakura, H., Iwakami, W., Furusawa, S., et al. 2018, ApJ, 854, 136, doi: 10.3847/1538-4357/aaac29
  * Nättilä et al. (2016) Nättilä, J., Steiner, A. W., Kajava, J. J. E., Suleimanov, V. F., & Poutanen, J. 2016, A&A, 591, A25, doi: 10.1051/0004-6361/201527416
  * Neustadt et al. (2021) Neustadt, J. M. M., Kochanek, C. S., Stanek, K. Z., et al. 2021, MNRAS, 508, 516, doi: 10.1093/mnras/stab2605
  * 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
  * Olson et al. (2016) Olson, J. P., Warren, M., Meixner, M., et al. 2016, arXiv e-prints, arXiv:1612.08992.  https://arxiv.org/abs/1612.08992
  * 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
  * Richers et al. (2017) Richers, S., Ott, C. D., Abdikamalov, E., O’Connor, E., & Sullivan, C. 2017, Phys. Rev. D, 95, 063019, doi: 10.1103/PhysRevD.95.063019
  * Roberts et al. (2012) Roberts, L. F., Shen, G., Cirigliano, V., et al. 2012, Phys. Rev. Lett., 108, 061103, doi: 10.1103/PhysRevLett.108.061103
  * Schneider et al. (2019) Schneider, A. S., Constantinou, C., Muccioli, B., & Prakash, M. 2019, Phys. Rev. C, 100, 025803, doi: 10.1103/PhysRevC.100.025803
  * Schneider et al. (2017) Schneider, A. S., Roberts, L. F., & Ott, C. D. 2017, Phys. Rev. C, 96, 065802, doi: 10.1103/PhysRevC.96.065802
  * 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
  * Shen et al. (2011) Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 2011, ApJS, 197, 20, doi: 10.1088/0067-0049/197/2/20
  * Shibata et al. (2011) Shibata, M., Kiuchi, K., Sekiguchi, Y., & Suwa, Y. 2011, Progress of Theoretical Physics, 125, 1255, doi: 10.1143/PTP.125.1255
  * 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
  * Suwa et al. (2013) Suwa, Y., Takiwaki, T., Kotake, K., et al. 2013, ApJ, 764, 99, doi: 10.1088/0004-637X/764/1/99
  * 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
  * Togashi et al. (2017) Togashi, H., Nakazato, K., Takehara, Y., et al. 2017, Nucl. Phys. A, 961, 78, doi: 10.1016/j.nuclphysa.2017.02.010
  * Typel et al. (2010) Typel, S., Röpke, G., Klähn, T., Blaschke, D., & Wolter, H. H. 2010, Phys. Rev. C, 81, 015803, doi: 10.1103/PhysRevC.81.015803
  * Yakunin et al. (2010) Yakunin, K. N., Marronetti, P., Mezzacappa, A., et al. 2010, Classical and Quantum Gravity, 27, 194005, doi: 10.1088/0264-9381/27/19/194005
  * Yasin et al. (2020) Yasin, H., Schäfer, S., Arcones, A., & Schwenk, A. 2020, Phys. Rev. Lett., 124, 092701, doi: 10.1103/PhysRevLett.124.092701
  * Zhang et al. (2018) Zhang, N.-B., Li, B.-A., & Xu, J. 2018, ApJ, 859, 90, doi: 10.3847/1538-4357/aac027