How Turbulence Enables Core-collapse Supernova Explosions

Understanding how massive stars end their lives still remains an important astrophysical problem. Observations indicate that most massive stars likely explode as core-collapse supernovae (Horiuchi et al. 2011; Li et al. 2011), but the theoretical details of the explosion mechanism are still uncertain. The most certain part of the theory is that the Fe core collapses, bounces at nuclear densities, forming a proto-neutron star, and launches a shock wave (Janka et al. 2016). However, this shock wave quickly stalls (Hillebrandt & Mueller 1981; Mazurek 1982; Mazurek et al. 1982). For over three decades, the most prominent mechanism has been the delayed-neutrino mechanism in which neutrinos reheat the matter below the stalled shock and relaunch it into an explosive blast wave (Bethe & Wilson 1985). For all but the least massive stars, the neutrino mechanism fails in one-dimensional, spherically symmetric simulations (Liebendörfer et al. 2001a, 2001b, 2005; Rampp & Janka 2002; Buras et al. 2003, 2006; Thompson et al. 2003; Kitaura et al. 2006; Radice et al. 2017), but multidimensional simulations do explode (Herant et al. 1994; Burrows et al. 1995, 2007; Janka & Müller 1995, 1996; Dolence et al. 2015; Melson et al. 2015; Bruenn et al. 2016; Müller 2016; Roberts et al. 2016), and initial analyses suggest that multidimensional instabilities and turbulence aid the neutrino mechanism toward explosion (Murphy & Burrows 2008; Marek & Janka 2009; Murphy & Meakin 2011; Hanke et al. 2012; Murphy et al. 2013; Couch & Ott 2015; Melson et al. 2015; Radice et al. 2016). Therefore, understanding the explosion mechanism requires an understanding of the conditions between failed and successful explosions and how turbulence aids the explosion. There appears to be a critical condition for explosion (Burrows & Goshy 1993; Murphy & Burrows 2008; Murphy & Dolence 2017), and the critical condition for explosion is easier to obtain in multidimensional simulations (Herant et al. 1994; Burrows et al. 1995; Janka & Müller 1995, 1996; Murphy & Burrows 2008; Melson et al. 2015). In this paper, we propose an analytic model for turbulence and investigate how it reduces the condition for explosion.

For a thorough review on the status and problems of CCSN theory, see Müller (2017) and Janka et al. (2016). Here, we motivate our work with some of the most salient points. We know that the core collapses, but not yet how the collapse reverses into explosion. Whether or not a massive star explodes, the Fe core collapses at the end of the massive star's life. Prior to collapse, the overlying Si-burning layer adds Fe onto the iron core. As the iron core grows in size, the electrons which supply much of the electron degeneracy pressure become more and more relativistic. As the core nears the Chandrasekhar mass limit, the relativistic electron degeneracy pressure becomes less effective at supporting the core against gravitational collapse. Meanwhile, neutrino losses reduce the lepton number, decreasing the number of electrons available to supply pressure. The iron core becomes gravitationally unstable and contracts down to a sea of nucleons, forming a proto-neutron star, and is thus supported mostly by the strong force. At these nuclear densities, the equation of state for the core stiffens, and the core abruptly bounces, slamming into the rest of the star, which is collapsing supersonically onto the bouncing proto-neutron star. This creates an outward moving shock wave. As the shock wave propagates outward, it loses energy via photodissociation of Fe, electron capture, and neutrino losses. The shock stalls into an accretion shock, but if the star is to explode, this stalled accretion shock must relaunch into an explosive blast wave.

One proposed solution to the relaunch of the blast wave is the delayed-neutrino mechanism (Bethe & Wilson 1985). During the stalled-shock phase, the neutron star is cooling with a neutrino luminosity of a few ×10⁵² erg s⁻¹, but only about 10% of this luminosity is recaptured in a net heating region, the gain region. This volume is above the proto-neutron star but below the shock. If neutrino heating were the only effect, then this would be sufficient to relaunch the explosion. However, the region below the shock is in sonic contact, and so the structure satisfies a boundary-value problem. Although the neutrino heating adds heat that would drive explosion, the neutrino cooling at the base and the ram pressure of matter accreting through the shock keep the shock stalled. The hope has been that for high-enough neutrino luminosities or low accretion rates, the neutrinos might overwhelm the ram pressure and relaunch the explosion. Alas, spherically symmetric simulations show that this mechanism fails in in all but the least massive progenitors (Liebendörfer et al. 2001a, 2001b, 2005; Rampp & Janka 2002; Buras et al. 2003, 2006; Thompson et al. 2003; Kitaura et al. 2006; Müller et al. 2017; Radice et al. 2017).

While one-dimensional simulations fail, many multidimensional simulations seem to succeed (sometimes weakly; Herant et al. 1994; Burrows et al. 1995, 2007; Janka & Müller 1995, 1996; Dolence et al. 2015; Melson et al. 2015; Bruenn et al. 2016; Burrows et al. 2018; Müller 2016; Roberts et al. 2016), and the best indication is that turbulence plays an important role in aiding the delayed-neutrino mechanism toward explosion (Bethe & Wilson 1985; Janka & Müller 1996; Marek & Janka 2009; Murphy & Meakin 2011; Murphy et al. 2013). Therefore, to truly understand the explosion mechanism of massive stars, we need to identify the conditions for explosion and how turbulence affects these conditions.

