Effect of the Nuclear Equation of State on Relativistic-Turbulence Induced Core-Collapse Supernovae

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.


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 Corresponding author: Luca Boccioli lbocciol@nd.edu 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) 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 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 = M = 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) 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 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 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: where y is the proton fraction and δ(y) = 1 − 2y. Hence, the symmetry energy is defined as: It is also common to encounter an alternative definition of the symmetry energy, which can be derived by expanding B (n, y) around its value for symmetric matter (i.e. B (n, y = 1/2)), yielding: 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;) 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 ∼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 ∼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 nonuniform 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 subnuclear 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.  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.
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 opensource 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 × 10 4 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).  Table 1. Parameters characterizing the equations of state utilized in this study: n0 is the nuclear saturation density, J is the symmetry energy at saturation density, L is the slope of the symmetry energy, K0 is the nuclear incompressibility, m * n is the nucleon effective mass and mn 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.

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. 2016Radice et al. , 2018Mabanta & 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: where the t − t component α 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: where v turb is the turbulent velocity and D = Xρ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 K = α K v turb Λ mix is a diffusion coefficient. The two main quantities characterizing the model are the mixing length and the Brunt-Väisälä frequency where φ = ln α reduces to the Newtonian potential in the non-relativistic limit. The quantity h = 1+ +P/ρ is the relativistic enthalpy, where is the internal energy, P is the pressure, and c s is the sound speed.
The above equations depend on two parameters: α MLT and α K . Additional diffusion coefficients analogous to D 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 α 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 α MLT , which determines the magnitude of the mixing length and therefore the scale of convection.

CALIBRATION OF THE TURBULENCE MODEL
The size of the convective cells in our model depends upon the mixing length Λ MLT , which is proportional to the mixing length parameter α MLT . The larger α 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 α MLT by comparing to realistic convection in 3D simulations. As pointed out in Boccioli et al. (2021), the calibrated value of α 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 α MLT we utilize the 20 M solarmetallicity progenitor of Farmer et al. (2016), and we employ the SFHo equation of state. We then compare our simulations using 6 different values of α 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 ∼ 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 ∼ 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: where we set R * = 30 km to exclude the convection in the PNS. From the comparison shown in Figure 3, the value of α MLT that most closely matches the 3D results appears to be in the range between 1.45 and 1.5. Specifically, for α MLT = 1.48, the integrated turbulent energy at radii r 150 km is in excellent agreement with the 3D case (see the right panel of Figure 3). In other words, at ∼ 135 ms after bounce, our simulation with α 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 ∼ 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 (∼ 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 neutrinoheated 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 α 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 α MLT : 1.45, 1.48, and 1.5. The outcome of these simulations is shown in Figure 4. The explodability with α 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 α 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. The explosion fraction for α MLT = 1.48 lies just outside the 90% confidence interval of the most recent results, while for α MLT = 1.5 it is fully consistent with observations. However, since our hydrodynamical quantities are better reproduced by α MLT = 1. 48  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 α MLT = 1.5. After calibrating the value of α MLT we simulate four progenitors with masses of 9, 15, 20 and 25 M for all of the EOSs listed in Section 2.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;Müller et al. 2012;Bruenn et al. 2016) as the integral of the binding energy E bind in regions where E bind > 0.
In the GR case, following Müller et al. (2012), one can define E bind in terms of the lapse function α, the density ρ, the specific internal thermal energy th , the pressure P and the Lorentz factor W : 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 , hence making the definition of the total energy of the fluid problematic. A way around this is to define th as the difference between the internal energy -taken from the EOS at a given density, temperature and electron fraction -and the internal energy 0 , defined as the internal energy at zero temperature and at the same density and electron fraction 1 .
From this definition, it is clear that 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: where dV 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   shows the integrated turbulent energy as a function of radius for a snapshot at ∼135 ms post bounce. Notice that, given its definition in Eq. 8, E turb (r) is zero below 30 km and then it's constant for radii larger than ∼ 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 α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.
comparison among different EOSs, which have different gain and shock radii, we define a dimensionless quantity: such that ξ gain = 0 when r = r gain and ξ gain = 1 when r = r shock . The profiles of the turbulent velocity in the gain region are plotted in Figure 7 for two different time snapshots. At ∼ 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 ∼ 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 , respectively. The shaded regions represent the 90% confidence interval. Notice that in those papers they actually report the fraction of failed SN fFSN, therefore here we assume that the fraction of explosion is 1 − fFSN 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 SROtype 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.

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 ∼ 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 α 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.
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 . 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 "SROtype". 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) 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.
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 . 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 ∼ m * n T c /ρ 2/3 c (Baym & Pethick 1991), where S c , T c and ρ 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 anal- . Turbulent velocity v 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 rgain and 1 at r 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.
ysis 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×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.

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 Yakunin et al. 2010;Müller et al. 2013;Andresen et al. 2017;Morozova et al. 2018) and trigger the socalled Lepton-number Emission Self-sustained Asymmetry (LESA) (Tamborra et al. 2014). However, its impact on the explosion is still unclear, and more highresolution 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.    (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.
The quantity chosen to represent the "vigor" of convection is the total turbulent energy in the PNS, defined as: The top panels in Figure 10 show the evolution of E pns 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.
It is also interesting to analyze each progenitor separately, and study how E pns turb correlates with the strength of the explosion. For each progenitor, as can be seen in Figure 10, the hierarchy of E pns 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.
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 pns 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 α MLT exhibit larger turbulent velocities inside the PNS. This indicates that, for large values of α 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.

EOS dependence of the neutrino signal
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 , 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 heavylepton 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. 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.
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  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,  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.
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  Table 1, with α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 νe, electron antineutrinoνe and heavy lepton neutrino νx. The densities of the electron and heavy-lepton neutrinos have been multiplied by 2 and 0.5, respectively.  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 νX has been multiplied by a factor of 2 to present them on the same scale as the other neutrinos.  Table 1, with α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. fective 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.

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, α 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  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 CC-SNe 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 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 nonparametric 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. ¡

ACKNOWLEDGMENTS
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  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. found in the other progenitors (see Figure 15). To understand the origin of this discrepancy, it's useful to define the quantity.
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.
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 54 Fe, 40 Ca and 44 Ti, 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.