Including Neutrino-driven Convection in the Force Explosion Condition to Predict Explodability of Multidimensional Core-collapse Supernovae (FEC+)

Most massive stars end their lives with core collapse. However, it is not clear which explode as a core-collapse supernova (CCSN), leaving behind a neutron star, and which collapse to a black hole, aborting the explosion. One path to predict explodability without expensive multidimensional simulations is to develop analytic explosion conditions. These analytic conditions also provide a deeper understanding of the explosion mechanism and they provide some insight into why some simulations explode and some do not. The analytic force explosion condition (FEC) reproduces the explosion conditions of spherically symmetric CCSN simulations. In this follow-up manuscript, we include the dominant multidimensional effect that aids explosion—neutrino-driven convection—in the FEC. This generalized critical condition (FEC+) is suitable for multidimensional simulations and has potential to accurately predict explosion conditions of two- and three-dimensional CCSN simulations. We show that adding neutrino-driven convection reduces the critical condition by ∼30%, which is consistent with previous multidimensional simulations.


INTRODUCTION
At the end of their lives, massive stars undergo core collapse.Some explode as core-collapse supernovae (CCSNe), leaving behind a neutron star (NS) (Li et al. 2011;Horiuchi et al. 2011), some of them fail to explode and collapse to a black hole (BH) (Fischer et al. 2009;O'Connor & Ott 2011), and a small subset may explode and produce a BH through fallback (Ertl et al. 2016;Sukhbold et al. 2016;Ebinger et al. 2019).Understanding how and which massive stars explode has been a major challenge for decades.Predicting which massive stars explode is key to fully understanding nucleosynthesis (Woosley et al. 2002;Diehl et al. 2021), neutron star and black hole formation (Baade & Zwicky 1934;Burrows & Lattimer 1986;Fischer et al. 2009;O'Connor & Ott 2011), and neutron star and black hole distributions.
The CCSN mechanism is nonlinear and involves hydrodynamics, general relativity, radiation transport, and nuclear physics.Given this complexity in physics, one approach to understand the CCSN mechanism is to run expensive multi-dimensional simulations (Müller et al. 2017;Radice et al. 2017;Vartanyan et al. 2021).
Another approach is to develop analytic explosion conditions; these can provide a deeper understanding of the explosion mechanism and quantify how close simulations are to explosion.Over the years, many have tried to develop a theory for successful explosions (Thompson 2000;Pejcha & Thompson 2012;Raives et al. 2018Raives et al. , 2021;;Summa et al. 2016Summa et al. , 2018;;Yamasaki & Yamada 2005;Janka 2012;Janka et al. 2016;Keshet & Balberg 2012;O'Connor & Ott 2011;Wang et al. 2022;Boccioli et al. 2023).One of the first was the semi-analytic critical curve of Burrows & Goshy (1993).Burrows & Goshy (1993) assumed that the CCSN problem can be considered as a steady-state, boundary value problem and solved the ODEs with the inner boundary being at the NS and the outer boundary at the shock.They considered a two-dimensional parameter space of the neutrino luminosity L ν and the mass accretion rate, Ṁ and found that there is a critical curve in this parameter space below which there are stalled-shock solutions, but above this curve they found no stalled-shock solutions.They suggested that the region above this curve is the explosive region.Later, Murphy & Burrows (2008) showed that the semi-analytic critical curve of Burrows & Goshy (1993) is consistent with 1D simulations.To explain the physical origin of the critical curve of Burrows & Goshy (1993), Murphy & Dolence (2017) found that explosions occur when the integral of the momentum equation, Ψ is greater than 0. From this observation, Murphy & Dolence (2017) suggested a path toward deriving an analytic explosion condition.
Inspired by this suggestion, Gogilashvili & Murphy (2022) analytically derived an analytic critical condition for spherically symmetric explosions, called it the Force Explosion Condition (FEC) due to the fact that it is indeed the imbalance of forces that determine the explodability.Strikingly, they found that the explosion condition only depends upon two dimensionless parameters.Gogilashvili & Murphy (2022) first tested the consistency of the FEC with the semi-analytic solutions of Burrows & Goshy (1993).Then they checked the validity of the FEC using the light-bulb simulations (Gogilashvili & Murphy 2022) as well as using 1D radiation hydrodynamics simulations with the appropriate neutrino transport (Gogilashvili et al. 2023).In both cases, they find that the FEC accurately predicts the explodability of 1D simulations and is a good diagnostic of how far away a simulation is from explosion.While these successes of the FEC are encouraging, this simple model assumes spherical symmetry and is valid only for spherical explosions.Yet, a significant body of work shows that multi-dimensional effects such as neutrino-driven convection aids neutrinos in driving the explosion (Bethe 1990;Janka & Mueller 1996;Foglizzo et al. 2007;Mabanta & Murphy 2018).
In this followup manuscript, we generalize the FEC by including multi-dimensional effects in the model.In particular, we incorporate neutrino-driven convection and turbulent dissipation into the derivation of Gogilashvili & Murphy (2022).We show that adding neutrino-driven convection reduces the critical condition around ∼ 30% which is consistent with previous results from multi-dimensional simulations (Murphy & Burrows 2008;Murphy & Meakin 2011;Mabanta & Murphy 2018).
The structure of the manuscript is as follows: In section 2.1, we review the spherical model of the FEC of Gogilashvili & Murphy (2022).In section 2.2, we derive the generalized FEC (FEC+).In section 3, we review the applications of the FEC+ for different convection models.In particular, section 3.1 derives applications of the FEC+ for the integral convection model of Mabanta et al. (2019) and section 3.2-for the mixing-length-theory model of Supernova Turbulence In Reduced-dimensionality (STIR) (Couch et al. 2020).Lastly, in section 4, we discuss and summarize the main conclusion and provide next steps for using the FEC+.

A Brief Review of the Force Explosion Condition (FEC)
Gogilashvili & Murphy (2022) derived an analytic model for the explosion condition.The derivation is complete and self-consistent in that it begins with the fundamental equations of hydrodynamics and is explicit about all approximations.The force explosion condition (FEC) considers the balance of integrated forces and associated boundary conditions (Murphy & Dolence 2017).Using dimensional analysis, Gogilashvili & Murphy (2022) also demonstrate that the explosion condition only depends upon two dimensionless parameters...as opposed to four dimension-full parameters.
The force explosion condition is Lν τ g −aκ ≥ b.Lν τ g is the net neutrino power deposited in the gain region normalized by the accretion power: Lν τ g = L ν τ g R NS /(G Ṁ M NS ).κ is the dimensionless neutrino opacity that parameterizes the neutrino optical depth in the accreted matter near the neutron star surface: κ = κ Ṁ / √ GM NS R NS .The coefficient a is the difference of the post-shock pressure and the ram pressure of Ṁ , and b is related to gravity.Gogilashvili & Murphy (2022) analytically estimated these coefficients to be a ∼ 0.05 and b ∼ 0.52.Gogilashvili & Murphy (2022) tested the accuracy of the FEC in two ways.The first approach showed that the FEC is a more general expression of the semi-analytic condition of Burrows & Goshy (1993).Furthermore, this comparison led to a numerical fit for a = 0.06 and b = 0.38, which are remarkably close to the analytic estimates.The second test compares, the FEC with one-dimensional (1D) simulations that use local neutrino heating and cooling approximations (1D light-bulb simulations).The FEC accurately predicts the explodability of 1D Light-Bulb simulations (Figure 10 of Gogilashvili & Murphy (2022)).
In a followup manuscript, Gogilashvili et al. (2023) compared the FEC with 1D simulations that use actual neutrino transport.To make this comparison, they also proposed small modifications to the FEC.For example, replacing L ν τ g with the net neutrino heating rate in the gain region, Qν , is a more accurate measure of neutrino power for neutrino transport.They tested the FEC using GR1D, which simulates the spherically symmetric CCSN using GR and neutrino transport (O'Connor & Ott 2010;O'Connor 2015).The FEC accurately predicts the explosion condition for these spherically symmetric CCSN simulations that use neutrino transport (Figure 5 of Gogilashvili et al. (2023)).
To further generalize the FEC for two-and three-dimensional simulations, Gogilashvili & Murphy (2022) and Gogilashvili et al. (2023) proposed a way to add multi-dimensional effects such as neutrino-driven convection and turbulent dissipation to the analytic model of the FEC.Section 2.2 presents the derivation of FEC+, which is the FEC with neutrino-driven convection.

The FEC+
To develop a FEC for multi-dimensional simulations, one needs to include multi-dimensional effects in the simple spherical model.Since neutrino-driven convection is the dominant multi-dimensional instability that aids explosion (Mabanta & Murphy 2018), the following derivation incorporates neutrino-driven convection.Before adding a specific convection model, this section presents a more general derivation of FEC+ that accommodates any convective or turbulent, mean-field model.Much of the following derivation follows the procedure of Mabanta et al. (2019) (section 3.1).The derivation starts with the fundamental equations of hydrodynamics; using Einstein notation, they are: where ρ is the mass density, u is the fluid velocity, P is the pressure, h is the enthalpy, Φ is the gravitational potential, and q is net heating, which includes local specific neutrino heating, q ν = Lν κ 4πr 2 , and cooling, Reynolds decomposition separates the fluid variables into the background (subscript 0) and convective parts (superscript '), χ = χ 0 + χ ′ , where χ is a fluid variable.Following the same technique and assumptions, the Reynolds decomposed equations are (see details in appendix A and in Mabanta et al. ( 2019)): where R i j = u ′i u ′ j is the Reynolds stress, F i I = ρεu ′i is the perturbed internal energy flux, and W b = − ρ ′ u ′j Φ ,j is the work done by buoyant driving.To find an analytic explosion condition given these equations, we follow the prescriptions outlined in Murphy & Dolence (2017) and Gogilashvili & Murphy (2022).The first step is to assume steady-state; this transforms the evolution equations into a boundary value problem.The lower boundary is the NS "surface" or electron neutrinosphere, and the upper boundary is the shock.The second assumption is to assume spherical symmetry for both the background quantities and the convective flow.To facilitate analytic and algebraic solutions, we integrate the steady-state equations from the NS surface to the shock.
Boundary conditions and an equation of state (EoS) complete these equations.To facilitate analytic solutions, we assume that the pre-shock material is in free fall and pressureless (P 0+ = 0).We also assume that Ṁ is constant.Moreover, we assume that the Reynolds stress for pre-shock material is negligible (⟨R r r ⟩ ≈ 0).The resulting shock jump conditions are: A γ−law EoS, P = (γ − 1)ρε further facilitates analytic solutions.Eqs. ( 4)-( 6) along with the EoS and boundary conditions (eqs.( 7)-( 9)) fully describe the evolution of spherically symmetric background flow with neutrino-driven convection.
Gogilashvili & Murphy (2022) derived the FEC condition from the (im)balance of integral forces in the momentum equation.Here, we derive the integral forces including Reynolds-decomposed terms.The resulting balance of forces, Ψ, is The first three terms on the right side of eq. ( 10) are the same terms as in Gogilashvili & Murphy (2022).Integrating the last two terms, on the other hand, requires additional assumptions for the density profile and Reynolds stress.We scale the density to the density at the surface of the NS, ρ = ρ NS f (x), where x = r/R NS .Previous simulations show that the average radial component of the Reynolds stress is zero at the surface of the NS, ⟨R r r ⟩ NS ≈ 0 (see Murphy et al. (2013)) and for simplicity, it is reasonable to assume that the Reynolds stress is constant in the gain region, ⟨R r r ⟩ = ⟨R r r ⟩ − .With these assumptions, the last two terms are: and where x s and x g are dimensionless shock radius and gain radius respectively.Including the jump conditions (eqs.( 7)&(8)), denoting the shock compression ratio as β = ρ − /ρ + , and dropping subscripts "0" and "r", Ψ becomes: Murphy & Dolence (2017) and Gogilashvili & Murphy (2022) showed that Ψ ≥ 0 corresponds to the critical condition for explodability.Therefore taking Ψ = 0 in eq. ( 13) leads to the following critical condition in the dimensionless form (see details in Appendix A) where εNS is the dimensionless specific internal energy measured at R NS and Rr r = R r r * R NS /(GM NS ) is the dimensionless Reynolds stress.Coefficient a gives the relative scale of the difference in the pre-and post-shock pressures, coefficient b gives the relative scale of the Reynolds Stress term, and coefficient c gives the relative scale of gravity in this critical condition.Appendix A shows the analytic derivation for each coefficient; we estimate the scales to be a ≈ 0.04, b ≈ 0.38, and c ≈ 0.4.
To express the critical condition in terms of neutrino power, we next integrate the energy equation, eq. ( 6).There are two additional terms due to convection, one is related to the buoyant force, W b , and other to the perturbed internal energy flux, F I (Mabanta et al. 2019).Murphy & Meakin (2011) where E k =

Rs
Rg ρϵ k dV is the total power of turbulent dissipation.Integrating energy conservation, eq. ( 6), leads to an expression for the the specific internal energy in terms of the neutrino luminosity, L ν : This expressions uses several additional constraints: 1) it uses the analytic boundary condition, eq. ( 9), 2) B(x s ) is the dimensionless Bernoulli integral at the shock, 3) and in the cooling region, neutrino cooling is balanced by neutrino heating and gravitational heating (Gogilashvili & Murphy 2021).Ẽk = E k R NS /(G Ṁ M NS ) is the dimensionless turbulent energy dissipation rate.Finally, plugging eq. ( 16) into eq.( 14) gives the FEC in terms of neutrino luminosity: where a ′ = aγ, b ′ = bγ, and c ′ = cγ − B(x s ) − 1 xg .Eq. ( 17) represents an analytic force explosion condition (FEC+) for non-spherical explosions.Compared to the 1D FEC of Gogilashvili & Murphy (2022), FEC+ includes the effects of neutrino-driven convection through turbulent dissipation and turbulent ram pressure.Since these convective terms are positive definite, the convection model reduces the neutrino power required for explosion.
To illustrate the quantitative effect of the convection model, we estimate the scale of each term in eq. ( 17) as follows.To provide a basis for comparison, the neutrino power near the explosion condition is Lν τ g ∼ 0.4 (Gogilashvili & Murphy 2022), the difference-in-pressures term is a ′ κ ∼ 0.05.Multi-dimensional simulations indicate that the dimensionless turbulent dissipation is, Ẽk ∼ 0.12 and the dimensionless Reynolds stress is Rr r ∼ 0.02 (Murphy et al. 2013).Based on these estimates, turbulent dissipation is ∼ 30% of the net neutrino power while turbulent ram pressure term is only ∼ 5% (see appendix B).In summary, the convection model reduces the critical condition mostly through turbulent dissipation, and these analytic calculations quantitatively match the numerical results of (Mabanta & Murphy 2018).