There have been many attempts to characterize these conditions, some empirical (O'Connor & Ott 2011; Ott et al. 2013; Ertl et al. 2016), some heuristic (Janka & Keil 1998; Thompson 2000; Thompson et al. 2005; Buras et al. 2006; Müller et al. 2016), and others attempting to derive a condition from fundamentals (Burrows & Goshy 1993; Pejcha & Thompson 2012; Murphy & Dolence 2017). Of these, the most illuminating to date has been the critical neutrino-luminosity and accretion-rate curve (Burrows & Goshy 1993). During the stalled-shock phase, Burrows & Goshy (1993) noted that one may derive steady-state solutions for the stalled-shock structure. The governing equations describe a boundary-value problem in which the lower boundary is set by the neutron star surface (the neutrinosphere) and the outer boundary is the shock. They parameterized the problem in terms of the accretion rate, $\dot{{ \mathcal M }}$, onto the shock, and the neutrino luminosity, L_ν, emanating from the core. In this two-dimensional parameter space, they found steady-state solutions below a critical curve. Above this curve, they did not find steady-state stalled solutions, and suggested, but did not prove, that the solutions above the critical neutrino-luminosity and accretion-rate curve are dynamic and explosive.

Murphy & Dolence (2017) take a step closer to proving that the solutions above the curve are explosive by showing that the only steady solutions above this curve indeed have a positive shock velocity (${v}_{s}\gt 0$). Burrows & Goshy (1993) focused on just $\dot{{ \mathcal M }}$ and L_ν, but there are five parameters that define the steady solutions. These are the neutrino luminosity (L_ν), mass-accretion rate ($\dot{{ \mathcal M }}$), neutron star mass (M_NS), neutron star radius (R_NS), and the neutrino temperature (T_ν). Murphy & Dolence (2017) point out that the critical condition is not a critical curve but a critical hypersurface; most importantly, they find that this critical hypersurface is described by one dimensionless parameter, Ψ. In essence, Ψ is an integral condition related to the balance of pressure and gravity behind the shock. For a given set of the five parameters, Ψ may be negative, zero, or positive, which correspond to v_s < 0, v_s = 0, and v_s > 0, respectively. There is always a minimum Ψ, and if ${{\rm{\Psi }}}_{\min }\lt 0$, then there is always a stable steady-state stalled solution such that Ψ = 0. The critical condition is where Ψ_min = 0. Above this critical condition, Ψ_min > 0, and all steady solutions have v_s > 0. Assuming that these v_s > 0 steady solutions correspond to explosion, they use Ψ_min to define an explodability parameter.

Using one- and two-dimensional simulations, Murphy & Burrows (2008) empirically confirmed that criticality is a useful condition for explosion in core-collapse simulations. Furthermore, they found that the critical condition is ∼30% lower in two-dimensional simulations. Subsequently, others have confirmed these findings and that the reduction is similar in three-dimensional simulations (Hanke et al. 2012; Dolence et al. 2013; Handy et al. 2014; Fernández 2015).

Initial indications show that turbulence causes this reduction, but these investigations were mostly suggestive and not conclusive. Decades ago, Bethe (1990) recognized the potential importance of neutrino-driven convection aiding explosion. This initial investigation suggested that turbulent ram pressure behind the material would push against the shock. In the early 1990s, the first two-dimensional simulations with crude neutrino transport exploded while the one-dimensional simulations did not. These investigators speculated that neutrino-driven convection aided the explosion (Burrows et al. 1995; Janka & Müller 1996; Janka 2001; Colgate & Herant 2004).

Blondin et al. (2003) identified a new instability that can also drive turbulence: the standing accretion shock instability (SASI). Linear analyses suggest that this instability results from an advective–acoustic feedback cycle (Foglizzo & Tagger 2000; Foglizzo et al. 2006; Sato et al. 2009a, 2009b; Guilet et al. 2010), and subsequent investigations considered the possibility that the SASI aids the delayed-neutrino mechanism toward explosion (Marek & Janka 2009; Hanke et al. 2012, 2013; Fernández et al. 2014). Instead, Murphy et al. (2013) found that in simplified simulations, convection dominates just before and during neutrino-driven explosions. Further analyses with less simplistic neutrino approaches found that the SASI does dominate at times, but convection likely dominates for most but not all explosion conditions (Murphy & Meakin 2011; Burrows et al. 2012; Murphy et al. 2013; Müller 2016; Radice et al. 2016).

In one attempt to find the reason why turbulence aids explosion, Murphy & Burrows (2008) found that entropy is higher in the multidimensional case. At the time, they as well as others suggested that a longer dwell time in the gain region causes this higher entropy (Thompson et al. 2005; Buras et al. 2006; Marek & Janka 2009); however, these investigations did not show that these dwell-time distributions actually led to enhanced entropy in the gain region. Couch & Ott (2015) revisited the idea of turbulent ram pressure being the main multidimensional contribution, but speculated it as an effective ram pressure, not distinguishing between multidimensional effects. At this point, many of these suggestions seem plausible, but it is not clear which of these, if any, can explain why turbulence aids explosions. In fact, we will show that none of these explanations truly captured the role of turbulence. However, we do note that the higher entropy profiles of multidimensional turbulence should have hinted that turbulent dissipation (TD) plays an important role.

Part of the reason why it was difficult to assess how turbulence aids explosion is the complexity of multidimensional simulations. In these large, nonlinear simulations, it is difficult to isolate the causes and effects of turbulence. In this paper, we propose a different, more illuminating approach. We model turbulence in the context of the critical condition. Because we have direct control over how turbulence affects the equations, we can directly assess the causes and effects of our turbulence model in reducing the critical condition.

To do this, we extend the explodability parameter of Murphy & Dolence (2017) to multiple dimensions by including a turbulence model. Currently, a turbulence model exists for neutrino-driven convection (Murphy & Meakin 2011; Murphy et al. 2013), but one does not exist yet for the SASI. Thus, the only viable avenue for an analytic investigation of turbulence is through neutrino-driven convection. Investigations on how turbulence driven by the SASI affects the condition will have to wait until we have a valid turbulence model for the SASI. Therefore, we use the neutrino-driven convection turbulence model of Murphy et al. (2013) in the critical condition of Murphy & Dolence (2017).

In Section 2, we revisit criticality, present the neutrino-driven turbulence model, derive the equations, and describe the solution method. Furthermore, we derive an analytic upper limit on the Reynolds stress. In Section 3, we present how turbulence modifies the structure of the post-shock flow and how it affects the critical condition. Finally, in Section 4, we conclude and discuss implications for the core-collapse mechanism and future investigations. Our main conclusion are that turbulent ram pressure is not the primary turbulent term aiding explosions, but rather, it is one of a few, and that the dominant turbulent effect aiding explosions is TD.

In this section, we outline the method for deriving the critical condition including the effects of neutrino-driven convection. Our primary goal is to incorporate a turbulence model in the integral explosion condition of Murphy & Dolence (2017), which is a generalization of the foundational work of Burrows & Goshy (1993).

Our attempt is not the first to include turbulence in calculations of critical conditions. Yamasaki & Yamada (2005, 2006) considered several effects that might reduce the explosion condition, including rotation and convection. However, their attempt to model turbulence only included a term representing turbulent enthalpy flux. A more self-consistent derivation of a turbulence model should include three terms at a minimum: the turbulent enthalpy flux, the turbulent ram pressure, and TD. Without including all of these effects, it is unclear how turbulence would actually affect the critical condition for explosion. In this paper, we incorporate the turbulence model of Murphy et al. (2013), which includes these three terms and has been validated with core-collapse simulations.

First, in Section 2.1, we revisit criticality and present the relevant equations and assumptions. Next, in Section 2.2, we present our methods for incorporating turbulence into these equations. To find the solutions for both the average background flow and the turbulent flow, we decompose the equations into the average and turbulent quantities, commonly called Reynolds decomposition. By decomposing the variables into average and turbulent flows, we introduce extra variables. To close the system of equations, we motivate a turbulence model in Section 2.3, and we specify our assumptions. In Section 2.4, we describe our method for solving this system of equations. As it turns out, there is a maximum allowable Reynolds stress in the stalled-shock solutions. Even though this upper limit does not affect the critical condition, we derive an analytic expression for it in Section 2.5 and include it in our critical condition calculations. Finally, in Section 2.6, we discuss some of our assumptions and parameters that need to be verified by multidimensional simulations

2.1. The Critical Curve

Two explosion conditions which start from first principles are the critical curve of Burrows & Goshy (1993) and a generalization of this condition, the Ψ condition of Murphy & Dolence (2017). We will modify these critical conditions to include turbulence, so first we revisit what constitutes a critical condition.

Burrows & Goshy (1993) introduced a critical curve that divides steady-state solutions from no solutions in the ${L}_{\nu }-\dot{{ \mathcal M }}$ plane. Below this curve, one can find steady-state stalled-shock solutions that satisfy all boundary conditions. Above this curve, there are no stalled-shock solutions that satisfy all of the boundary conditions. In detail, the main discriminant for finding steady-state solutions is whether or not they could match the density at the neutrinosphere with the stalled-shock solution. They suggested, but did not prove, that the solutions above the ${L}_{\nu }-\dot{{ \mathcal M }}$ curve are explosive.

Murphy & Dolence (2017) took one step closer in proving that the solutions above the critical curve are explosive. They showed that the only steady-state solutions above the curve have v_s > 0, lending support to the supposition of Burrows & Goshy (1993). They did this by connecting the discriminant of Burrows & Goshy (1993; the τ = 2/3 condition) to a new dimensionless parameter Ψ. This Ψ parameter is a measure of overpressure compared to hydrostatic equilibrium and directly indicates whether the shock would move out or in, or stay stationary. Because it is more straightforward, we do our calculations in the Burrows & Goshy (1993) formalism, using the neutrino density discriminant, but we will also present our results in the context of Ψ, the overpressure. Although the method of Burrows & Goshy (1993) is simpler, Murphy & Dolence (2017) provides a more direct connection to v_s, and thus to explosion.

The first step in finding the critical curve is finding steady-state solutions to the Euler equations. The steady-state equations are

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\rho u)=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\rho u\otimes u)=-{\rm{\nabla }}P-\rho {\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot \left[\rho u\left(h+\displaystyle \frac{{u}^{2}}{2}\right)\right]=-\rho u\cdot {\rm{\nabla }}{\rm{\Phi }}+\rho q.\end{eqnarray} \tag{ 3 }$$

where ρ is the mass density, u is the velocity, P is the pressure, Φ is the gravitational potential, h is the enthalpy, and q is the total heating. In general, heating and cooling by neutrinos is best described by neutrino transport (Janka 2017; Tamborra et al. 2017); we simplify neutrino transport by invoking a simple light-bulb prescription for neutrino heating and a local cooling (Janka 2001),

$$\begin{eqnarray}&&q=\displaystyle \frac{{L}_{\nu }\kappa }{4\pi {r}^{2}}-{C}_{0}\,{\left(\displaystyle \frac{T}{2\mathrm{MeV}}\right)}^{6}.\end{eqnarray} \tag{ 4 }$$

That is, we have adopted a spherically symmetric neutrino source. r is the radius from this source, L_ν is the neutrino luminosity emitted from the core of the star, and

$$\begin{eqnarray}&&\kappa =7.5\times {10}^{-28}\,{\left(\displaystyle \frac{{T}_{\nu }}{4\mathrm{MeV}}\right)}^{2}({Y}_{n}+{Y}_{p})[{\rm{g}}\cdot {\mathrm{cm}}^{2}]\end{eqnarray} \tag{ 5 }$$

is the neutrino opacity (Murphy et al. 2013), where Y_n and Y_p are the neutrino and proton fractions per baryon. κ is related to the optical depth by

$$\begin{eqnarray}&&\tau ={\int }_{{r}_{g}}^{{r}_{s}}\rho \kappa {dr},\end{eqnarray} \tag{ 6 }$$

which dictates the absorptivity of the material in the gain region; r_s and r_g are the positions of the shock and gain radius, respectively. T is the matter temperature, and C₀ is the cooling factor (1.399 × 10²⁰ erg g⁻¹ s⁻¹).

Following in the steps of Burrows & Goshy (1993), Equations (1)–(3) represent a boundary-value problem with the boundaries being the neutron star surface and the shock. At the shock, we want the pre-shock, inflowing material to match the post-shock material through the Rankine–Hugoniot jump conditions. For the moment, let us assume that v_s is zero, and in this case, the jump conditions become

$$\begin{eqnarray}&&{\rho }_{1}{u}_{1}={\rho }_{2}{u}_{2},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rho }_{1}{u}_{1}^{2}+{P}_{1}={\rho }_{2}{u}_{2}^{2}+{P}_{2},\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&{\rho }_{1}{u}_{1}\left({e}_{1}+\displaystyle \frac{1}{2}{u}_{1}^{2}+\displaystyle \frac{{P}_{1}}{{\rho }_{1}}\right)={\rho }_{2}{u}_{2}\left({e}_{2}+\displaystyle \frac{1}{2}{u}_{2}^{2}+\displaystyle \frac{{P}_{2}}{{\rho }_{2}}\right),\end{eqnarray} \tag{ 9 }$$

