General Relativistic Neutrino-Driven Turbulence in One-Dimensional Core-Collapse Supernovae
DDraft version February 16, 2021
Typeset using L A TEX twocolumn style in AASTeX63
General Relativistic Neutrino-Driven Turbulence in One-Dimensional Core-Collapse Supernovae
Luca Boccioli , Grant J. Mathews , and Evan P. O’Connor Center for Astrophysics, Department of Physics, University of Notre Dame, 225 Nieuwland Science Hall, Notre Dame, IN 46556, USA The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden
Submitted to ApJABSTRACTConvection and turbulence in core-collapse supernovae (CCSNe) are inherently three-dimensional innature. However, 3D simulations of CCSNe are computationally demanding. Thus, it is valuable tomodify simulations in spherical symmetry to incorporate 3D effects using some parametric model. Inthis paper, we report on the formulation and implementation of general relativistic (GR) neutrino-driven turbulent convection in the spherically symmetric core-collapse supernova code
GR1D . This isbased upon the recently proposed method of Supernova Turbulence in Reduced-dimensionality (
STIR )in Newtonian simulations from Couch et al. (2020). When the parameters of this model are calibratedto 3D simulations, we find that our GR formulation of
STIR requires larger turbulent eddies to achievea shock evolution similar to the original
STIR model. We also find that general relativity may alterthe correspondence between progenitor mass and successful vs. failed explosions.
Keywords:
Supernovae, Supernova dynamics, General Relativity, Gravitational Collapse, RelativisticAstrophysics INTRODUCTIONDespite the fact that core-collapse supernovae (CC-SNe) have been the focus of computational simulationsfor more than 50 years, many fundamental details of themechanism driving the explosion of massive stars remainunknown. Nevertheless, our knowledge has significantlyimproved since the first theoretical (Colgate & White1966; Arnett 1966) and observational (Sonneborn et al.1987) efforts. This has been aided by an exponentialgrowth of computing power.In the early ’80s, one-dimensional (1D) sphericallysymmetric simulations revealed the crucial role of neu-trinos in triggering the explosion through the so-calleddelayed neutrino-heating mechanism (Bethe & Wilson1985; Bruenn 1985). A few years later, the first two-dimensional (2D) axi-symmetric numerical investiga-tions of CCSNe were also performed (Miller et al. 1993;
Corresponding author: Luca [email protected]
Herant et al. 1994), while more recently the extensionto three dimensions (3D) has finally become feasible(Fryer & Warren 2002). Since then, multi-D simula-tions have continued to improve, and have shed lighton the details of the explosion mechanism (cf. reviewin Ref. M¨uller (2016)). For example, the importance ofneutrino-driven turbulent convection is now widely rec-ognized (Radice et al. 2016, 2018; Mabanta & Murphy2018). Moreover, new phenomena such as the StandingAccretion Shock Instability (SASI) (Blondin et al. 2003),the Lepton-number Emission Self-sustained Asymmetry(LESA) (Tamborra et al. 2014), Collective Neutrino Os-cillations (CNO) (Duan et al. 2010) and others, havebeen revealed in multi-D simulations. However, the im-pact of these phenomena on the explosion is still a topicof active research.After many years of theoretical and observationalwork, it is now well established that spherically symmet-ric 1D simulations only lead to self-consistent explosionsof low mass (typically around 9-11 M (cid:12) ) progenitors.These stars develop an oxygen-neon-magnesium core, a r X i v : . [ a s t r o - ph . H E ] F e b Boccioli et al. whose collapse is triggered by electron captures on Neand Mg. Contrary to these so-called Electron CaptureSupernovae (ECSNe), the more common Fe-core CCSNedo not self-consistently explode in spherically symmet-ric simulations that employ a modern treatment of thenuclear equation of state (EOS) and neutrino interac-tions, with the possible exception of a zero-metallicity,9.6 M (cid:12) progenitor (Melson et al. 2015; Betranhandy &O’Connor 2020). Nevertheless, there have been a num-ber of simulations in 2D and 3D that have led to suc-cessful explosions (M¨uller et al. 2012; Lentz et al. 2015;Janka et al. 2016; Bruenn et al. 2016; O’Connor & Couch2018a; M¨uller et al. 2019; Burrows et al. 2020). How-ever, axial symmetry has recently been shown (Couch &Ott 2015) to favor explosions by artificially enhancingneutrino-heating behind the shock through an inverseturbulence cascade that is not present in 3D. Therefore,it is generally believed that the final explanation as towhat causes the explosion must come from simulationsin three spatial dimensions.Unfortunately, despite the technological improve-ments of the last few decades, three-dimensional core-collapse supernova simulations still pose a difficult com-putational challenge, even for modern supercomputers.In addition to this, 3D simulations from different groupsonly agree for a few tens of milliseconds after bounce(Cabez´on et al. 2018). Afterward, they begin to signif-icantly deviate from each other, even when they startfrom the same initial conditions. The more accessible2D simulations, that require a relatively affordable CPUtime, are also still complicated enough to the point thatdifferent groups often find different outcomes for the ex-plosion (see Table 1 from O’Connor & Couch (2018a)),although some promising benchmarking work has al-ready been done (Pan et al. 2019).In comparison, modern 1D simulations are muchfaster to run and more consistent across different codes(O’Connor et al. 2018). Hence, when starting fromthe same initial conditions, most groups obtain simi-lar results (while older 1D simulations, with outdatedsets of neutrino interactions, lead to some inconsisten-cies [cf. Mathews et al. (2020)]). This guarantees asomewhat solid foundation and makes 1D simulationsan ideal bench-marking tool for studying how differentinput physics can affect supernova explosions.Precisely due to the low computational cost associatedwith 1D simulations, there have been several attempts toartificially drive explosions in 1D codes that would oth-erwise yield failed supernovae. The advantage of havingmodels capable of artificially driving explosions is quitestraightforward: if the models are reliable , one can effi-ciently explore the parameter space of different physical phenomena. For example, one could explore the connec-tion between the explodability of CCSNe and the pro-genitor mass (O’Connor & Ott 2011; Ugliano et al. 2012;Pejcha & Thompson 2015; Ertl et al. 2016; Sukhboldet al. 2016; Ebinger et al. 2020; Couch et al. 2020), orthe impact of different input physics, such as the nuclearEOS (Schneider et al. 2019; Olson et al. 2016) and theweak interactions (Fischer et al. 2020; Guo et al. 2020;Betranhandy & O’Connor 2020)) on the explosion itself.Furthermore, one can quantify the dependence of vari-ous nucleosynthesis processes on the progenitor mass,or calculate the diffuse relic neutrino background [c.f.review in Mathews et al. (2020)]. Moreover, 1D simula-tions can efficiently study any aspect of supernovae thatrequires a large sample of progenitors and/or evolutionto late times post bounce.A number of artificially driven explosion mechanismshave been proposed over the past few decades, basedupon different parametric models (Blinnikov & Bar-tunov 1993; Woosley & Weaver 1995; Ugliano et al. 2012;Perego et al. 2015; Sukhbold et al. 2016; Couch et al.2020). For most of them, the parameters were chosenby comparing the explosion energies and other relevantobservables to known measured properties of SN 1987a(Sonneborn et al. 1987) and other SNe.The advantage of comparing with observations is thatone can directly test the results with respect to realdata. On the other hand, since these parametric modelsare not derived from a physical mechanism, one cannotmeaningfully compare these results to 3D simulationsand gain insight on the physics of the explosion. There-fore, a different (and complementary) way of approach-ing the problem is to develop a parametric model basedupon a known physical mechanism that is motivated bythe success of 3D simulations (M¨uller et al. 2016; Couchet al. 2020).An obvious candidate for such a mechanism isneutrino-driven turbulent convection, which plays animportant role in CCSNe (Couch & Ott 2015; Radiceet al. 2016, 2018; Mabanta & Murphy 2018). Recently,Couch et al. (2020) (hereafter CWO20) developed amodel called
STIR (Simulated Turbulence In Reduced-dimensionality), based upon a modified Mixing LengthTheory (MLT), that simulates features of 3D turbulencein spherically symmetric simulations.
STIR is a Newtonian model for turbulent convection,and the simulations from CWO20 only partially includegeneral relativistic effects through a General Relativis-tic Effective Potential (GREP) from Marek et al. (2006).However, we know that General Relativity (GR) playsan important role in the explosion of supernovae (Wil-
D Turbulence in GR
STIR model to a generalrelativistic treatment. We then analyze the differencesbetween full GR and GREP. In Section 2 we describethe essential physics of STIR and its generalization toGR. Results and a comparison to 3D simulations are de-scribed in Section 3, followed by a study of explodabilityvs. progenitor mass in Section 4. Conclusions are givenin Section 5. Throughout the manuscript, we chose toadopt natural units, i.e. G = c = M (cid:12) = 1. METHODS2.1.
The STIR model of Couch et al. (2020)
The importance of convection in CCSNe was recog-nized very early on, and the first attempts at applyingMLT to CCSN date back to Wilson & Mayle (1988)and Bruenn et al. (1995). We follow a more modernversion of those models, that can be found in Mabanta& Murphy (2018), Mabanta et al. (2019) and CWO20.We highlight that this formulation does consider bothpressure and composition gradients (see Appendix A ofM¨uller et al. (2016) for an explicit proof of this). How-ever, it doesn’t capture lepton number-driven convec-tion, relevant at the surface of the PNS, where neutrinosare trapped below the neutrinosphere and free streamingabove it, creating a gradient which could drive convec-tion. However, given the complicated thermodynamicsderivatives involved in the modeling of this type of con-vection (Wilson & Mathews 2003), we leave its treat-ment to a future work.As described in the aforementioned papers, the ef-fects of turbulence can be treated as a perturbation ona background fluid. This can be done by decompos-ing a generic fluid variable ψ into a mean component ψ and a perturbed component ψ (cid:48) : ψ = ψ + ψ (cid:48) , with (cid:104) ψ (cid:105) = ψ , where (cid:104) ... (cid:105) represents an average over a properspatial and time domain. Because of the chaotic natureof the motions that characterize turbulence, the averageshould vanish (i.e. (cid:104) ψ (cid:48) (cid:105) = 0). This decomposition, calledReynolds averaging, can be applied to the compressibleEuler equations to derive conservation laws that con-tain explicit contributions from the turbulent motion fortotal energy, momentum and mass. In the Newtonianlimit, these are: ∂ t (cid:104) ρ (cid:105) + ∇ · ( ρ v ) = 0 , (1) ∂ t (cid:104) ρ v (cid:105) + ∇ · ( ρ v × v + P I ) = − ρ g − ∇ · (cid:104) ρ R (cid:105) , (2) ∂ t (cid:104) ρ(cid:15) (cid:105) + ∇ · u ( ρ (cid:15) + P )= − ρ v · g − ∇ · v (cid:104) ρ R (cid:105) − ∇ · F e + ρ (cid:15) (cid:48) . (3)Here, ρ is the mass density, v is the velocity, P is thepressure, I is the identity matrix, g is the gravitationalacceleration, (cid:15) is the internal energy per unit mass, R isthe Reynolds stress tensor and F e is the energy flux dueto turbulence. In the above equations, we have neglectedterms that are proportional to the turbulent pressure.This is a good approximation in regimes of low Machnumbers (for more details on the derivation of equa-tions (1)-(3) see Murphy & Meakin (2011), Mabanta &Murphy (2018), and CWO20).In addition to Eqs. (1)-(3), one needs another equa-tion to describe the evolution of the specific turbulentkinetic energy K = Tr( R ) /
2. The derivation of thisequation can be found in Murphy & Meakin (2011) andCWO20. Here, we simply highlight the important stepsthat are relevant to our model.For isotropic turbulence in spherical symmetry, theevolution equation for the turbulent kinetic energy is: ∂ (cid:104) ρK (cid:105) ∂t + 1 r ∂∂r [ r ( (cid:104) ρK (cid:105) v r + (cid:104) ρKv (cid:48) r (cid:105) )]= − (cid:104) ρ R rr (cid:105) ∂v r ∂r + (cid:104) ρ (cid:48) v (cid:48) r (cid:105) g − ρ(cid:15) (cid:48) . (4)One can rewrite almost all of the unknown expressionsin (4) in terms of gradients and divergences of v turb ≡ v (cid:48) and (cid:15) turb ≡ (cid:15) (cid:48) , where: (cid:15) turb = v Λ (5)with Λ being the largest turbulent eddy size. The onlyquantities that remain to be determined are the buoy-ancy force (cid:104) ρ (cid:48) v (cid:48) r (cid:105) g and Λ. Therefore, to close the systemof equations, one needs an expression that relates thetransport of turbulent energy to the typical speed of abuoyant convective element. To do this, one refers toMLT, for which the buoyant force that makes the con-vective element rise through the stratified fluid is: (cid:104) ρ (cid:48) v (cid:48) (cid:105) g ≈ ρv turb ω bv Λ mix (6)with: Boccioli et al. Λ ≡ Λ mix = α mlt Pρg , (7) ω bv = g eff (cid:18) ρ ∂ρ∂r − ρc s ∂P∂r (cid:19) . (8)In the above equations, Λ mix is the mixing length, ω bv is the Brunt-V¨ais¨al¨a frequency, c s is the sound speed,and g eff is the magnitude of the local effective accelera-tion. For a fluid in hydrostatic equilibrium, g eff simplyreduces to the local gravitational acceleration g . Morein general, however, in the rest frame one should takethe acceleration of the fluid into account. Therefore, thetotal acceleration g eff can be expressed as: g eff = g − v r ∂v r ∂r , (9)as described in CWO20, where v r is the radial velocity ofthe fluid. During the pre-bounce and early post-bouncephases, both v r and ∂v r /∂r are small, so that g eff simplyreduces to g . Once the explosion sets in, however, thecontribution from the acceleration of the fluid is non-negligible, hence g eff becomes smaller than g . This ef-fectively shuts off turbulent convection in the explodingmaterial by driving ω bv to zero.The mixing length Λ mix is the average distance that aconvective element will travel before being mixed with(and increasing the internal energy of) the surroundingmaterial. The Brunt-V¨ais¨al¨a frequency ω bv is the rateat which the convective elements are rising. As one cannotice from Eq. (8), ω bv can be either positive or nega-tive: when ω bv > ω bv < α mlt , which scalesthe mixing length to the pressure scale height in Eq. (7).Typically α mlt ∼ O (1).After closing the system of equations derived from theReynolds decomposition, one obtains the following evo-lution equation for the turbulent kinetic energy: ∂ρv ∂t + 1 r ∂∂r [ r ( ρv v r − ρD K ∇ v )]= − ρv ∂v r ∂r + ρv turb ω bv Λ mix − ρ v Λ mix . (10)In the above equation, D K is a diffusion coefficient dueto turbulence, defined as: D K = α K v turb Λ mix . (11) Similar terms appear in the internal energy, electronfraction, and neutrino energy density evolution equa-tions (for the complete set of hydrodynamic equationsused in the model, see Eqs. 25-29 and 33 in CWO20).Therefore, strictly speaking, STIR has 4 additional freeparameters: α K , α e , α Ye , α ν . However, the convectivemotions are not very sensitive to the value of these pa-rameters, so we set them to 1 / STIR in General relativity
The first attempts to create a general relativisticmodel for convection date back to Thorne (1966). Wewill follow the same approach, but using a slightly dif-ferent formalism.All the simulations described in this paper were runwith the open-source, spherically symmetric, general rel-ativistic code
GR1D (O’Connor & Ott 2010). The Boltz-mann equation for neutrino transport is solved usingan M1-scheme, with opacity tables generated using theopen-source code
NuLib (O’Connor 2015).The metric evolved in
GR1D is Schwarzchild-like: ds = g µν x µ x ν = − α ( r, t ) dt + X ( r, t ) dr + r d Ω , (12)where α and X can be expressed as functions of a metricpotential φ (which reduces to the Newtonian potentialin the Newtonian limit) and the enclosed gravitationalmass M grav : α ( r, t ) = exp[ φ ( r, t )] ,X ( r, t ) = (cid:18) − M grav ( r, t ) r (cid:19) − / (13)For the present work, we first note that turbulence ismostly relevant far from the proto-neutron star (PNS)where GR effects can be treated as a perturbation.Therefore, one can simply make a few changes to theterms in Eq. (10) without having to re-derive the en-tire Reynold’s decomposition. The expression for ω bv ,however, must be carefully re-derived (a formal deriva-tion of ω bv for a fluid in hydrostatic equilibrium can befound in Appendix A). Far from the PNS, we invoke thefollowing:(i) replace the conserved variable ρ with its GR coun-terpart, i.e. D = W Xρ , where W = (1 − v ) − / and v = Xv r ;(ii) multiply the RHS of Eq. (10) by αX . D Turbulence in GR α/X (seeO’Connor & Ott (2010) for more details on thederivation of the hydrodynamic equations in
GR1D .We note, however, that GR effects are not small withinthe PNS, where interior lepton number gradients candrive convection and affect the neutrino emission (Wil-son & Mathews 2003). We leave a proper treatment oflepton number-driven convection to a future work.The expression of ω bv , with the assumption of thebackground fluid being in hydrostatic equilibrium, is: ω bv = α ρhX d φ d r (cid:18) ∂ρ (1 + (cid:15) ) ∂r − ρc s ∂P∂r (cid:19) . (14)For an accelerating fluid, the GR derivation becomesmuch less obvious, since non-inertial frames are involved.Fortunately, GR effects are very small in the outer lay-ers, where the fluid has a non-negligible acceleration.Hence, we can simply use the Newtonian expression (9)with ad hoc corrections for relativistic space-time cur-vature, so that ω bv becomes: ω bv = α X (cid:18) d φ d r − v ∂v∂r (cid:19) (cid:18) ρh ∂ρ (1 + (cid:15) ) ∂r − ρhc s ∂P∂r (cid:19) , (15)where now v = Xv r .The main difference between Eqs. (8) and (15) is theinclusion of ∂ρ(cid:15)/∂r in the latter. In the gain region theinternal energy decreases with radius, i.e. ∂ρ(cid:15)/∂r < ω bv and therefore theamount of turbulence that is generated. We will comeback to this in Section 3.22.3. A word of caution: conservation of energy
As briefly noted in CWO20, the equations of the
STIR model do not conserve energy. As pointed out by M¨uller(2019), the first and second term on the RHS of Eq.(4) break total energy conservation. Since turbulenceis driven by neutrinos interacting with the material inthe gain region, the rate of turbulent energy generationby buoyant driving is of the same order as the volume-integrated neutrino heating rate ˙ Q ν . This is also what(Mabanta et al. 2019) assume in their turbulence modelin the gain region, where the buoyant driving is givenby β ˙ Q ν .The extra energy coming from buoyant driving breaksthe energy-conserving structure of the hydrodynamicalequations. As pointed out in Mabanta & Murphy (2018)this can be thought of as energy extracted from the freeenergy stored in unstable compositional and thermody-namic gradients. Therefore, although the energy is com-ing from a physically motivated phenomenon, it cannot be consistently modeled in spherical symmetry, as rec-ognized by M¨uller (2019).If one looks at the problem of energy conservation inother artificially-driven explosions in one dimension, it isclear that energy is not conserved in any of them. Thisreflects the main issue with 1D models of CCSNe: theydo not self-consistently explode. Hence, generally speak-ing, to drive an explosion in one dimension one needsmore efficient neutrino heating. This extra-heating isinjected in the model following different prescriptions,that can be based on physical phenomena (Mabanta &Murphy 2018; Couch et al. 2020) and/or observationalconstraints (Ugliano et al. 2012; Sukhbold et al. 2016),or simply by well-reasoned tweaking of some parame-ters (O’Connor & Ott 2011; Perego et al. 2015). Un-til multi-dimensional models become computationallymore affordable, the above methods can give differentand complementary insights into the physics behind theexplosion mechanism of CCSNe. RESULTS: COMPARISON WITH 3DSIMULATIONS3.1.
Results using an effective potential
Following the approach of CWO20, we compare ourGREP model to the mesa20 LR v
3D simulation byO’Connor & Couch (2018b). In the original
STIR paper,the advantage of carrying out such a comparison is thatboth O’Connor & Couch (2018b) and CWO20 use thecode
FLASH , so they can directly compare their 1D and3D simulations with everything else unchanged. Despiteusing a different code, we still expect our GREP modelto show very little differences to the original
STIR , since
GR1D and
FLASH have been shown to yield similar resultsin 1D (O’Connor et al. 2018).In all the simulations presented in this section weused the same setup chosen by CWO20. That is, wesimulated the collapse of a 20 M (cid:12) progenitor fromFarmer et al. (2016), adopting the SFHo EOS (Steineret al. 2013) and assuming Nuclear Statistical Equilib-rium (NSE) everywhere. The neutrino radiation trans-port solver used in
FLASH closely resembles the one usedin
GR1D (O’Connor 2015), and we used the same set of
NuLib opacities adopted in CWO20. For the neutrinos,we used 12 energy groups geometrically spaced up to ∼
250 MeV.In The upper panels of Figure 1 we plot the shockradius versus time and the turbulent velocity profile at ∼
135 ms post-bounce for our GREP model (to be com-pared to Figures 1 and 2 of CWO20).
GR1D consistentlyyields larger values than
FLASH for the turbulent veloc-ity at a given α mlt . This translates into larger values ofthe shock radius at a given time. Except for these small Boccioli et al. differences, the two models agree very well, yielding ex-plosions for α mlt > .
2. It is worth noting, however, thatthe 3D results show features that are not present in ourMLT-like model. Specifically, the convective speed in3D has a slower-decaying tail at 50-80 km. As pointedout by CWO20, this is due to angular variations in the3D model, rather than presence of convection in the re-gion below the gain layer. As one would expect, in fact,convection in our model shuts off at 80 km, which isapproximately the location of the gain layer. More in-teresting is the lack of PNS convection at 25 km, notcaptured by our model. A possible explanation for thisis that we are not taking lepton number-driven convec-tion into account, and with a more careful treatment ofthis type of convection, this feature might be correctlyreproduced, but this goes beyond the scope of this pa-per. 3.2.
Results using GR
We approached the analysis of the simulations in fullGR from two different points of view: on the one hand,we can compare the GR and GREP models with every-thing else being unchanged; on the other hand, we cancompare the GR model to a 3D simulation in full GR.For the first part of the analysis, we show the resultsof our GREP model in the upper panels of Figure 1,while the bottom panels refer to the runs with full GR.The first and most important difference to notice is therange of α mlt used in the GREP and GR simulations. Toachieve turbulent velocities and shock radii that are sim-ilar to the GREP results, the model in full GR requiresan α mlt that is ∼
20% larger. To understand the originof this increase in α mlt , it’s useful to look at the form of ω bv . As already pointed out in Section 2.2, the inclusionof the internal energy gradient in eq. (15), compared toeq. (8), is in fact the main difference between the GREPand GR models.In the gain region, which is where turbulence is mostrelevant, the internal energy gradient is negative, andtherefore it decreases ω bv , making the fluid more stableagainst convection. The inclusion of ρ(cid:15) in the definitionof ω bv is therefore needed for a realistic characteriza-tion of turbulent convection. Indeed, if one takes theform of ω bv from GR and implements it in the GREPmodel, the value of α mlt needed to achieve an explosionincreases, becoming comparable to the one we use forthe GR model.The second part of the analysis was carried out bycomparing the full GR model to the 3D octant simu-lations of Schneider et al. (2019). For this part of theanalysis, we simulated a 20 M (cid:12) progenitor from Woosley& Heger (2007) with two of the EOS’s used in Schnei- der et al. (2019): SLy4 . (Chabanat et al. 1997) andLS220 (a version of the original EOS with incompress-ibility K ∞ = 220 MeV from Lattimer & Swesty (1991)recomputed using the SRO code from da Silva Schneideret al. (2017)). For the neutrino opacities, we adopted thesame NuLib tables used by Schneider et al. (2019), with24 energy groups logarithmically spaced up to 280 MeV.Like Schneider et al. (2019), at 20 ms after bounce weswitch from inelastic electron-neutrino scattering to asimpler elastic treatment, and we also shut off velocity-dependent terms in the transport equation. Unlikethem, we consistently keep 24 energy groups through-out the whole simulation, while they switch to 16 en-ergy groups at 20 ms post-bounce. Our choice ensuresa smooth transition in our code from inelastic to elasticelectron-neutrino scattering.A comparison between the 3D data and our 1D modelis presented in Figure 2. As one can immediately see, forthe LS220 case, both the shock radius and the turbulentvelocity profile show very good agreement between the3D and our 1D results. As already pointed out in section3.1, minor differences in the turbulent velocity profile atsmall radii can be explained as the consequence of angu-lar variations in the 3D simulation. The SLy4 case is abit more complicated, and the agreements is not as goodas in the LS220 case. After ∼
300 ms the shock trajec-tory in 1D starts deviating from the 3D case, suggestingthat turbulence in our MLT-like model might be slightlymore efficient than in 3D, hence generating a successfulexplosion. A similar feature can be seen in the LS220case, in which after ∼ α mlt ,the turbulence inside the PNS is weaker. This is consis-tent with the findings of CWO20. In their paper, they(semi-quantitatively) assert that the convection insidethe PNS is too fast to allow a build-up of turbulenceenergy. This is because the gradients are quickly flat-tened so there isn’t enough time for turbulence to de-velop. This is the reason why for small values of α mlt our model is able to generate some turbulence outsidethe PNS, but as convection becomes more efficient (i.e.for larger values of α mlt ) the gradients are flattened and ω bv becomes smaller, resulting in less turbulent energygenerated above the PNS.3.3. Another word of caution: Resolution
D Turbulence in GR Time post-bounce [s] S h o c k R a d i u s [ k m ] GREP α MLT = 1.03D FLASH 0.00.20.40.60.81.01.21.4 α M L T (a) Radius [k ] C o n v e c t i v e s p ee d [ k s − ] GREP α MLT = 1.03D FLASH 0.00.20.40.60.81.01.21.4 α M L T (b) Time post-bounce [s] S h o c k R a d i u s [ k m ] GR α MLT = 1.03D FLASH 0.00.20.40.60.81.01.21.41.61.8 α M L T (c) Radius [km] C n v e c t i v e s p ee d [ k m s − ] GR α MLT = 1.03D FLASH 0.00.20.40.60.81.01.21.41.61.8 α M L T (d) Figure 1.
The plots on the upper row were generated using our GREP model, while the plots on the bottom row were generatedusing full GR. Panels (a,c) show the time evolution of the shock radius for different values of the parameter α mlt , and can becompared to Figure 2 from Couch et al. (2020). Panels (b,d) show a snapshot at ∼
135 ms post bounce of v turb , and can becompared to Figure 1 from Couch et al. (2020). The dashed lines are from the 3D simulation of O’Connor & Couch (2018b). It is well known that, in 3D, low resolution simulationscan overestimate the amount of turbulence produced inthe gain region. The turbulence cascade is in fact arti-ficially truncated at larger scales than what would oth-erwise happen in higher resolution simulations (Radiceet al. 2016). Low resolution leads to an accumulationof turbulent energy at the largest scales, resulting in anadditional pressure support in the gain region, whichfacilitates an explosion. In an MLT approach, this cantypically be avoided, but other problems may arise (fora more detailed discussion see Appendix B). For the remainder of the paper, we chose to performsimulations at a higher spatial resolution and with moreenergy groups with respect to the results presented inSection 3. As explained in Appendix B, a higher resolu-tion in space and in neutrino energy groups can changehow much turbulent energy is injected in the model.However, it is quite unlikely that this will affect the qualitative results (i.e. explosion or fallback) of our sim-ulation, because we are using a parametric model whereone can choose how much turbulent energy will be in-jected in the model simply by changing the value of α mlt .Therefore, if a simulation run with higher spatial resolu- Boccioli et al.
Time post-bounce [s] S h o c k R a d i u s [ k m ] LS220 α MLT = 13D 0.00.20.40.60.81.01.21.4 α M L T (a) Radius [km] C o n v e c t i v e s p ee d [ k m s − ] LS220 α MLT = 13D 0.00.20.40.60.81.01.21.4 α M L T (b) Time post-bounce [s] S h o c k R a d i u s [ k m ] SLy4 α MLT = 13D 0.00.20.40.60.81.01.21.4 α M L T (c) Radius [km] C o n v e c t i v e s p ee d [ k m s − ] SLy4 α MLT = 13D 0.00.20.40.60.81.01.21.4 α M L T (d) Figure 2.
Comparison of our GR turbulent convection model (solid lines) as a function of α mlt with the 3D octant simulationsfrom Schneider et al. (2019) (dashed lines). The upper (lower) panels refer to the simulations using the LS220 (SLy4) EOS. Thepanels on the left show the shock evolution, while the panels on the right show the turbulent velocity profile at ∼
220 ms afterbounce. tion develops more (or less) turbulence, one can obtainthe same qualitative result in a lower resolution simula-tion by increasing (or decreasing) the value of α mlt .For the remainder of the paper, however, we choseto perform all of the simulations with a higher spatialresolution and with more neutrino energy groups withrespect to Section 3. PROGENITOR STUDYAfter validating our turbulent convection model bycomparing it to the 3D results of Couch & O’Connor(2014) and Schneider et al. (2019), we simulated the collapse and subsequent shock revival of 20 progenitorsfrom Sukhbold et al. (2016) for different values of α mlt ,chosen such that the fraction of successful explosions isapproximately between 25% and 80%. All of the simu-lations were run using a high spatial resolution, shownin the upper panel of Figure 4 as a dashed line, with 18neutrino energy groups geometrically spaced from 0 to280 MeV.From the analysis carried out in section 3, it is clearthat once an explosion sets in for a certain value α ∗ mlt ,that progenitor will explode for all α mlt > α ∗ mlt (see Fig-ures 1 and 2). Therefore, if a given progenitor explodes D Turbulence in GR α mlt , it will also explode for highervalues of α mlt . Analogously, if a given progenitor endsup in a failed explosion for a given value of α mlt , it willalso fail for lower values of α mlt . Hence, we did notrun simulations for those values of α mlt that could bealready predicted by the outcome of the simulations atother values of α mlt .The results of our GREP model agree with the onesobtained by CWO20, although shifted by ∆ α mlt (cid:39) . STIR results for a higher value of α mlt ,as Figure 1 might suggest, but it modifies the explosionpattern of CCSNe. To accurately characterize the dif-ferences between the explodability patterns of the GRand GREP models, a study with more progenitors anda finer resolution in α mlt would be desired. That, how-ever, is beyond the scope of the present work.We can conclude from Figure 3 that general relativ-ity changes how likelihood for certain progenitors to ex-plode. Focusing on the α mlt = 1 .
27 and α mlt = 1 . (cid:12) progenitors, andfailed explosions for the 18 M (cid:12) . In the GR model, theexplodability is the exact opposite.The explodability pattern of the GR model with α mlt = 1 .
48 is intermediate between the results ofCWO20 and Sukhbold et al. (2016): (i) like the former(and unlike the latter), we find that low mass progen-itors with M = 13 -15 M (cid:12) result in failed explosions;(ii) like the latter (and unlike the former) we find thathigher mass progenitors with M = 24 -25 M (cid:12) result infailed explosions.It is worth mentioning that the explodability patternfound in the GR model with α mlt = 1 .
48 cannot bereproduced using the GREP model. Focusing on theupper panel of Figure 3, one can notice how, in the α mlt = 1 .
25 case, the 25 M (cid:12) progenitor fails, but atthe same time all progenitors below 22 M (cid:12) fail as well.This shows that the GREP model cannot produce suc-cessful explosions for low mass progenitors and at thesame time failed explosions of the 24-25 M (cid:12) progenitors.Notably, the GR model with α mlt = 1 . α mlt = 1 .
27 (with the only exception of the 18 M (cid:12) progenitor). This tells us that: (i) the threshold between successful and failed explosions is a very steep functionof α mlt ; and (ii) GR can indeed reproduce the GREP re-sults. However, GR also shows another possible patternnot yet observed in the GREP model of CWO20.Multi-D simulations cannot yet produce results anal-ogous to Figure 3, but the number of 3D simulationsin the last 5 years has dramatically increased. In par-ticular, a (relatively) large set of progenitors has beenrecently simulated by Burrows et al. (2020). In theirwork, they simulate 14 progenitors from Sukhbold et al.(2016), 12 of which are in the range M = 9 −
20 M (cid:12) ,and the remaining two have M = 25 M (cid:12) and M = 60M (cid:12) . Moreover, they use the same SFHo EOS that weemployed for our progenitor study. However, like thevast majority of 3D simulations, they only include gen-eral relativistic effects in the hydrodynamics through theGREP model of Marek et al. (2006), which they also useto account for redshift effects in the neutrino transport.Hence, we can qualitatively compare the outcome of oursimulations with theirs.In Burrows et al. (2020) all of the progenitors explodeexcept for the 13, 14 and 15 M (cid:12) , which is exactly whatwe find in the GREP model with α mlt = 1 .
27. The GRmodel with α mlt = 1 .
48 finds that the 25 M (cid:12) shouldalso fail, and it is interesting to note that in Burrowset al. (2020) that progenitor is the one with the explo-sion occurring at the latest time, i.e. it is triggered at ∼
400 ms post bounce. Therefore, it can be considereda weak explosion, and a slightly less effective neutrinoheating would cause the supernova to fail. This is con-sistent with what the GR model is showing. That is,among all progenitors in the range 9-25 M (cid:12) , the 25 M (cid:12) progenitor is the one that requires the largest value of α mlt to explode (i.e. ∼ . CONCLUSIONSIn this article, we have shown how to extend
STIR , theMLT-like model of CWO20, to a full general relativisticformalism. After verifying that we can reproduce theresults of CWO20 using the same GREP model thatthey developed, we showed how a GR version of this0
Boccioli et al. = succesfull explosio = failed explosio
ZAMS mass [M ⊙ ]1.23 α M L T GR⊙P (a) = s ccesf ll explosion= failed explosion
ZAMS mass [M ⊙ ]1.45 α M L T ⊙R (b) Figure 3.
Explosion pattern of CCSN for the GREP (upperpanel) and GR (lower panel) models as a function of the ZeroAge Main Sequence mass. Green bands represent successfulexplosions (i.e. the shock has reached 500 km), while blackbands represent failed explosions. model needs a larger value of α mlt to reproduce similar shock dynamics. The reason behind this is the inclusionof the internal energy gradient in eq. (15), which reducesthe magnitude of ω bv , hence inhibiting the generation ofturbulence, so that larger values of Λ mix are needed todevelop convective mixing as strong as in a Newtonianformalism.The GR model reproduces very well the shock dynam-ics and turbulent energy generation in the octant 3Dsimulations of Schneider et al. (2019), performed withtwo different equations of state.After comparing our results to 3D results, we ran sim-ulations of 20 different progenitors from (Sukhbold et al.2016) for different values of α mlt , for both the GREPand GR models. We find that GR changes the explod-ability pattern of CCSNe. Specifically, the 24-25 M (cid:12) progenitors seem to be hared to explode with the GRmodel. This results in an explodability as a function ofprogenitor mass that is intermediate between the resultsof CWO20 and (Sukhbold et al. 2016). However, the GRmodel also shows, for α mlt = 1 .
5, an explodability thatis compatible with the results of CWO20.ACKNOWLEDGMENTSThe authors would like to thank Sean Couch, Andreda Silva Schneider and Mackenzie Warren for fruitfuldiscussions. Work at the University of Notre Dame sup-ported by the U.S. Department of Energy under NuclearTheory Grant DE-FG02-95-ER40934. EOC would liketo acknowledge Vetenskapsr˚adet (the Swedish ResearchCouncil) for supporting this work under award numbers2018-04575 and 2020-00452.APPENDIX A. DERIVATION OF THE BRUNT-V ¨AIS ¨AL ¨A FREQUENCYA derivation of ω bv in general relativistic hydrostatic equilibrium can be found in Appendix C of M¨uller et al. (2013).For clarity, here we go over a similar derivation, but for a different gauge, trying not to skip any non-obvious steps.In the fluid reference frame, the velocity is simply the convective velocity v c . Hence, the 4-velocity is: u µ = ( W/α, W v c , ,
0) (A1)Conservation of momentum requires that ∇ µ T µr = 0, where: T µν = ( ρh + δρh ) u µ u ν + g µν P (A2) T tt = ( ρh + δρh ) W − Pα , T tr = ( ρh + δρh ) W v c αX , T rr = ( ρh + δρh ) W v c + PX (A3) T θθ = Pr , T φφ = Pr sin θ (A4) D Turbulence in GR ∇ µ T µr = 0 = 1 √− g ∂ µ ( √− gT µr ) + Γ rρσ T ρσ = 1 √− g ∂ t ( √− gT tr ) + 1 √− g ∂ r ( √− gT rr ) + Γ rtt T tt + Γ rrr T rr + Γ rθθ T θθ + Γ rφφ T φφ + 2Γ rrt T rt = 1 αX ∂∂t (cid:2) ( ρh + δρh ) W v c (cid:3) + 1 αXr ∂∂r (cid:18) αr X (cid:2) P + ( ρh + δρh ) W v c (cid:3)(cid:19) + 1 X [( ρh + δρh ) W − P ] ∂ r φ + PX ∂ r X − rX Pr − r sin θX Pr sin θ + 2( ρh + δρh ) W v c αX ∂ t X (A5)The Christoffel symbols for the metric (12) can be found in Table A.1 from O’Connor & Ott (2010).By explicitly evaluating the derivative in the second term of Eq. (A5), it can be easily shown that:1 αXr ∂∂r (cid:18) αr X P (cid:19) = 1 X ∂P∂r + PX ∂φ∂r + 2 PrX − PX ∂ r X (A6)In addition to this, using Eqs. A.5-7 from O’Connor & Ott (2010) one has: v c ∂ t X = v c X r ∂ t m ∝ (cid:19)(cid:19)(cid:55) v c (A7)By combining Eqs. (A5) and (A6), and using ρh t = ρh + δρh we arrive at:1 X ∂P∂r + 1 αX ∂∂t [ ρh t W v c ] + 1 αXr ∂∂r αr X ρh t W (cid:19)(cid:19)(cid:55) v c = − ρh t W X ∂φ∂r . (A8)Here we have ignored terms of order v c . Now one can use the condition of hydrostatic equilibrium: ∂P∂r = − ρh ∂φ∂r . (A9)By considering subsonic motions in the gain region (hence v c << c s << c ⇒ W = 1), one obtains:1 αX ∂∂t ( ρh t v c ) = − δρhX ∂φ∂r = − X ∂φ∂r C L δr , (A10)where for C L we have used the Schwarzschild discriminant from Thorne (1966): C L = (cid:18) ∂ρ (1 + (cid:15) ) ∂r − c s ∂P∂r (cid:19) . (A11)By using conservation of energy and rest mass (i.e. equations A.13 and A.17 from O’Connor & Ott (2010)), andassuming pressure equilibrium, it can be shown that: v c ∂ρh t ∂t ∝ (cid:19)(cid:19)(cid:55) v c (A12)Therefore Eq. (A10) becomes: ∂v rc ∂t = − αρh t X ∂φ∂r C L δr . (A13)Since ¨ r = α ∂v rc ∂t (where the factor of α converts t to proper time), we finally have:¨ r = − α X ρh t ∂φ∂r C L δr = ⇒ ω = α ρh t X ∂φ∂r (cid:18) ∂ρ (1 + (cid:15) ) ∂r − c s ∂P∂r (cid:19) (A14)2 Boccioli et al.
Time post-bounce [s] S h o c k R a d i u s [ k m ] GRGREP 10 −1 Radius [km] d R [ k m ]
50 100 150 2000.00.51.01.52.0 (a)
Time post-bounce [s] S h o c k R a d i u s [ k m ] GRGREP 0 50 100 150 200 250
Energy [MeV] d E [ M e V ] (b) Figure 4.
Evolution of the shock radius for the mesa20 LR v progenitor by O’Connor & Couch (2018b) for the GR (blue lines)and GREP models (red lines). The GR models were all run for α mlt = 1 .
5, while the GREP models were all run with α mlt = 1 . D Turbulence in GR B. DEPENDENCE OF MLT-LIKE TURBULENCE ON SPATIAL AND ENERGY RESOLUTIONBecause turbulence is most effective behind the shock, it is crucial to have high resolution in the proximity of theshock, especially in the first few hundred milliseconds after bounce. In fact, in this time window, the shock is stallingat ∼
150 km, and turbulent energy can provide the extra-pressure necessary to trigger the explosion. In our numericaltreatment of equation (10) we use a 5-point interpolation formula to evaluate the gradient of v r , and to avoid numericalinstabilities we shut off turbulence within two zones from the shock. This ensures that the turbulence is not advectedthrough the shock, which, being a discontinuity in pressure and radial velocity, would make the evaluation of gradientsproblematic. However, this means that in our model there is no turbulence in the two zones right behind the shock.Therefore, if the resolution behind the shock is poor (typically ∆ r ∼ / . r ∼ /
150 km), one ends upshutting off turbulence in the most important regions (i.e. close to the shock).We ran some tests by changing the spatial resolution for the GREP and GR model (see upper panel of Fig. 4) forthe same mesa20 LR v progenitor used in Section 3. From these tests it appears that, with the increased resolution inthe GR model, the shock is pushed to larger radii at early times. This allows a widening of the gain region, generatinga stronger explosion. This feature isn’t present in the GREP model, which on the contrary shows weaker explosionsat higher spatial resolutions. We also ran some tests by changing the resolution in energy. For this case the GR andGREP models showed the same behaviour, i.e. higher energy resolution yields weaker explosions (see lower panel ofFig. 4).The dependence of the shock dynamics on resolution is not as straightforward as what might be expected based onsimple arguments. However, the effects are relatively small, limited to a few percents in α mlt , and therefore do notsignificantly affect our analysis. REFERENCES Arnett, W. D. 1966, Canadian Journal of Physics, 44, 2553,doi: 10.1139/p66-210Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14,doi: 10.1086/163343Betranhandy, A., & O’Connor, E. 2020, PhRvD, 102,123015, doi: 10.1103/PhysRevD.102.123015Blinnikov, S. I., & Bartunov, O. S. 1993, A&A, 273, 106Blondin, J. M., Mezzacappa, A., & DeMarino, C. 2003,ApJ, 584, 971, doi: 10.1086/345812Bruenn, S. W. 1985, ApJS, 58, 771, doi: 10.1086/191056Bruenn, S. W., Mezzacappa, A., & Dineva, T. 1995, PhR,256, 69, doi: 10.1016/0370-1573(94)00102-9Bruenn, S. W., Lentz, E. J., Hix, W. R., et al. 2016, ApJ,818, 123, doi: 10.3847/0004-637X/818/2/123Burrows, A., Radice, D., Vartanyan, D., et al. 2020,MNRAS, 491, 2715, doi: 10.1093/mnras/stz3223Cabez´on, R. M., Pan, K.-C., Liebend¨orfer, M., et al. 2018,A&A, 619, A118, doi: 10.1051/0004-6361/201833705 Chabanat, E., Bonche, P., Haensel, P., Meyer, J., &Schaeffer, R. 1997, NuPhA, 627, 710,doi: 10.1016/S0375-9474(97)00596-4Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626,doi: 10.1086/148549Couch, S. M., & O’Connor, E. P. 2014, ApJ, 785, 123,doi: 10.1088/0004-637X/785/2/123Couch, S. M., & Ott, C. D. 2015, ApJ, 799, 5,doi: 10.1088/0004-637X/799/1/5Couch, S. M., Warren, M. L., & O’Connor, E. P. 2020, TheAstrophysical Journal, 890, 127,doi: 10.3847/1538-4357/ab609eda Silva Schneider, A., Roberts, L., & Ott, C. 2017,PhRvC, 96, doi: 10.1103/PhysRevC.96.065802Duan, H., Fuller, G. M., & Qian, Y.-Z. 2010, AnnualReview of Nuclear and Particle Science, 60, 569,doi: 10.1146/annurev.nucl.012809.104524 Boccioli et al.
Ebinger, K., Curtis, S., Ghosh, S., et al. 2020, ApJ, 888, 91,doi: 10.3847/1538-4357/ab5dcbErtl, T., Janka, H. T., Woosley, S. E., Sukhbold, T., &Ugliano, M. 2016, ApJ, 818, 124,doi: 10.3847/0004-637X/818/2/124Farmer, R., Fields, C. E., Petermann, I., et al. 2016, TheAstrophysical Journal Supplement Series, 227, 22,doi: 10.3847/1538-4365/227/2/22Fischer, T., Typel, S., R¨opke, G., Bastian, N.-U. F., &Mart´ınez-Pinedo, G. 2020, PhRvC, 102, 055807,doi: 10.1103/PhysRevC.102.055807Fryer, C. L., & Warren, M. S. 2002, ApJL, 574, L65,doi: 10.1086/342258Guo, G., Mart´ınez-Pinedo, G., Lohs, A., & Fischer, T. 2020,PhRvD, 102, 023037, doi: 10.1103/PhysRevD.102.023037Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate,S. A. 1994, ApJ, 435, 339, doi: 10.1086/174817Janka, H.-T., Melson, T., & Summa, A. 2016, AnnualReview of Nuclear and Particle Science, 66, 341,doi: 10.1146/annurev-nucl-102115-044747Lattimer, J. M., & Swesty, D. F. 1991, NuPhA, 535, 331,doi: 10.1016/0375-9474(91)90452-CLentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, ApJL,807, L31, doi: 10.1088/2041-8205/807/2/L31Mabanta, Q. A., & Murphy, J. W. 2018, ApJ, 856, 22,doi: 10.3847/1538-4357/aaaec7Mabanta, Q. A., Murphy, J. W., & Dolence, J. C. 2019,ApJ, 887, 43, doi: 10.3847/1538-4357/ab4bccMarek, A., Dimmelmeier, H., Janka, H. T., M¨uller, E., &Buras, R. 2006, A&A, 445, 273,doi: 10.1051/0004-6361:20052840 Mathews, G. J., Boccioli, L., Hidaka, J., & Kajino, T. 2020,Modern Physics Letters A, 35, 2030011,doi: 10.1142/S0217732320300116Melson, T., Janka, H.-T., & Marek, A. 2015, ApJL, 801,L24, doi: 10.1088/2041-8205/801/2/L24Miller, D. S., Wilson, J. R., & Mayle, R. W. 1993, ApJ,415, 278, doi: 10.1086/173163M¨uller, B., Janka, H.-T., & Marek, A. 2013, TheAstrophysical Journal, 766, 43,doi: 10.1088/0004-637x/766/1/43M¨uller, B. 2016, PASA, 33, e048, doi: 10.1017/pasa.2016.40—. 2019, MNRAS, 487, 5304, doi: 10.1093/mnras/stz1594M¨uller, B., Janka, H.-T., & Marek, A. 2012, ApJ, 756, 84,doi: 10.1088/0004-637X/756/1/84M¨uller, B., Viallet, M., Heger, A., & Janka, H.-T. 2016,ApJ, 833, 124, doi: 10.3847/1538-4357/833/1/124M¨uller, B., Tauris, T. M., Heger, A., et al. 2019, MNRAS,484, 3307, doi: 10.1093/mnras/stz216Murphy, J. W., & Meakin, C. 2011, ApJ, 742, 74,doi: 10.1088/0004-637X/742/2/74O’Connor, E. 2015, ApJS, 219, 24,doi: 10.1088/0067-0049/219/2/24O’Connor, E., & Ott, C. D. 2010, Classical and QuantumGravity, 27, 114103,doi: 10.1088/0264-9381/27/11/114103—. 2011, ApJ, 730, 70, doi: 10.1088/0004-637X/730/2/70O’Connor, E., Bollig, R., Burrows, A., et al. 2018, Journalof Physics G: Nuclear and Particle Physics, 45, 104001,doi: 10.1088/1361-6471/aadeaeO’Connor, E. P., & Couch, S. M. 2018a, ApJ, 854, 63,doi: 10.3847/1538-4357/aaa893—. 2018b, ApJ, 865, 81, doi: 10.3847/1538-4357/aadcf7
D Turbulence in GR15