The Integral Convection Model
Eq. ( 17) represents the FEC+ when considering convection models in general.This section considers the integral, mean field model of Murphy et al. (2013).While multi-dimensional effects such as neutrino-driven convection and turbulent dissipation are important in the shock revival, simulating two-and three-dimensional CCSN simulations are computationally very expensive, making it difficult to perform systematic studies.Therefore, as it is often the case, the multi-dimensional effects are in some way incorporated into 1D simulations (Murphy & Meakin 2011;Murphy et al. 2013;Perego et al. 2015;Mabanta et al. 2019;Couch et al. 2020).One can apply these specific convection models to further simplify the FEC+.
The integral convection model of Murphy et al. (2013) starts with fundamental equations of hydrodynamics and Reynolds decomposes physical variables into the background and turbulent parts.Based on multi-dimensional simulations, they show that many of the turbulent correlations are negligible and the ones that are important are the Reynolds stress, R, the work done by buoyant force, W b , and the turbulent energy flux, F I .To close the system of equations, one needs to know the relation between these quantities.There are several relationships that close the equations.For one, multi-dimensional simulations (Murphy et al. 2013) indicate that buoyant driving roughly balances turbulent dissipation, W b ≈ E k .Another relationship is that the turbulent energy dissipation and the turbulent flux is proportional to the neutrino power absorbed in the gain region: L max e ≈ αL ν τ and E k ≈ ζL ν τ (Murphy et al. 2013).Mabanta et al. (2019) incorporated this integral convection model into 1D simulations and showed that these 1D+ simulations accurately reproduce the explosion conditions of two-and three-dimensional simulations.To relate the Reynolds stress and turbulent dissipation, Mabanta et al. (2019) assume that where L is the largest turbulent eddy size and E k = rs rg ϵ k ρdV .On the other hand ϵ k can be approximated as where M g is the mass in the gain region.Eqs. ( 18) & ( 19) provide a relationship between the Reynolds stress and turbulent dissipation: (20) If one applies these assumptions, the FEC+ (eq.17) becomes: where the ratio of integrated turbulent dissipation to neutrino power is