where e is the internal energy and the subscripts 1 and 2 indicate the downstream and upstream flows at the jump, respectively. Normally, one uses Equation (7) in Equation (9) to eliminate the mass flux in the Hugoniot–Rankine jump condition. Here, we explicitly include it because this term cannot be ignored when we Reynolds-decompose and derive the jump conditions including the turbulent terms. One then integrates inward to the neutron star surface where the density profile must match the neutron star surface such that the neutrino optical depth is 2/3. If the neutrino optical depth is not 2/3, then one searches for a new r_s so that the shock and neutron star boundary conditions are met. In practice, Yamasaki & Yamada (2005) noticed that the density at the neutrinosphere has about the same value in most situations. For the opacities used in this paper, Murphy & Burrows (2008) calculated the density in one-dimensional and two-dimensional simulations to be about 7 × 10¹⁰ g cm⁻³ at the neutron star "surface." Since this density condition is faster to integrate, we use it instead of the neutrino optical depth condition.

Although Burrows & Goshy (1993) provide an elegant one-dimensional explosion condition, it does not accurately diagnose realistic supernova explosions in multiple dimensions (Murphy & Burrows 2008). However, subsequent multidimensional work suggests that one might be able to augment the technique for finding the critical curve to include turbulence (Murphy & Meakin 2011; Hanke et al. 2012; Couch 2013; Murphy et al. 2013; Radice et al. 2016). To do this, we build on the work of Murphy & Meakin (2011) and Murphy et al. (2013), and use the Reynolds-decomposed continuity equations to find a multidimensional critical curve.

2.2. Reynolds-decomposed Equations

A standard method to incorporate turbulence is through Reynolds decomposition. The primary goal is to derive mean-field, steady-state equations for turbulence. The first step is to Reynolds-decompose the variables into background and perturbed components, i.e., $u={u}_{0}+{u}^{{\prime} }$, where 0 denotes the background component and the prime indicates the turbulent term. Next, one substitutes these terms into the time-dependent Navier–Stokes equations. To obtain the mean-field correlations for turbulence, one averages the equations both in time and solid angle. For simplicity, we denote both of these averages by the operator $\langle \cdot \rangle$. Since $\langle u\rangle ={u}_{0}$ and $\langle {u}^{{\prime} }\rangle =0$, terms that involve one component of a turbulent variable are zero. All non-zero turbulent terms are higher-order correlations of turbulent variables.

Technically, the turbulence represents a time-dependent fluctuation. However, the mean-field variables or turbulent correlations are time-averaged correlations and can be in steady state. For core-collapse simulations, the turbulent correlations are effectively in steady state, so one may drop the time derivatives in the Reynolds-averaged equations. The resulting equations represent steady-state equations for the background flow and the mean-field turbulent correlations. In this manuscript, we highlight the important correlations and steady-state equations. For a more thorough derivation of the equations from the Navier–Stokes equations, see Meakin & Arnett (2007) or Murphy & Meakin (2011).

The three dominant Reynolds turbulent correlations are the Reynolds Stress (R), TD (${\epsilon }_{k}$), and turbulent luminosity (L_e; Murphy & Meakin 2011; Murphy et al. 2013). The Reynolds stress is the turbulent fluctuation in momentum stress, the turbulent luminosity is the transport of turbulent internal energy, and TD is the viscous conversion of mechanical energy to heat. These terms are

$$\begin{eqnarray}&&{R}_{{ij}}={u}_{i}^{{\prime} }{u}_{j}^{{\prime} },\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{L}_{e,i}=4\pi {r}^{2}{\rho }_{0}\langle {u}_{i}^{{\prime} }e^{\prime} \rangle =4\pi {r}^{2}{F}_{e,i},\end{eqnarray} \tag{ 11 }$$

and in the limit of small viscosity, TD is

$$\begin{eqnarray}&&{\boldsymbol{\epsilon }}=2\nu ({\rm{\nabla }}u^{\prime} )\cdot ({\rm{\nabla }}u^{\prime} ).\end{eqnarray} \tag{ 12 }$$

The turbulent kinetic energy dissipation is ${\epsilon }_{k}={tr}({\boldsymbol{\epsilon }})/2$. Note that this definition is slightly different from the Murphy & Meakin (2011) definition, which presented a confusing sign and needlessly included turbulent diffusion in the TD term. Here, we take the limit of small viscosity so that the turbulent diffusion term goes away, but TD remains. Furthermore, we corrected the sign so that a positive corresponds to taking energy from the kinetic energy equation and putting it in the internal energy equation. Since this term requires higher-order correlations, we model it using Kolmogorov's assumptions (see Section 2.3 or refer to the result of Canuto 1993 for a more robust description).

The resulting steady-state, Reynolds-decomposed equations are

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot ({\rho }_{0}{{\boldsymbol{u}}}_{0}+\langle \rho ^{\prime} {\boldsymbol{u}}^{\prime} \rangle )=0,\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\langle \rho {\boldsymbol{u}}\rangle \cdot {\rm{\nabla }}{{\boldsymbol{u}}}_{0}=-{\rm{\nabla }}{P}_{0}+{\rho }_{0}{\boldsymbol{g}}-{\rm{\nabla }}\cdot \langle \rho {\boldsymbol{R}}\rangle ,\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&\langle \rho u\rangle \cdot {\rm{\nabla }}{e}_{0}+\langle {P}_{0}{\rm{\nabla }}\cdot {u}_{0}\rangle +\langle P^{\prime} {\rm{\nabla }}\cdot u^{\prime} \rangle \\ &&\quad =\,-{\rm{\nabla }}\cdot {F}_{e}+{\rho }_{0}q+{\rho }_{0}{\epsilon }_{k}.\end{eqnarray} \tag{ 15 }$$

To see the exact equation, please refer to Meakin & Arnett (2007). There, they have fully expanded the above equation into their background and perturbed components. The internal energy flux and Reynolds stress are

$$\begin{eqnarray}&&{F}_{e}=\langle \rho {ue}^{\prime} \rangle \end{eqnarray} \tag{ 16a }$$

and

$$\begin{eqnarray}&&\langle \rho R\rangle =\langle \rho {uu}^{\prime} \rangle .\end{eqnarray} \tag{ 16b }$$

Alternative and common definitions are ${F}_{e}={\rho }_{0}\langle u^{\prime} e^{\prime} \rangle$ and $\langle \rho R\rangle ={\rho }_{0}\langle u^{\prime} u^{\prime} \rangle$. We use Equation 16(a) because it gives the same result but is much simpler and cleaner to calculate in numerical simulations. Expanding Equations 16(a) and (b) gives

$$\begin{eqnarray}&&{F}_{e}={\rho }_{0}\langle u^{\prime} e^{\prime} \rangle +\langle \rho ^{\prime} u^{\prime} e^{\prime} \rangle +{u}_{0}\langle \rho ^{\prime} e^{\prime} \rangle \end{eqnarray} \tag{ 17a }$$

and

$$\begin{eqnarray}&&\langle \rho R\rangle ={\rho }_{0}\langle u^{\prime} u^{\prime} \rangle +\langle \rho ^{\prime} u^{\prime} u^{\prime} \rangle +{u}_{0}\langle \rho ^{\prime} u^{\prime} \rangle .\end{eqnarray} \tag{ 17b }$$

Within the convective region, Murphy & Meakin (2011) and Murphy et al. (2013) found that the first term is the dominant term. However, just using the first term creates a large spike at the aspherical shock, which has nothing to do with convection and everything to do with the jump conditions across the aspherical shock. However, using Equation 16(a) mitigates this problem and gives the correct turbulent energy flux within the convective region.

The decomposed boundary conditions are

$$\begin{eqnarray}&&{\rho }_{1}{u}_{1}^{2}+{P}_{1}+{\rho }_{1}{R}_{{rr}}={\rho }_{2}{u}_{2}^{2}+{P}_{2}\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{1}}{{\rho }_{1}}+{e}_{1}+\displaystyle \frac{{L}_{e}}{\dot{{ \mathcal M }}}+\displaystyle \frac{{L}_{k}}{\dot{{ \mathcal M }}}+\displaystyle \frac{1}{2}{u}_{1}^{2}+{R}_{{rr}}=\displaystyle \frac{{P}_{2}}{{\rho }_{2}}+{e}_{2}+\displaystyle \frac{1}{2}{u}_{2}^{2},\end{eqnarray} \tag{ 19 }$$

where $\dot{{ \mathcal M }}=4\pi {r}^{2}\rho u$ is the mass-accretion rate, ${L}_{e}\,=4\pi {r}^{2}{\rho }_{0}\langle u^{\prime} e^{\prime} \rangle$ is the internal energy luminosity, and ${L}_{k}(={\rho }_{0}{u}_{0}\langle u{{\prime} }^{2}/2\rangle )$ is the kinetic energy luminosity. Now that we have introduced three new turbulent variables, we have a total of six unknown variables and only three equations, necessitating more equations.

2.3. Turbulence Models

Including the turbulent components, the Reynolds-decomposed conservation equations, Equations (13)–(15), now have more variables than equations. These extra variables are the Reynolds stress, turbulent luminosity, and TD. Therefore, to find a solution to these equations, we need a turbulence closure model. Turbulence depends upon the bulk macroscopic flow, so the equations for turbulence represent a boundary-value problem that depends upon the specifics of the background flow. For this reason, Murphy et al. (2013) developed a global turbulence model for neutrino-driven convection. The steady-state equations require local turbulent expressions and derivatives. Therefore, to use the global model, we must make some assumptions and translate the global model into a local model.

In the core-collapse problem, there may be two sources of turbulence: convection and SASI (Bethe 1990; Blondin et al. 2003). In principle, to correctly model turbulence in the core-collapse problem, we need a turbulence model that addresses both driving mechanisms: convection and SASI. Although there are nonlinear models to describe turbulent convection, there are no nonlinear models yet to describe SASI turbulence. Thus, we proceed with a convection-based analysis of turbulence.

There are five turbulent variables (${\boldsymbol{R}},{L}_{e},{\epsilon }_{k}$), and three of them are Reynolds stress terms (R_rr, ${R}_{\phi \phi }$, and ${R}_{\theta \theta }$); therefore, we invoke five constraints. Our five global constraints are as follows. First, we eliminate the tangential components of the Reynolds stress. In neutrino-driven convection, the radial Reynolds stress is in rough equipartition with both of the tangential components (Murphy et al. 2013),

$$\begin{eqnarray}&&{R}_{{rr}}\sim {R}_{\phi \phi }+{R}_{\theta \theta }.\end{eqnarray} \tag{ 20 }$$

Similar simulations showed that the transverse components are roughly the same scale,

$$\begin{eqnarray}&&{R}_{\phi \phi }\sim {R}_{\theta \theta }.\end{eqnarray} \tag{ 21 }$$

Using Kolmogorov's (1941) hypothesis, we relate the Reynolds stress to the TD,

$$\begin{eqnarray}&&{\epsilon }_{k}\approx \displaystyle \frac{u{{\prime} }^{3}}{{ \mathcal L }}=\displaystyle \frac{{R}_{{rr}}^{3/2}}{{ \mathcal L }},\end{eqnarray} \tag{ 22 }$$

where ${ \mathcal L }$ is the largest TD scale. From Murphy & Meakin (2011), we note that buoyant driving roughly balances TD,

$$\begin{eqnarray}&&{W}_{b}\approx {E}_{k},\end{eqnarray} \tag{ 23 }$$

where the buoyant driving is the total work done by buoyant forces in the convective region,

$$\begin{eqnarray}&&{W}_{b}={\int }_{{r}_{g}}^{{r}_{s}}\langle \rho ^{\prime} {u}_{i}^{{\prime} }\rangle {g}^{i}{dV},\end{eqnarray} \tag{ 24 }$$

and the total power of the dissipated turbulent energy is

$$\begin{eqnarray}&&{E}_{k}={\int }_{{r}_{g}}^{{r}_{s}}\rho {\epsilon }_{k}{dV}.\end{eqnarray} \tag{ 25 }$$

Lastly, three-dimensional simulations from Murphy et al. (2013) show that the source of neutrino-driven convection, the neutrino power, is related to the TD and the turbulent luminosity via

$$\begin{eqnarray}&&{L}_{\nu }\tau \approx {E}_{k}+{L}_{e}^{\max }.\end{eqnarray} \tag{ 26 }$$

Together, Equations (20)–(23) and (26) represent our turbulence closure model.

Now that we have combined a series of global conditions to close our global model, we must relate these back to local functions in order to incorporate them into our conservation equations. We do this by making assumptions about the local profile for each term and scaling them by one parameter for each turbulent term. In translating from global to local, we introduce three parameters of the turbulent region: a constant Reynolds stress (R), a constant dissipation rate (${\epsilon }_{k}$), and a maximum for the turbulent luminosity (L_e^max); the corresponding local terms are ${\rm{\nabla }}\cdot \langle \rho {\boldsymbol{R}}\rangle$, ${\rho }_{0}{\epsilon }_{k}$, and ${\rm{\nabla }}\cdot \langle {{\boldsymbol{L}}}_{e}\rangle ,$ respectively (see Equations (14)–(15)). Thus, the final solution for turbulence boils down to finding these three parameters.

To find the three parameters, we insert the profiles for the turbulent terms and their parameters into the global conditions, Equations (20)–(23) and (26). This then leads to a set of equations for the parameters, and we use simple algebra to solve for the scale of the parameters that satisfies those global conditions. We solve for these in the order they are presented: the Reynolds stress, TD, and finally the turbulent luminosity.