Supernova Turbulence In Reduced-dimensionality (STIR)
Another convection model is the STIR model of Couch et al. (2020).The STIR method is similar in intent to the integral method in Mabanta et al. (2019), but there are some important differences.Both STIR and the integral method in Mabanta et al. (2019) relate the Reynolds stress and turbulent dissipation via eq.( 18).One major difference is that STIR uses a mixing length prescription to close the convection model.In particular, STIR assumes that the higher order Reynolds correlations are proportional to background gradients.For example, STIR assumes that the buoyant force is proportional to the Brunt-Väisä frequency: where ω BV is the Brunt-Väisälä frequency and Λ mix = α λ p/ρg is the mixing length with mixing length parameter α λ .
Assuming the Ledoux criterion, the Brunt-Väisälä frequency is where g eff = − ∂Φ ∂r +u r ∂ur ∂r is effective gravitational acceleration and c s is the adiabatic sound speed.Another difference is that STIR includes turbulent production from the divergence of the background flow: ∇ • (u 0 ⟨ρ(T rR)⟩) which is ∼ u ′2 .This term is ignored both in the integral convection model of Mabanta et al. (2019) as well as in the derivation of the FEC+ in this manuscript (2).Another major difference is that Mabanta et al. (2019) assumes steady state and presented an integrated convection model.In contrast, Couch et al. (2020) evolve the turbulent energy equation.
These differences lead to a different formulation for FEC+.In the following, we include ∇ • (u 0 ⟨ρ(T rR)⟩) term and solve for the FEC+.This additional term in the energy conservation, eq. ( 3) leads to an additional term in the FEC+ of ⟨R r r ⟩.Finally, since STIR solves for the Reynolds stress, a form that is more appropriate for STIR simulations is: 4. DISCUSSION AND CONCLUSIONS Understanding the explosion mechanism of CCSNe is important in understanding which massive stars collapse to black holes and which explode.The 1D FEC of Gogilashvili & Murphy (2022) accurately predicts the explodability of 1D radiation-hydrodynamic simulations (Gogilashvili et al. 2023).However, CCSN explosions are not spherical, and in fact, multi-dimensional effects such as neutrino-driven convection are important in aiding explosion.To be be able to predict the outcome of actual CCSN explosions in Nature, we generalize the FEC to include neutrino-driven convection.Eq. ( 17) presents this more general form, which we call FEC+.Including convection and turbulent dissipation generates two additional terms compared to the 1D FEC1 : the first is the dimensionless turbulent dissipation, Ẽk ; the second is due to turbulent ram pressure, b ′ Rr r .At first glance, there appear to be four dimensionless parameters in FEC+, two additional from the convection model.However, Kolmogorov's theory for turbulent dissipation relates the Reynolds stress to turbulent dissipation.This in turn reduces the number of dimensionless parameters.Therefore, the FEC+ effectively depends upon three dimensionless parameters: 1) the dimensionless net neutrino heating deposited in the gain region, 2) the dimensionless neutrino opacity, and 3) the dimensionless Reynolds stress (or turbulent velocity).
These additional terms due to convection and turbulent dissipation reduce the net neutrino heating required for explosion, and this reduction is consistent with published studies.The major contributor of this reduction is the turbulent dissipation rate which is ∼ 30% of net neutrino heating deposited in the gain region.The turbulent ram pressure has a smaller effect on the reduction as it is only ∼ 5% of total net neutrino heating.Taken together, we estimate that the convective model reduces the neutrino heating required for explosion by 26%.This analytic estimate is consistent with numerical results.For one, Mabanta & Murphy (2018) also found that turbulent dissipation dominates the reduction of the critical condition for explosion.Moreover, numerical simulations also show that convection reduces the critical condition by ∼30% (Murphy & Meakin 2011;Murphy et al. 2013;Mabanta & Murphy 2018).While these comparisons with previous studies are encouraging, it is important to check the predictions of FEC+ with multidimensional CCSN simulations.
If FEC+ proves to accurately predict the explosion conditions of multi-dimensional simulations, then this theory could significantly improve our understanding of how massive stars explode.For example, the FEC+ could quantify how close a multi-dimensional simulation is to explosion; it could also help to quantify the relative impact of the numerous physical processes involved.Additionally, the FEC+ suggests a clear path to predicting which stars explode from the progenitor structure alone, without performing expensive simulations.The dimensionless parameters depend upon dimension-full parameters such as L ν , M NS , R NS , and Ṁ .If one can predict the evolution of these quantities from the progenitor structure alone, then one could determine whether a progenitor structure will explode or not.Reducing the question of whether real stars explode in Nature to just three dimensionless numbers illustrates the elucidating potential of the FEC+.
Integrating to R s + ϵ includes the shock conditions in the equations.To eliminate the jump conditions, we then split these integrals into the integral up to shock and at the shock: This allow us to eliminate the boundary terms.
Therefore, the resulting Ψ parameter is Compared to Ψ for the spherically symmetric case (Gogilashvili & Murphy 2022), there are only two additional terms due to convection, the third and fourth terms in eq. ( A6).CCSN simulations show that in neutrino-driven convection, there is an equipartition between the radial direction and both tangential directions of Reynolds stress, R rr ∼ R θθ + R ϕϕ .In addition, the off-diagonal terms of the Reynolds stress are small (Murphy et al. 2013).This relation between the radial and tangential components of Reynolds stress allows us to further simplify the additional terms in the momentum equation, in particular: − 1 r ρ 0 R θ θ + ρ 0 R ϕ ϕ ≈ − 1 r ⟨ρ 0 R r r ⟩ (Murphy et al. 2013).These same simulations also show that ⟨R r r ⟩ is roughly constant in the gain region and quickly tapers to zero at the base of the gain region such that ⟨R r r ⟩ vanishes at the NS "surface".At the upper end of the convective region, ⟨R r r ⟩ remains nonzero all the way up to the shock.The density profile scale by the density at the NS surface is ρ = ρ NS f (x).These assumptions allow us to rewrite the integrals in terms of relevant scales such as ρ NS and an average Reynolds stress in the gain region.This in turn facilitates an algebraic and analytic condition.The nuances in the profile are absorbed into an integral of the dimensionless function f (x).Combing these assumptions, the third and fourth terms are: Dropping subscripts "0" and "r" and using boundary conditions, the third and fourth terms become: We normalize these terms in the same way as Gogilashvili & Murphy (2022) (eqs.( 7)-( 12) in Gogilashvili & Murphy (2022)).We also make the same assumptions as Gogilashvili & Murphy (2022); for example, we assume that Ṁ is For a typical mass of the NS, M NS ∼ 1.4M ⊙ , the radius of the NS, R NS ∼ 50km, and mass accretion rate, | Ṁ | ∼ 0.05M ⊙ , the dimensionless Reynolds stress and dimensionless turbulent dissipation are, Rr r ≲ 0.02 and Ẽk ∼ 0.12.Therefore, Ẽk / Lν τ g ∼ 30% while Rr r / Lν τ g ≲ 5%.Hence, the the reduction of the critical neutrino power is mostly due to turbulent dissipation.The contribution of the turbulent ram pressure is small.