Our first task is to reduce the three Reynolds stress terms down to one. In neutrino-driven convection, there is a preferred direction (i.e., in the direction of gravity) and simulations show that there is an equipartition between the radial direction and both of the tangential directions (Murphy et al. 2013). Simulations also show that the two tangential directions have the same scale (Murphy et al. 2013). Evaluating these assumptions in Equations (20)–(21) reduces the representation of the Reynolds stress from three variables down to one,

$$\begin{eqnarray}&&R=2{R}_{{rr}}.\end{eqnarray} \tag{ 27 }$$

Kolmogorov's theory of turbulence predicts that the TD rate scales as the perturbed velocity cubed over the characteristic length of the instabilities, Equation (22). Numerical simulations suggest that the largest dissipation scale in convection is the size of the convective zone (Couch & O'Connor 2014; Fernández 2015; Foglizzo et al. 2015) or the gain region in the core-collapse case. Moreover, Equations (23)–(26) have a global definition of ${\epsilon }_{k}$, and so we assume that ${\epsilon }_{k}$ is roughly constant over the gain region. Therefore, from Equation (25),

$$\begin{eqnarray}&&{E}_{k}={\int }_{{r}_{g}}^{{r}_{s}}\rho {\epsilon }_{k}{dV},\end{eqnarray} \tag{ 28 }$$

and we have

$$\begin{eqnarray}&&{\epsilon }_{k}\approx \displaystyle \frac{{E}_{k}}{{M}_{\mathrm{gain}}}=\displaystyle \frac{{W}_{b}}{{M}_{\mathrm{gain}}}.\end{eqnarray} \tag{ 29 }$$

The final missing piece is the turbulent luminosity, L_e, and we connect this to the TD by rewriting the buoyant work in terms of the turbulent luminosity. We combine the observation that buoyant work is approximately equal to the TD power (Murphy et al. 2013) and that this energy can be converted to heat (Murphy & Meakin 2011). Ignoring compositional perturbations, the density perturbation in terms of the perturbed energy and pressure is

$$\begin{eqnarray}&&\rho ^{\prime} =e^{\prime} {\left(\displaystyle \frac{\partial \rho }{\partial e}\right)}_{P}+P^{\prime} {\left(\displaystyle \frac{\partial \rho }{\partial P}\right)}_{e}.\end{eqnarray} \tag{ 30 }$$

Convective flows are generally dominated by buoyant perturbations and not pressure perturbations. Even for high Mach number convection, buoyancy tends to dominate; see Murphy et al. (2013). Therefore, one may express the density perturbation in terms of the energy perturbation alone. Applying the above step and some algebra to Equation (24), we obtain

$$\begin{eqnarray}&&{W}_{b}=\int \displaystyle \frac{(\gamma -1){L}_{e}{\rm{\Phi }}\rho }{P}\displaystyle \frac{{dr}}{r}.\end{eqnarray} \tag{ 31 }$$

Now that we have an equation for the buoyant driving as a function of the turbulent luminosity (L_e), we now need the radial profile for L_e to complete the connection between turbulent luminosity and buoyant driving. Initial numerical investigations of L_e suggest that it rises from zero at the gain radius to a nearly constant value above the gain radius until the shock (Murphy et al. 2013). Therefore, we develop an ansatz for L_e which roughly satisfies the shape seen in simulations,

$$\begin{eqnarray}&&{L}_{e}={L}_{e}^{\max }\tanh \left(\displaystyle \frac{r-{r}_{g}}{{\boldsymbol{h}}}\right).\end{eqnarray} \tag{ 32 }$$

where h is the distance it takes for the turbulent luminosity to increase from zero to roughly L_e^max. Thus, our buoyant driving power becomes

$$\begin{eqnarray}&&{W}_{b}=(\gamma -1){L}_{e}^{\max }\int \tanh \left(\displaystyle \frac{r-{r}_{g}}{h}\right)\displaystyle \frac{{\rm{\Phi }}\rho }{P}\displaystyle \frac{{dr}}{r}.\end{eqnarray} \tag{ 33 }$$

Finally, relating this back to our global conditions (26) and (29), we can find the maximum turbulent luminosity,

$$\begin{eqnarray}&&{L}_{e}^{\max }\approx \displaystyle \frac{{L}_{\nu }\tau }{1+(\gamma -1)\int \tanh \left(\tfrac{r-{r}_{g}}{h}\right)\tfrac{{\rm{\Phi }}\rho }{P}\tfrac{{dr}}{r}}.\end{eqnarray} \tag{ 34 }$$

Within the context of our assumptions, Equation (34) is our final algebraic expression which gives the scale of the turbulence in terms of the background structure and the driving neutrino power. The second term in the denominator is a weak function of the shock radius and is of order unity. Therefore, the turbulent luminosity is a fraction of the neutrino-driving power ${L}_{e}^{\max }\approx {L}_{\nu }\tau /(1+f({R}_{s}))$, and the TD is also a fraction of the driving neutrino power ${E}_{k}\approx {L}_{\nu }\tau f({R}_{s})/(1+f({R}_{s}))$.

Though we have successfully closed our set of equations, we did employ several assumptions about the local structure of turbulence. Many of the assumptions made in this section have been verified by simulations individually (Murphy et al. 2013), and since we have been careful in being self-consistent, should have equal validity when combined. However, there are some assumptions that require further verification, and in Section 2.4, we discuss these details.

2.4. Solution Method Including Turbulence

2.5. A New Upper Bound on the Reynolds Stress

We have discovered that there is an upper bound on the Reynolds stress in the core-collapse problem. This comes from the fact that the Reynolds stress applies a turbulent ram pressure at the shock, and for large enough shock radii, there is an upper limit to the Reynolds stress that allows solutions to the stalled-shock jump conditions. Of course, the Reynolds stress depends upon the driving neutrino power, but otherwise we find that as one increases the shock radius, the Reynolds stress only goes up slightly; however, the scale of the ram pressure at the shock is set by the gravitational potential energy, which decreases as 1/r. Eventually, the ratio of the Reynolds stress to the gravitational potential becomes so large that there are no longer solutions to the steady-state jump conditions. We now derive the analytic upper bound to the Reynolds stress and describe our implementation of this limit into our solution method.

To quantify this upper bound, we start with our three jump conditions, Equations (7)–(9), that include turbulence, and make some assumptions so that we can derive an analytic expression. Our first assumption is that the fluxes of the internal energy and kinetic energy at the shock are small relative to the other terms (Meakin & Arnett 2007, 2010; Murphy & Meakin 2011). Our new decomposed boundary conditions are thus

$$\begin{eqnarray}&&{P}_{1}+{\rho }_{1}{u}_{1}^{2}+{\rho }_{1}{R}_{{rr}}={P}_{2}+{\rho }_{2}{u}_{2}^{2}\end{eqnarray} \tag{ 35 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{1}}{{\rho }_{1}}+{e}_{1}+\displaystyle \frac{1}{2}{u}_{1}^{2}+{R}_{{rr}}\approx \displaystyle \frac{{P}_{2}}{{\rho }_{2}}+{e}_{2}+\displaystyle \frac{1}{2}{u}_{2}^{2}.\end{eqnarray} \tag{ 36 }$$

Here we have omitted the 0 subscript, where all non-perturbed terms are implied to be the background. Since all of the perturbed components have either canceled or been defined, there is no need to differentiate with a "₀" or a "'".

Using a γ-law equation of state, we can solve for a solution of the ratio of densities,

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\gamma +1}{\gamma +{{ \mathcal M }}_{2}^{-2}-\sqrt{{(1-{{ \mathcal M }}_{2}^{-2})}^{2}-(\gamma +1)\tfrac{2{R}_{{rr}}}{{u}_{2}^{2}}}},\end{eqnarray} \tag{ 37 }$$

where β is the compression factor, $\beta ={\rho }_{1}/{\rho }_{2}$. For physical solutions of β, the term under the radical needs to be positive. This sets an upper limit on the second term, which is an upper limit on the Reynolds stress.

We now make some approximations to derive a simple analytic limit for the Reynolds stress. If we assume that the velocity of infalling matter onto the shock is roughly in free fall,

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{v}_{2}^{2}\approx {\rm{\Phi }}({r}_{s}),\end{eqnarray} \tag{ 38 }$$

an upper limit on our Reynolds stress becomes

$$\begin{eqnarray}&&{R}_{{rr}}\leqslant \displaystyle \frac{{(1-{{ \mathcal M }}_{2}^{-2})}^{2}}{\gamma +1}{\rm{\Phi }}({r}_{s}),\end{eqnarray} \tag{ 39 }$$

or in dimensionless form,

$$\begin{eqnarray}&&{ \mathcal R }=\displaystyle \frac{(\gamma +1){R}_{{rr}}}{{\rm{\Phi }}({r}_{s}){(1-{{ \mathcal M }}^{-2})}^{2}}\leqslant 1.\end{eqnarray} \tag{ 40 }$$

Above this limit, there can be no stalled-shock solutions, and since our method takes this assumption as an intermediate step, above this limit, we cannot find these quasi-steady solutions. Thus, we terminate our Ψ curve at the radius for which our ${ \mathcal R }$ parameter crosses this threshold. In practice, we had numerical difficulty when ${ \mathcal R }$ got close to one. So to avoid that numerical difficulty, we set a cap of 0.6 on this value (see Figure 4).

This upper limit on the Reynolds stress could have affected the critical curve; however, it does not. To reiterate: the critical neutrino luminosity curve is determined by the point of the minimum of the Ψ curve. Theoretically, if the imposed upper limit on ${ \mathcal R }$ were to be at a shock radius smaller than our Ψ_min, an entirely new critical curve might have to be defined in order to ensure that the above constraint was not being violated. Luckily, in all cases where an actual cap is necessary, this happens to be at a shock radius greater than Ψ_min. Thus, the cutoff point that we use, ∼60%, is synonymous with the condition of Equation (39). Though imposing an upper limit on R_rr would ostensibly mitigate its effect on the critical curve, we have just shown that all pseudo-solutions after Ψ_min are irrelevant. Hence, we are only constraining the Reynolds stress for pseudo-solutions, which are already non-steady state. That said, it is possible that once we start looking at the explodability of a set of initial conditions, our upper limit may start to interfere with the predictions. This upper limit will be treated accordingly with more rigor in future publications.

2.6. Discussion, Parameters, and Limitations of the Method

Our primary objective is to understand how turbulence affects the conditions for explosion. To fully understand this influence, we also need to understand how turbulence affects the background structure, so we first show how the turbulent terms affects the density and temperature profiles. We then show that turbulence raises the Ψ parameter, implying an easing of the explosion condition. We then consider how this affects the critical hypersurface. To connect with previous investigations, we focus on the neutrino-luminosity and accretion-rate slice of this critical hypersurface. We find that this reduction in the critical curve is ∼30%, in concordance with multidimensional simulations. To investigate how turbulence reduces the condition for explosion, we calculate the critical condition with each term included and omitted. Lastly, we provide evidence that our upper limit on the Reynolds stress does not affect the actual reduction of the critical curve.

In Figure 1, we illustrate how turbulence affects the density profile; in particular, we show the neutrino and convective terms in the derivative for the density. Since the density profile most closely matches a power law, we plot terms of d ln ρ/d ln r to compare how neutrinos and turbulence affect the slope in the log. To reduce the clutter, we do not plot all of the terms; for the full ODE, see Equation (48). In general, the missing terms give a slope of about −3. In the gain region, both neutrinos and convective terms make the slope shallower. In the cooling region, neutrino cooling makes the slope steeper. TD and luminosity terms ultimately originate from the energy Equation (3). The turbulent ram pressure term comes from the momentum equation. While turbulent ram pressure does modify the density structure, the two turbulent terms from the energy equation, TD and turbulent luminosity, have a considerably larger effect on the structure.

Zoom In Zoom Out Reset image size

Figure 2 shows how turbulence affects the density and temperature profiles of the steady-state stalled-shock solutions. The net effect of turbulence is to make the profiles much shallower. In part, turbulent ram pressure explains some of the difference, but for the most part, turbulence provides extra heating through dissipation in the convective region. One consequence is that the temperature (and entropy) are higher with turbulence. This is consistent with the entropy profiles of multidimensional simulations (Murphy & Meakin 2011; Murphy et al. 2013). In agreement with simulations (Murphy et al. 2013; Abdikamalov et al. 2016), Figure 2 also shows that the solution including turbulence has a larger shock radius. The shock radius is a monotonic function of L_ν, and since $\rho {R}_{{rr}}$ and $\rho \epsilon$ effectively add energy in the same fashion as the luminosity, we intuitively retrieve a larger shock radius. Note that this larger shock radius is not just a consequence of turbulent ram pressure. Instead, the post-shock profile is shallower, pushing out the shock.

Zoom In Zoom Out Reset image size

Now that we understand how turbulence affects the density and temperature profiles, we now present how turbulence affects the critical curve in three figures. One, Figure 3 shows how turbulence raises the dimensionless overpressure parameter, Ψ, in Murphy & Dolence (2017). Two, Figure 5 shows how turbulence reduces the five-dimensional critical hypersurface for explosion Finally, Figure 7 shows that the dominant turbulent term in reducing the critical condition is TD.

Zoom In Zoom Out Reset image size

In Figure 3, we plot Ψ to show that the increase in the dimensionless parameter due to turbulence is roughly equivalent to increasing the neutrino luminosity by 30%. To clarify, the Ψ parameter is a measure of the overpressure compared to the hydrostatic equilibrium. This integral condition is normalized by the pre-shock ram pressure. This dimensionless parameter maps directly to the shock velocity in that when Ψ is positive, the shock expands; when Ψ is negative, the shock stalls; and when Ψ = 0, the shock has zero velocity. Note that there is always a minimum Ψ, and if this ${{\rm{\Psi }}}_{\min }$ is greater than zero, then there are no stalled-shock solutions, only steadily expanding shocks. For the case where the minimum of Ψ is exactly zero, this set of solutions defines our critical explosion condition for all parameters. Clearly, turbulence raises the minimum Ψ, and therefore would affect the critical curve.

In Figure 4, we show how the Reynolds stress upper limit affects the explodability parameter and the critical curve. We suspected that, at high-enough neutrino luminosities, the additional ram pressure at the shock would prevent finding solutions to the boundary conditions. Figure 4 demonstrates that we consistently encounter this cap, but only for pseudo-solutions that have already found a critical L_ν. We present several sets of solutions at different neutrino luminosities to emphasize that our ${L}_{\nu }-\dot{{ \mathcal M }}$ critical curve is in fact not affected by the upper limit on ${ \mathcal R }$, even at unrealistic luminosities. The sole determination of the critical curve is on ${{\rm{\Psi }}}_{\min }={\rm{\Psi }}({r}_{{\rm{s}}}^{\mathrm{crit}})$, and since this threshold is only reached when ${r}_{{\rm{s}}}\gt {r}_{{\rm{s}}}^{\mathrm{crit}}$, there is no effect on the critical curve.

Zoom In Zoom Out Reset image size

Figure 5 shows how turbulence modifies the critical hypersurface. Murphy & Dolence (2017) point out that the critical condition for explosion is not a critical curve, but a critical hypersurface in a five-dimensional space that is defined by one dimensionless condition: ${{\rm{\Psi }}}_{\min }=0$. The neutrino-luminosity–accretion-rate curve is one slice of this critical hypersurface. In Figure 5, we show six slices of this hypersurface. Note that in all panels except the top-left panel (${M}_{\mathrm{NS}}$ versus $\dot{{ \mathcal M }}$), the region of explosion is above the curve. For ${M}_{\mathrm{NS}}$ versus $\dot{{ \mathcal M }}$, it is below the curve; a lower-mass neutron star has a lower potential to overcome to explode. Because turbulence raises the Ψ curves, it also reduces the critical hypersurface for explosion.

Zoom In Zoom Out Reset image size

In Figure 6, we make a case for the validity of using a convection-based turbulence model. The χ parameter, first introduced by Foglizzo et al. (2006), is a measure of the linear stability of the convective region in the presence of advection. For $\chi \lt 3$, advection stabilizes the flow, and convection does not manifest. The assumption is that since convection is suppressed, SASI dominates the turbulence. The χ parameter is defined as

$$\begin{eqnarray}&&\chi \equiv {\int }_{{r}_{g}}^{{r}_{s}}\left|\displaystyle \frac{{\omega }_{\mathrm{buoy}}}{u}\right|\,{dr},\end{eqnarray} \tag{ 41 }$$

where ${\omega }_{\mathrm{buoy}}$ is the Brunt–Väisälä frequency, defined as

$$\begin{eqnarray}&&{\omega }_{\mathrm{buoy}}^{2}\equiv \left[\displaystyle \frac{1}{p}{\left(\displaystyle \frac{\partial p}{\partial S}\right)}_{\rho ,{Y}_{e}}\displaystyle \frac{{dS}}{{dr}}+\displaystyle \frac{1}{p}{\left(\displaystyle \frac{\partial p}{\partial {Y}_{e}}\right)}_{\rho ,S}\displaystyle \frac{{{dY}}_{e}}{{dr}}\right]\displaystyle \frac{g}{{{\rm{\Gamma }}}_{1}}\,\end{eqnarray} \tag{ 42 }$$

and

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1}\equiv {\left(\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}\rho }\right)}_{S,{Y}_{e}}.\end{eqnarray} \tag{ 43 }$$

In Figure 6, above the χ = 3 dashed line, all stalled solutions have χ > 3 (convection dominated), and below this dashed line, all stalled solutions have χ < 3 (SASI dominated). To calculate this line, we use an approach similar to the derivation of the critical luminosity curves: we input all of our parameters, and numerically solve for the luminosity, which gives a value of χ = 3. According to these models, convection dominates for all scenarios near explosion. Recent simulations seem to support this conclusion (Burrows et al. 2012; Ott et al. 2013; Couch & O'Connor 2014; Iwakami et al. 2014; Abdikamalov et al. 2015; Lentz et al. 2015; Roberts et al. 2016). The results of Figure 6 validate our exploration of how convection-based turbulence affects the critical condition. However, we do encourage future explorations of SASI and χ in simulations to validate whether convection is indeed dominant near the critical condition.

Zoom In Zoom Out Reset image size

In Figures 3 and 5, we considered the overall effect of turbulence on the critical condition, but this does not illuminate how turbulence reduces the condition for explosion. In Figure 7, we focus exclusively on the ${L}_{\nu }-\dot{{ \mathcal M }}$ critical curve and explore the effect of each term in this slice of the critical condition. We find that TD within the gain region acts as an even greater driving force for explosion than the turbulent ram pressure. That is not to say that the Reynolds stress is negligible; both terms are indeed needed to obtain the critical curve reduction predicted by multidimensional simulations. Although this result relies upon some assumptions in the turbulence model, we suspect that the qualitative outcome will persist: TD cannot be dismissed.

Zoom In Zoom Out Reset image size

3.1. Buoyancy-driven Heating

We argue that TD adds another heat source, which aids explosion. Since neutrinos are the ultimate source of energy, it may seem that we are asking neutrinos to do twice the work. However, this is not the case. Using a simple convective model, we propose that in one dimension, some of the neutrino energy goes into heating the material, and some of it goes into creating a higher-potential profile. In multiple dimensions, this higher-potential profile is unstable and goes to a lower-potential profile, and the excess energy goes into kinetic energy, which in turn dissipates as heat via TD.

Rayleigh–Bénard convection is a simple convective model that can clearly demonstrate this conversion from potential to kinetic to dissipated internal energy. Figure 8 illustrates the fundamental physics of convection. First, neutrinos provide a source of heating and drive a convective instability. In this cartoon model, we consider two parcels; the lower one receives more neutrino heating, and has a higher entropy and lower density compared to its surroundings. The parcel at a greater height has a lower entropy and higher density compared to its surroundings. If one switches the positions of these two parcels, then one finds that the gravitational potential is lower. Therefore, neutrino heating causes a higher-potential structure that is unstable to convective overturn. The difference in potentials between the two states gets converted into the kinetic energy of the parcels. Kolmogorov's hypothesis suggests that the dissipation of this kinetic energy is ${R}^{3/2}/{ \mathcal L }$ and happens on the order of one turnover timescale. Burrows et al. (1995) also considered this idea in which two layers of varying densities are swapped, inducing a buoyant work being done on the system. This energy is then converted into kinetic energy in the form of eddies, and is in turn dissipated into heat. Therefore, not only do neutrinos heat the gain region, they also create a higher-potential system. This higher potential gets converted to kinetic energy and is consequently dissipated as internal energy.

Zoom In Zoom Out Reset image size

To illustrate this more clearly, consider the energy Equation (3). It is more illuminating if we rewrite the equation considering a constant mass-accretion rate,

$$\begin{eqnarray}&&\dot{{ \mathcal M }}{\partial }_{i}\left[h+\displaystyle \frac{{u}^{j}{u}_{j}}{2}-{\rm{\Phi }}\right]=\rho q.\end{eqnarray} \tag{ 44 }$$

Neutrinos heat the convective region of the star, changing the enthalpy and gravitational potential. This gives rise to an entropy and density gradient such that ${S}_{2}\gt {S}_{1}$ and ${\rho }_{2}\lt {\rho }_{1}$ (thus satisfying the Schwarzschild condition for convection). We then treat the potential energy of two parcels in a similar manner to Burrows et al. (1995). Before the exchange of parcels, the gravitational potential energy is

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{1}={gV}({\rho }_{1}{h}_{1}+{\rho }_{2}{h}_{2}).\end{eqnarray} \tag{ 45 }$$

But after the top parcel sinks and the bottom parcel buoyantly rises, the potential is then

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{2}={gV}({\rho }_{2}{h}_{1}+{\rho }_{1}{h}_{2}).\end{eqnarray} \tag{ 46 }$$

Thus, the change in potential is

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Phi }}={gV}({\rho }_{2}-{\rho }_{1})({h}_{1}-{h}_{2}).\end{eqnarray} \tag{ 47 }$$

Since ${h}_{1}\gt {h}_{2}$ and ${\rho }_{1}\gt {\rho }_{2}$, ΔΦ < 0. Conservation of energy suggests that K.E. ≈ −ΔΦ, and via TD, this energy from buoyant driving is converted into heat. In summary, neutrinos heat the gain region and set up a higher-potential profile. Turbulence allows a lower-potential state, and the kinetic energy is converted into thermal energy via TD, aiding explosion.

A major result of core-collapse theory is that one-dimensional simulations fizzle for all but the least massive stars, while multidimensional simulations explode for the most part. Even though there are still some differences between the simulations, there is a general consensus that turbulence spells the difference between a fizzled outcome and a successful explosion. To explain how turbulence aids explosion, we develop a turbulence model for neutrino-driven convection (Murphy et al. 2013) and investigate how turbulence reduces the critical condition for explosion (Burrows & Goshy 1993; Murphy & Dolence 2017). Our turbulence model reduces the critical condition for explosion by ∼30%, in general agreement with simulations. By modeling each turbulent term, we are able to investigate the effect of each turbulent term on the critical condition. We find that although ram pressure plays a role in reducing the critical curve, it is not the dominant term. The dominant term in reducing the critical curve is TD. Furthermore, we are not the first to suggest that TD is important in aiding neutrino-driven explosions. Thompson et al. (2005) suggested that MRI-driven turbulence for very rapidly rotating proto-neutron stars could add significant heating and aid explosion.

In the turbulence model, we include all turbulent terms, both in the background solution and in the boundary conditions. The three main turbulent terms are TD, Reynolds stress, and turbulent luminosity. Overall, we find that all three play an important role in modifying the background structure and the critical condition for explosion. However, it is the terms in the energy equation, the combined turbulent luminosity and TD, that give the largest reduction in the critical curve. Of these two, the TD provides the largest effect in reducing the critical condition.

The ultimate source of power for convection and TD is the neutrino-driving power. This may seem as if we are double-counting the power supplied by neutrinos. However, we are not. Instead, we propose that neutrinos heat the gain region and set up a higher-potential, convectively unstable structure. In multidimensional simulations, convection converts this higher-potential structure to a lower-potential structure. The change in potential energy is converted to kinetic turbulent energy, which in turn dissipates as heat. We suggest that CCSN modelers check the structures of their one-dimensional and multidimensional simulations to test this supposition.

The turbulence model is a global, integral model, and provides little constraint on the local structure of turbulence, but a local model is necessary in solving the background equations and deriving a critical condition. Therefore, we made some fairly straightforward assumptions to translate the global model into a local model. For example, we assumed that the specific TD rate is constant throughout the gain region. We also had to assume a specific spatial profile for the turbulent luminosity. These assumptions probably do not affect our qualitative results. However, our results have large implications for how turbulence aids explosion. In particular, we propose that TD is a key contributor to reviving a stalled shock to a successful supernova. Therefore, these assumptions should be verified with multidimensional simulations and treated more rigorously in future investigations.

To verify the predictions of this manuscript, we identify at least three open questions that multidimensional simulations should address. One, does TD actually lead to significant heating? In multiple dimensions, the entropy profile should be a result of neutrino heating and cooling, turbulent entropy flux, and TD. One can easily calculate the turbulent entropy flux, and the heating and cooling by neutrinos. Whatever is left should be equal to the expected entropy generated by TD. Second, do the one-dimensional profiles have a higher potential compared to their multidimensional counterparts? Third, do the local details of the turbulent model matter? We suspect that the local details do not change the qualitative result. However, one should compare our local assumptions with multidimensional simulations, and assess how (if at all) the quantitative results vary from this work.

In summary, combining a turbulence model and critical condition analyses helps to illuminate how turbulence aids explosions. Specifically, by modeling each turbulent term, we are able to assess the effect of each term on the conditions for explosion. Contrary to prior suppositions, we find that turbulent ram pressure is not the dominant effect. Rather, each of the three terms have quite large effects, with TD being the largest. At present, these conclusions are qualitative. To be predictive, the community will need to verify these conclusions with multidimensional simulations. Eventually, we may be able to use these turbulence models to make one-dimensional simulations explode under conditions similar to multidimensional simulations, thereby enabling rapid and systematic exploration of the explosion of massive stars.

We would like to thank Joshua C. Dolence for illuminating discussions on the effect of turbulence on the critical condition. We are also thankful to David Radice for challenging us to explain the source of energy for turbulent dissipation. This material is based upon work supported by the National Science Foundation under grant No. 1313036.

Here, we present the full set of ordinary differential equations that we solve to find pseudo steady-state solutions. The equation for density is

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}\rho }{d\mathrm{ln}r}\\ &&=\,\displaystyle \frac{\tfrac{2{v}^{2}}{{\rm{\Phi }}}-1+{ \mathcal R }{ \mathcal S }+\tfrac{y}{l}\tfrac{\tfrac{\partial \mathrm{ln}y}{\partial \mathrm{ln}T}}{\tfrac{\partial \mathrm{ln}l}{\partial \mathrm{ln}T}}\left(1+{ \mathcal H }-{ \mathcal C }+{ \mathcal T }{ \mathcal D }-{ \mathcal T }{ \mathcal L }+\tfrac{2{v}^{2}}{{\rm{\Phi }}}\right)}{y-\tfrac{{v}^{2}}{{\rm{\Phi }}}+y\tfrac{\partial \mathrm{ln}y}{\partial \mathrm{ln}\rho }+\tfrac{\tfrac{\partial \mathrm{ln}y}{\partial \mathrm{ln}T}}{\tfrac{\partial \mathrm{ln}l}{\partial \mathrm{ln}T}}\left(\tfrac{{{yv}}^{2}}{{\rm{\Phi }}l}-\tfrac{\partial \mathrm{ln}l}{\partial \mathrm{ln}\rho }\right)},\end{eqnarray} \tag{ 48 }$$

and in terms of the density equation, the temperature ODE is

$$\begin{eqnarray}\displaystyle \frac{\partial \mathrm{ln}T}{\partial \mathrm{ln}r} & = & \displaystyle \frac{\left(1+{ \mathcal H }-{ \mathcal C }+{ \mathcal T }{ \mathcal D }-{ \mathcal T }{ \mathcal L }+\tfrac{2{v}^{2}}{{\rm{\Phi }}}\right)}{l\tfrac{\partial \mathrm{ln}l}{\partial \mathrm{ln}T}}\\ & & +\tfrac{\tfrac{\partial \mathrm{ln}\rho }{\partial \mathrm{ln}r}}{\tfrac{\partial \mathrm{ln}l}{\partial \mathrm{ln}T}}\left(\tfrac{{v}^{2}}{{\rm{\Phi }}l}-\tfrac{\partial \mathrm{ln}l}{\partial \mathrm{ln}\rho }\right).\end{eqnarray} \tag{ 49 }$$

In this form, these equations seem somewhat unwieldy, but they would be even more so if we had not made the following shorthand for the important physics. To further help illustrate the important scales in the problem, we present each important physics in terms of a dimensionless variable. The Reynolds stress (or ram pressure) appears in the above equations as

$$\begin{eqnarray}&&{ \mathcal R }{ \mathcal S }=\displaystyle \frac{2{x}^{2}{R}_{{rr}}}{({x}_{s}-{x}_{g}){{\rm{\Phi }}}_{1}}.\end{eqnarray} \tag{ 50 }$$

In these expressions, the subscript 1 indicates the base of the solution. In our particular case, that is the neutrinosphere radius. Neutrino heating and cooling are

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \frac{{L}_{\nu }\kappa \rho {r}_{1}{x}^{2}}{\dot{{ \mathcal M }}{{\rm{\Phi }}}_{1}},\end{eqnarray} \tag{ 51 }$$

and

$$\begin{eqnarray}&&{ \mathcal C }=\displaystyle \frac{{C}_{0}{\left(\tfrac{T}{{T}_{0}}\right)}^{6}4\pi {r}_{1}^{3}{x}^{4}\rho }{\dot{{ \mathcal M }}{{\rm{\Phi }}}_{1}}.\end{eqnarray} \tag{ 52 }$$

The two turbulent terms from the energy equation are the TD and turbulent luminosity:

$$\begin{eqnarray}&&{ \mathcal T }{ \mathcal D }=\displaystyle \frac{{L}_{e}^{\max }{r}_{1}^{3}\rho {w}_{b}4\pi }{{M}_{g}\dot{{ \mathcal M }}{{\rm{\Phi }}}_{1}}\end{eqnarray} \tag{ 53 }$$

and

$$\begin{eqnarray}&&{ \mathcal T }{ \mathcal L }=\displaystyle \frac{{L}_{e}^{\max }{r}_{1}{x}^{2}}{\cos {h}^{2}\left(\tfrac{({r}_{1}(x-{x}_{g}))}{h}\right)h\dot{{ \mathcal M }}{{\rm{\Phi }}}_{1}}.\end{eqnarray} \tag{ 54 }$$

Two dimensionless measures of the pressure and enthalpy are

$$\begin{eqnarray}&&y=\displaystyle \frac{P}{\rho {\rm{\Phi }}}\end{eqnarray} \tag{ 55 }$$

and

$$\begin{eqnarray}&&l=\displaystyle \frac{e+\tfrac{P}{\rho }}{{\rm{\Phi }}},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{\left.{P}_{,\rho }=\displaystyle \frac{\partial \mathrm{ln}P}{\partial \mathrm{ln}\rho }\right|}_{T}.\end{eqnarray} \tag{ 57 }$$

The normalized radius is

$$\begin{eqnarray}&&x=\displaystyle \frac{r}{{r}_{\mathrm{NS}}}.\end{eqnarray} \tag{ 58 }$$

A dimensionless measure of the buoyant driving is

$$\begin{eqnarray}&&{w}_{b}={\int }_{{x}_{g}}^{{x}_{s}}\displaystyle \frac{\tan h\left(\tfrac{{r}_{1}(x-{x}_{g})}{h}\right)}{3{yx}}{dx}.\end{eqnarray} \tag{ 59 }$$

To observe and distinguish the effects of each turbulent term, we add the capability to turn each term on or off in the equations. In doing so, we are able to investigate the effect of each term on the solutions and critical curves.