On the treatment of phenomenological turbulent effects in one dimensional simulations of core-collapse supernovae

We have developed a phenomenological turbulent model with one-dimensional (1D) simulation based on Reynolds decomposition. Using this method, we have systematically studied models with different effects of compression, mixing length parameters, and diffusion coefficient of internal energy, turbulence energy and electron fraction. With employed turbulent effects, supernova explosion can be achieved in 1D geometry, which can mimic the evolution of shock in the 3D simulations. We found that enhancement of turbulent energy by compression affects the early shock evolution. The diffusion coefficients of internal energy and turbulent energy also affect the explodability. The smaller diffusion makes the shock revival faster. Our comparison between the two reveals that the diffusion coefficients of internal energy has a greater impact. These simulations would help understand the role of turbulence in core-collapse supernovae.


INTRODUCTION
Massive stars undergo core-collapse supernova explosions in the final stage of their evolution.The energy source of the explosion is the enormous gravitational potential energy released from the implosion by their self-gravity.Despite intensive and extensive studies in the past, the mechanism of the explosion has not yet been fully elucidated (see Kotake et al. 2012;Burrows 2013;Foglizzo et al. 2015;Zhou et al. 2019;Burrows & Vartanyan 2021;Müller 2020, for the recent progress).Complex phenomena such as neutrino radiative transport and convection play an essential role in the explosion, and their numerical simulations require extremely sophisticated schemes.
Thanks to the enlarged computational resources, recently, selfconsistent multi-dimensional simulations have been performed.Some of the models show the shock-revivals and explosions since neutrino-driven convection substantially enhances the neutrino heating (Takiwaki et al. 2012;Hanke et al. 2013;Murphy et al. 2013;Abdikamalov et al. 2015;Lentz et al. 2015;Summa et al. 2016;Müller et al. 2016;O'Connor & Couch 2018;Pan et al. 2018;Kuroda et al. 2018;Vartanyan et al. 2019;Nagakura et al. 2019a;Melson et al. 2020;Burrows et al. 2020).Although those simulations are outstanding achievements of this field, we have to point out that most of them have not been able to reproduce the typical explosion energy or the amount of synthesized 56 Ni (Murphy et al. 2019;Suwa et al. 2019, however, see Bruenn et al. 2013;Nakamura et al. 2016;Bollig et al. 2021 for a few exceptions).Still, we have not captured the nature of core-collapse supernovae, and we have missed something important in the simulations.
★ E-mail: sasaki.shunsuke.astro@gmail.comEncouraged by the partial success of the 3D simulations, a variety of phenomenological 1D approaches have been used to mimic the 3D models and have linked the character of the progenitor and the explodability (O'Connor & Ott 2011; Ugliano et al. 2012;Ertl et al. 2016;Sukhbold et al. 2016;Ebinger et al. 2020).The authors enlarge the neutrino luminosity and hence neutrino-heating to obtain the successful explosions since spherically symmetric 1D simulations cannot revive shocks (e.g., Sumiyoshi et al. 2005).The studies found that several parameters govern the explodability: e.g., compactness,   ;  4 and  4 1 .Those approaches are useful to compare the distribution of the observed properties (Horiuchi et al. 2014;Pejcha & Thompson 2015) and to estimate diffusive supernova neutrino background integrating all neutrino flux from the past supernovae (Lunardini & Tamborra 2012;Horiuchi et al. 2018Horiuchi et al. , 2021;;Kresse et al. 2021;Ashida & Nakazato 2022;Ashida et al. 2023;Ekanger et al. 2022).
One current hot topic in this field is a discrepancy between the prediction in the phenomenological and self-consistent models.Burrows et al. (2020) have performed that 3D simulations, changing the progenitor mass systematically.That and the following works (Wang et al. 2022;Tsang et al. 2022) claim that the compactness does not govern explodability (see also Nakamura et al. 2015;Summa et al. 2016;O'Connor & Couch 2018, for 2D studies).The previous phenomenological studies just enlarged the neutrino luminosities by hand.Though this treatment would be motivated by the fact that turbulent effects make a longer residency timescale in the gain layer (Murphy & Burrows 2008;Hanke et al. 2012;Nordhaus et al. 2010), the convective effects of the self-consistent simulations would not be the same as the phenomenological one.
To solve the discrepancy between the self-consistent and phenomenological simulations, we need to theoretically understand the turbulent effect and make more sophisticated phenomenological models.Here we introduce efforts to understand the turbulent convection in supernovae.Murphy & Meakin (2011) have derived and summarized the equations for convective effects in this context.Yamasaki & Yamada (2006) have estimated how the shock revival becomes easier when employing strong energy transport due to convection.Couch & Ott (2015) and Radice et al. (2016) have concluded that the turbulent pressure pushes the shock outward and causes the explosion.On the other hand, Mabanta & Murphy (2018) have claimed that the increased thermal pressure by the dissipation of turbulence makes a more substantial effect than the turbulent pressure.Kazeroni et al. (2018) have taken into account the size of the initial perturbations and the effect of the drag force.Müller (2019) has focused on the effect of turbulent viscosity.We are under the middle way to the complete understanding of convective turbulence in core-collapse supernovae.We should learn a lot from numerical studies in stellar evolution (e.g., Arnett et al. 2015).
Very recently, time-dependent simulations with phenomenological convective turbulence have been performed (Mabanta et al. 2019;Couch et al. 2020;Müller 2019;Boccioli et al. 2021Boccioli et al. , 2022Boccioli et al. , 2023)).Those try to model the convective turbulence in multi-dimensional simulations, and install the effects in one-dimensional simulations with several parameters, while the previous phenomenological studies just increased neutrino luminosity to trigger the explosion.Following Mabanta et al. (2019), we refer this method as 1D+2 .Note that the Mabanta et al. (2019) try to establish global modeling of the convective turbulence where the turbulent properties are not determined by local hydrodynamic variables.On the other hand, other studies formulate it by local modeling (Couch et al. 2020;Müller 2019;Boccioli et al. 2021).In the following, we basically focus on the local modeling, however, we must stress the importance of the global modeling and further investigation on this topic.
These approaches not only would reproduce the dynamics close to the results of 3D simulations, but also are useful to determine the progenitor dependence of the explodability (Couch et al. 2020;Boccioli et al. 2021;Zapartas et al. 2021), taking advantage of its low computational cost as in other phenomenological approaches.Predicting the neutrino luminosity, the frequency of the gravitational waves, and the optical light curves, 1D+ is also used in the context of multi-messenger observation of core-collapse supernova (Warren et al. 2020;Barker et al. 2022).
The purpose of this study is to understand the effects of turbulence in a current 1D+ simulation.In particular, we perform a parameter survey to mimic the evolution of shock in a 3D simulation using the same progenitor model.We show that compression affects prompt convection in our 1D+ simulations.We also analyze the effect of turbulence parameters and show that the effect of diffusion is different from the effect proposed in previous studies.The weaker the diffusion, the more likely the explosion.Our simulation method is introduced in Section 2 and Appendix A. In Section 3.1, we show that the treatment of compression of turbulent energy has impacts on convection by comparing our results with a 3D simulation.In Section 3.2, we show the dependence of mixing length parameter on the shock evokution comparing it to 3D simulation.In Section 3.3, we show the dependence of diffusion parameters on supernova explodability.

METHOD
We employ the neutrino radiation hydrodynamic code of 3DnSNe-IDSA (Takiwaki et al. 2016) and have newly implemented phenomenological convection effects to incorporate convection effects into 1D spherically symmetric simulations (e.g., Mabanta et al. 2019;Müller 2019;Couch et al. 2020).Based on Reynolds decomposition, the governing equations are written as follows: where ρ, v, êtot , P, Ŷ , g are density, radial velocity, total energy, pressure, electron fraction, and gravitational acceleration, respectively.We follow the notation of Müller (2019), and Â or ⟨⟩ mean angle-average of the variable, .Prime is used for the deviation from the mean value.The radial divergence, 1  2    2 , is simply written as ∇  .Here êtot = ρ ε + ρ v2 /2 +  turb , where ε is specific internal energy and  turb is turbulent energy.In the first term of RHS in Eq. (2), sum of  and  components of stress tensor appears3 (e.g.see Table B.1 in González et al. 2007).The term comes from the divergence of tensor in the spherical coordinate.Eq. ( 130) in Stone & Norman (1992) denotes the derivation.
We explain the difference between the turbulence models of previous studies and ours, where Here  is the turbulent components of the Reynolds tensor.The constant  determines how much fraction of the turbulent tensor contributes to the radial-radial component of turbulent pressure tensor.The relation between  turb and  turb is  turb = 2 turb .The pressure due to turbulence is assumed as  turb =  turb with reference to Couch et al. (2020) and Müller (2019).Couch et al. (2020) adopted the same closure of Mabanta et al. (2019) as Reynolds tensor, which by our definition is   = 2   = 2   .We assume those for turbulent pressure and set  = 0.5 in order to be consistent with notation of Couch et al. (2020) and Müller (2019).Similar to the pressure tensor, we also have turbulent pressure in RHS of Eq. ( 2), .This term is ignored in the previous studies 4 .
The effects of neutrino interaction are written as   , Γ  , and   .We calculate   , Γ  , using Isotropic Diffusion Source Approximation (IDSA, Liebendörfer et al. 2009;Takiwaki et al. 2014;Kotake et al. 2018).  is ignored in the approximation since the mass of nucleons are high compared to the energy of the neutrinos, and the momentum of the matter is hardly changed by the neutrino scattering.On the other hand, the effect of the collision is included in the neutrino transport.
This paper argues the impact of  HT in Eq.(4) (HT is abbreviation of Hydrodynamic enhancement of the Turbulent energy).Müller (2019) and Couch et al. (2020) use  HT =  turb 1  2    2 v and  HT =  turb   v, respectively.Though those looks similar, the meanings are slightly different.The former indicates the enhancement of the turbulent energy is due to the compression and the latter indicates that the shear is important.The origin of this term is actually the shear (see Eq. ( 7) of Murphy & Meakin (2011), the term is called shear production) 5 .However, the total energy conservation is not achieved by the latter expression and the former expression ensures the energy conservation (see Appendix A).Suming up  • ∇ turb and  turb ∇ •  in Eq. ( A8) and (A9), we obtain ∇ • ( turb ) in Eq. ( 3).If  HT ≠  turb 1  2    2 v, consequently Eq. ( 3) should be acoordingly modified.In this paper, we tolerate this inconsistency.We need to compare these pros and cons employing this term.See Section 3.1 for the detail.
We estimate multi-dimensional turbulent effect using mixing length theory as follows.The effect is introduced from 10 ms after core bounce.At first, the mixing length Λ is where   is scale height of pressure and  Λ is dimensionless parameter, which is also used in Couch et al. (2020) and Müller (2019).
In the governing equations,   ,  = ,   ,  are the turbulent flux of internal energy, turbulent energy and electron fraction,   , respectively.To evaluate them, we take the gradient diffusion approximation.In the mixing length theory, the diffusion coefficients are defined by the mixing length, Λ and the turbulent velocity,  turb : where   is the parameter for each of the three diffusion coefficients.Approximately, the turbulent velocity is evaluated as The three fluxes due to turbulence can be written in terms of their 4 In Couch et al. (2020), (   +    )/ =  turb / is ignored in Eq. ( 26).
The source term of the turbulent energy in Eq. ( 4) is expressed as where  BV is the Brunt-Väisälä frequency, which is defined as: where  L represents the density difference from the surroundings when the fluid element is adiabatically displaced by d.g eff is effective gravitational acceleration, This term represents the apparent gravitational acceleration in a fluid stationary system mentioned in Eq.( 16) of (Couch et al. 2020).In our notation,  2 BV takes a positive value in the convectively unstable region where the entropy, , or   gradients are negative6 .Since small turbulent seed is required to turbulent energy source, we assume that the minimum value of turbulent energy source following Eq.( 23) in Couch et al. (2020) is Note that Müller (2019) use a different handling, see their Eq.( 47).
In this study, we focus on the role of turbulence in the gain region and switch off this term at ρ > 10 11 g/cm 3 not to include the convection in PNS.In order to suppress turbulence generation at the shock discontinuity,  BW is set to zero near the shock by hand.Note that this term naturally appears in Eq. ( 3).Since we ignore ⟨ ′  ′ ⟩ in Eq. ( 1), the treatment is not consistent and is criticised by Müller (2019).The total energy including gravitational potential energy, does not conserve in this formalism.However, as discussed in Appendix A, the treatment of Müller (2019) invokes another problem.In addition, Mabanta and Murphy have stated in private communication that the addition of this term will not have a significant impact.We prefer to stick the current formalism.
Eq. ( 14) is a major assumption of our study.Murphy et al. (2013) claim that the generation term is proportional to the neutrino luminosity, not the negative entropy gradient itself, which is significantly altered by turbulent mixing.Since our simulations solve timedependent fluid dynamics with neutrino radiative transport, we adopt Eq. ( 14) for simplicity as a first step.This treatment should be updated in the future.We discuss this issue again in Section 3.3.
The turbulent energy is injected at the largest scales and cascades to smaller scales.Consequently, the largest scale of the eddy, Λ diss , governs the rate of the dissipation: Following Müller (2019), we incorporates the effect of the overshooting: where  = 10 −10 cm.
The details of numerical treatment on the governing equations are written in O' Connor et al. (2018).The numerical flux is calculated by a HLLC solver (Toro et al. 1994).A piecewise linear method with the geometrical correction of Mignone ( 2014) is used to reconstruct variables at the cell edge, where a modified van Leer limiter is employed to satisfy the condition of total variation diminishing (TVD).The computational grid is comprised of 512 logarithmically spaced, radial zones that cover from the center up to the outer boundary of 5 × 10 8 cm.The effect of the GR potential is included by Case A of Marek et al. (2006).
The setup for the microphysics is the same to set-all of Kotake et al. ( 2018), which employ the latest neutrino reaction rate, e.g., the medium effect is included (Horowitz et al. 2017).We use the EOS of Lattimer & Swesty (1991)

RESULTS
This section presents the outcomes of hydrodynamic simulations, which utilize two turbulence models.The first model, the noncompressive case, ignores the enhancement of turbulent energy due to compression, while the second model, the compressive case, takes this enhancement into account.We systematically explore the impact of changes in the mixing length parameter,  Λ , and the three diffusion coefficients,   ,   , and    .For all models, 12 ⊙ progenitor model of Woosley & Heger (2007) is used as the initial condition.
The initial part of the section examines the effect of compression on turbulence, displaying the contrasts in the shock evolution between the non-compressive and compressive cases.We also present the time-dependent variation of the source term of turbulence energy and the distribution of turbulent velocity to illustrate the influence of compression in our turbulence models.
The subsequent section employs the non-compressive models and showcases the explosion dynamics while altering the mixing length parameter.We exhibit the differences in shock evolution for varying mixing length parameters and present the variations in the velocity distribution of turbulence to emphasize the impact of this parameter.
Finally, the last section illustrates how the dynamics change as a result of diffusion parameters.Analogous to the previous section, we present the differences in shock evolution and velocity distribution of turbulence to demonstrate the effects of diffusion parameters.

Compressive enhancement of turbulence
In Section 2, we introduced several types of terms that appear in the turbulent energy equation, Eq.( 4).Here we just write it again for the readers' convenience.
We compare two turbulent models, non-compressive case and compressive case.For the non-compressive case, we assume  HT = 0, while for the compressive case, we use  HT =  turb 1  2    2 v.We implement this term as follows, where the subscript  means -th grid.The radial velocity, v, should be evaluated in the cell edge (represented by  + 1/2 or  − 1/2).Since we use 2nd-order interpolation to obtain the numerical flux, we utilize the interpolated value.The subscript  and  denotes the direction of the interpolation.Figure 1 displays the evolution of the shock for the model without compression (non-compressive case) and with compression (compressive case) .The vertical axis is the shock radius and the horizontal axis represents the time after core bounce.The solid line corresponds to the non-compressive case and the dashed line does compressive case.The black line denotes the time evolution of the average shock radius for the 3D simulation using the same progenitor.The angular grid number of the 3D model is 64 × 128 (see Takiwaki et al. 2016, for the numerical setup.The emplyoed progenitor is different from this paper).In the 3D model, the turbulence naturally occurs, and the phenomenological effect is not included.We define the shock radius based on Eq. ( 72) of Martí & Müller (1996), i.e., a discontinuous surface where the pressure difference is greater than the value of the pressure in the neighborhood, and the velocity is converging.The value of entropy is also used to make it more smooth in time.
Based on the results, we can conclude that the non-compressive case (red solid) is closer to the evolution of the 3D model (black bold) than the compressive case (red dashed).Compressive case shows different behavior from the 3D model from 50 ms, and radius of the shock is significantly deviates from the 3D model as time passes; for instance, 50 km of the difference appears in 160 ms.The turbulent energy is considerably increased by the compression term.To reproduce the early shock trajectory, accurately, it is preferable to ignore the effect of  HT .One may expect that the compressive case with smaller  Λ can fit to the 3D model.If we use the parameter that makes the shock radius at 100 ms consistent with 3D, we cannot find shock revival and the later evolution is not consistent with 3D (see Figure 4 for the evolution in the later epoch).
To understand the detail of the two models, we compare the generation and dissipation term of the turbulent energy.Figure 2 shows neutrino heating rates and the source terms in the governing equation for turbulent energy in the gain regime.The vertical axis is source terms in the unit of 10 51 erg/s and the horizontal axis is postbounce time in ms.The top panel and bottom panels represent non-compressive case, and compressive case, respectively.Orange, blue, and green indicate the turbulent energy source rate, turbulent energy dissipation rate, and turbulent energy increase rate due to turbulence compression in the gain region, respectively.The definition of each term is as follows In the figure the red curve is the integrated neutrino heating rate, i.e., In the governing equations,  turb gain increase the turbulent energy and  dis gain decrease the energy.We expect that  turb gain is larger than the  dis gain before the system becomes equilibrium.In the top panel, we confirm  turb gain >  dis gain .Murphy et al. (2013) reported that the turbulence in the gain region is maintained by the neutrino heating and   >  turb gain is also expected.In the top panel,   >  turb gain is found before 170 ms and those terms remains comparable after that time.These features are roughly consistent with Murphy et al. (2013).Here we add a caveat.In our model, the source term,  turb gain indirectly related to the neutrino heating.The neutrino heating is stronger at the bottom of the gain region and is weaker near the shock.This makes negative entropy gradient and that determines Brunt-Väisälä frequency, see Eqs. ( 15)-( 16).Essentially, the neutrino heating accounts for the generation of the turbulent energy.Our treatment is valid in the liner phase of turbulent growth.In the nonlinear phase, the entropy gradients is flatten by the turbulence and the our treatment may underestimate the generation term.To avoid the problem, the source term in Mabanta et al. (2019) is directory related to the neutrino heating.We may try to implement this method in the future.
In contrast to the non-compressive case, the compressive case poses several issues.In the bottom panel of Figure 2, the contribution of source terms in the gain region are shown in the compressive case.First, it is unlikely that the dissipation term,  dis gain , is larger than the source term,  turb gain in all time, 0-200 ms.Instead of  turb gain , the compression term  gain plays dominant production role in the initial phase, 0-50 ms.Here the turbulent energy of the seed turbulence is amplified by compression near the shock.In this phase, the energy gain is much larger than the neutrino heating.It is unexpected and contrary to the situation observed in multi-dimensional simulations.After 50 ms,  gain becomes smaller but still it is larger than  turb gain until 100 ms.Therefore the turbulent energy generation is overestimated in the compressive case until 100 ms.The generated turbulent enegy and turbulent pressure strongly pushes the shock and the shock trajectory in the compressive case is overestimated compared to the 3D model (see Figure 1).Our governing equations for the compressive case is similar to a model in Müller (2019), the model is called non-conservative, no turbulent viscosity.Surprisingly, their calculation did not exhibit an overestimation of the shock position, and we are uncertain about the underlying reasons for this.It could potentially be due to the implementation of the term, see Eq. ( 22) for our particular case.
To understand that turbulent energy enhancement is different in the compressive and non-compressive cases, we conducted a detailed investigation into the radial profile of turbulence in the early phase.Figure 3 shows turbulent velocity distribution at 30 ms postbouce.The horizontal axis is the radius, and the vertical axis is the turbulent velocity.The solid red line corresponds to the non-compressive case, and the dashed red line is the compressive case.The black line is the result for the same time in 3D, which is derived by Reynolds decomposition, which is a way to decompose physical quantities into average and perturbations,  = ⟨⟩ +  ′ , where ⟨⟩ is the solid angular average and is defined as ⟨⟩ = 1 4  ∬ (, , ) sin dd.The velocity of the 3D turbulence shown in the figure is defined as where    = (  − ⟨  ⟩)(  − ⟨  ⟩) is Reynolds tensor.The factor 2 in RHS comes from Eqs. ( 6) and (10), i.e.,  turb = √︁ tr()/2.The turbulent velocity of the compressive case (dashed red) is clearly higher than that of 3D because the turbulent energy is amplified near the shock due to compression.The seed turbulence is ∼ 10 5−6 cm/s outside the shock and undergoes enhancement due to the shock-compression in the early phase, as described in Figure 2. Comparatively, the turbulent pressure is excessively amplified in the compressive models, thereby driving the shock forward.The resultant radial profile diverges significantly from the 3D models.To mitigate this issue, we can omit the compression term.In the noncompressive case the source terms are reasonable (see Figure 2), and the shock evolution is more consistent with 3D (see Figure 1).As an inherent trade-off of this artificial hudling, the turbulent velocity near the shock is underestimated.In Figure 3, the velocity at 80-100 km in the non-compressive case (solid red) is smaller than the 3D case.On the evolution of the shock, this underestimate of the turbulent velocity is less harmful than the overestimate in the compressive case.In the future, we should investigate the method that can mimic 3D profile.The inconsistency between 1D+ model and 3D model starts from the seed turbulent velocity at the pre-shocked region.As discussed in Kazeroni et al. (2018), the growth rate of the convective velocity also depends on the radial velocity and wavenumber of the perturbation.The rate has the maximum in some wavenumber with zero velocity, that is Brunt-Väisälä frequency used in 1D+ model, i.e., the growth rate is overestimated in 1D+ model.If we can keep the seed turbulent velocity in the pre-shocked region as same as that of 3D model, the compressed turbulent energy could become the level of the 3D simulation, that does not significantly alter the shock evolution.Then the compressive model of 1D+ would provide more consistent result with 3D.Since there is no obvious way to keep the seed turbulent velocity small in the pre-shocked region, we prefer using the non-compressive case in order to mimic shock evolution of 3D.We keep this issue in mind when developing the code in the future.
Based on the above argument, we consider the non-compressive case, which is unaffected by the compression, as the fiducial model that can mimic 3D shock evolution.In the following, we present the results of the non-compressive case by systematically changing the turbulence parameters.In Couch et al. (2020),  HT is described as shear production and they reported that the term is negligible.The treatment would be practically similar to ours.Obviously, the current treatments are not perfect.Other sources of the turbulence and the amplification near the shock should be investigated in the future (see Abdikamalov et al. 2016, and references therein).

Dependence of Mixing length parameter
Here, we systematically vary the mixing length parameter in the noncompressive setup and present their dynamics.First, we show the dependence of the shock evolution on the mixing length parameter.Next, we show the velocity distribution of the turbulence and illustrate the effect of the mixing length parameter.
Since the source term of the turbulence is governed by the mixing length and proportional to  Λ in Eq. ( 14), this parameter should play an essential role in the dynamics. Λ significantly affects the shock propagation.In the models with larger  Λ , the stalled shock more easily revive.Here, we explain those issues one by one.
Figure 4 shows the shock evolution of the models of  Λ = 0.0, 1.1-1.3.The horizontal axis is the time from the core bounce, and the vertical axis is the shock radius.The positions of the averaged shock radius of 3D are plotted by black line.The model of  Λ = 0.0 (gray line) is a model that does not employ multi-dimensional effects.In the figure, the later time evolution is depicted, which is not shown in Figure 1.
The model of  Λ = 1.2 can mimic the 3D model well.The shock revival becomes faster if  Λ becomes larger.In the model of  Λ < 1.2, the shock stagnates and fails to explode.On the other hand, for models  Λ ≥ 1.2, the shock revives and propagates out of the iron core.The  Λ monotonically enhances the explodability of the massive stars.This tendency is roughly consistent with the previous study (see Figure 2 of Couch et al. 2020).
Next, we compare the distribution of turbulent velocity of 1D+ with the 3D model.Figure 5 shows the radial profile of turbulent velocity for 110 ms after core bounce for three models  Λ = 1.1-1.3.The horizontal axis is radius.The color of the plot is green for 1.1, orange for 1.2, and red for 1.3 as same as Figure 4.The black line is the turbulent velocity distribution in 3D.
The turbulent velocity of 1D+ peaks near gain radius, 100 km (solid color curves).The profile of the 3D model near the gain radius is similar to the 1D+ models while the width of the turbulence is wider than 1D+.As expected, the turbulent velocity increases monotonically as  Λ increases.This trend is similar to Couch et al. (2020).See their Figure 1.The corresponding peak turbulent pressure,  turb in our model is ∼ 2 × 10 27 erg/cm 3 at 100-120 km.That is larger than that in Figure 2 of Mabanta et al. (2019), ∼ 0.8 × 10 27 erg/cm 3 .Their assumed neutrino luminosity is 2.1 × 10 52 erg/s and is smaller than that of our model, 5 × 10 52 erg/s.Murphy et al. (2013) argues the strength of turbulence is positively correlated with neutrino heating.The difference from their model is reasonable.Note that we use normal Reynolds decomposition while previous study may use different definitions.Our treatment in 3D may exaggerate the peak near the shock (see footnote 2 of Murphy & Meakin 2011).

Diffusion parameters
The effective turbulence model has other parameters than  Λ .The effect of the parameters has not been adequately investigated in the previous works (see Section 1).We systematically change the diffusion constant   ( = , ,   ).The definition of the parameters is given in Eq. ( 9).We fix  Λ = 1.2 in the non-compressive setup as in Section 3.1.
Before going to the details, let us summarize the results in this section.Interestingly,   and   also affect the hydrodynamic evolution.In particular,   has a significant impact on the explodability.Figure 6 compiles the shock evolution in different values of diffusion parameters.The horizontal axis is the time after core bounce and the vertical axis is radius.The top panel shows the shock propagation for the fiducial model and the parameters   (solid In the top panel, the green line is the fiducial model (   ,   ,   ) = (1/6, 1/6, 1/6).We change   (solid line) and   (dashed line) and the red line corresponds   = 1/12, and the blue line is   = 1/3.Green dashed line and solid line are overlapped.
The bottom panel is the same to the top panel but   is changed.All lines are overlapped.
line) and   (dashed line).The green line is the fiducial model (  ,   ,    ) = (1/6, 1/6, 1/6), the red line is 1/12, and the blue line is 1/3.As shown in the top panel of Figure 6, the shock evolution strongly depends on   and   .Surprisingly, the model of   = 1/3 (blue line) fails to explode.The dynamics might be more sensitive to   than   .The shock revival time of   = 1/12 is earlier than   = 1/12.The time is delayed in   = 1/3 and   = 1/3 cannot show the shock revival.On the other hand,    does not change the shock dynamics (see the bottom panel).Note that we switch off the turbulent effect in PNS by hand.In this setup, the diffusion of    does not play any role.

Diffusion of internal energy
Here, we try to interpret why the smaller diffusion parameter,   , make the shock revival easier.This is unexpected and has not been discussed in the literature.Yamasaki & Yamada (2006) showed that the strong energy transport due to the turbulence reduces the critical From top to bottom panel, we show the radial profiles of turbulent velocity, entropy and radial velocity at 100 ms core bounce.We change the parameter   = 1/3, 1/6, and 1/12 and the colors of the curves are blue for 1/3, green for 1/6, and red for 1/12 as same as Figure 6.neutrino luminosity that is necessary for the explosion, i.e., a larger diffusion parameter is favorable for the explosion.Why are our results contrary to this naive expectation?
The earlier shock revival is correlated with longer advection time,  adv (e.g., Buras et al. 2006).In the following paragraphs, we elucidate the underlying reasons why the advection time is longer when the diffusion constant becomes smaller.Figure 7 shows the radial profile of turbulent velocity, entropy and velocity for 100 ms after core bounce for three models   = 1/3, 1/6, and 1/12.The hori- zontal axis is radius.The color of the curve is blue for 1/3, green for 1/6, and red for 1/12 as same as Figure 6.
Examining the top panel and the middle panel, it becomes evident that reducing the diffusion term leads to an increase in turbulent velocity within the range of 60 km to 140 km.This is because turbulent energy generation depends on the Brunt-Väisälä frequency, which is associated with the entropy gradient.Since the diffusion weakens the entropy gradient, a larger turbulent energy generation is expected in the lower diffusion constant.
Consequently, the enhanced turbulent velocity, or stronger turbulent pressure, contributes to reducing the post-shock radial velocity.As observed in the bottom panel, specifically around the region of 120 km to 150 km, the absolute value of the radial velocity is smaller for the red curve (corresponding to   = 1/12) compared to the green and blue curves (representing   = 1/6 and 1/3, respectively).Although the difference is subtle in this particular snapshot, it becomes more pronounced over time.
As a consequence of that, the gain mass becomes larger in the less diffusive models.Figure 8 show the evolution of gain mass [ ⊙ ] (top) and advection timescale (bottom),  adv =  gain /  acc .The horizontal axis is time after core bounce.The color of the curve is blue for 1/3, green for 1/6, and red for 1/12 as same as Figure 6.
In the top panel, the difference in the gain mass is apparent from 100 ms.The detailed radial profile at that time is shown in Figure 7.The model with   = 1/12, which has a slow radial velocity profile, has a larger gain mass.On the other hand, the   = 1/3 model has a smaller gain mass.The larger gain mass corresponds to longer advection timescale as shown in the bottom panel.
To summarize, the lower diffusion constant of internal energy keeps the entropy gradients, leading to efficient turbulence generation.The high turbulent pressure decelerates the radial velocity in the post shock region and consequently the gain mass becomes larger and advection timescale longer.That results in the earlier shock revival.

Diffusion of turbulent energy
We examine the impact of turbulent energy diffusion,   , in a similar vein to the previous section, which discusses the internal energy diffusion.Figure 9 shows the radial profile of turbulent velocity (top), entropy (middle), and velocity (bottom) for 100 ms after core bounce for three models   = 1/3, 1/6, and 1/12 (see Figure 7 for the case of   ).The horizontal axis is radius.The color of the curve is blue for 1/3, green for 1/6, and red for 1/12 as same as Figure 6.
From the middle panel, it is evident that reducing the diffusion constant,   , does not alter the entropy profile.Moving to the top panel, the amplitude of turbulent velocity is higher if the diffusion constant is smaller and the width becomes narrower in the range of 80 km to 150 km.The diffusion blunts the turbulent energy distribution.While the diffusion of internal energy,   , influences the turbulence generation by changing the entropy distribution as explained by Figure 7, the diffusion of turbulent energy does not affect the generation of turbulence, but just alters the turbulent velocity distribution.
The profile of the radial velocity is affected by the turbulent pressure.In the bottom panel, the less diffusive model (red curve) shows the more decelerated profile.The model has the steepest gradients in the middle panel and leads to the highest turbulent pressure gradient force which efficiently slow down the radial velocity.
Consequently, the gain mass becomes larger in the less diffusive model as same as effects of diffusion of internal energy,   .Figure 10 shows Evolution of gain mass [ ⊙ ] (top) and advection timescale (bottom).The horizontal axis is time after core bounce.The color of the curve is blue for 1/3, green for 1/6, and red for 1/12 as same as Figure 6.In the bottom panel of Figure 9, the model with   = 1/12, which has a slow radial velocity profile, has a larger gain mass and longer advection timescale.On the other hand, the   = 1/3 model has a smaller gain mass and shorter advection timescale.Compared to the bottom panel of Figure 7, the variation in gain mass is small.
The diffusion of turbulent energy profoundly impacts gain mass, advection timescale, and the evolution of the shock.The stronger diffusion makes the radial profile of the turbulent velocity flatter leading to smaller turbulent pressure gradient force.The radial velocity of the post shock region is less affected and maintained at higher values.As a result, the gain mass becomes smaller and the advection timescale shorter.

SUMMARY AND DISCUSSION
We have performed 1D simulations with phenomenological turbulent effects (so-called 1D+ simulation).The impact of several turbulent terms are investigated with a 12 ⊙ progenitor model.First the effect of the compression or shear production is focused and then the effects of the mixing length parameter  Λ and the diffusion parameters   ,   ,    are investigated.
We have compared the evolution of the shock among the models with the compression term,  HT and without it.The model without it agree with the 3D simulations.The enhancement of turbulent energy The turbulent velocity, entropy and radial velocity at 100 ms core bounce are shown in the top, middle and bottom panels, respectively.The colors of the curves are blue for   = 1/3, green for 1/6, and red for 1/12 as same as Figure 6.due to the compression, which is defined by Eq. ( 25), is too large at early phase of shock propagation in our 1D+ simulations.
We have adjusted the mixing length parameter so that can mimic the evolution of a 3D shock in the non-compression case.Our model of  Λ = 1.2 mimics well the shock evolution.The turbulent velocity profile also shows good agreement with the 3D turbulent velocity profile, especially near the gain radius where neutrino heating is strong.
We analyze the effect of diffusion parameters   ,   ,    .The standard model is the non-compressive case with  Λ = 1.2.Among them, the diffusion of internal energy,   , most significantly affect the evolution of shock.The diffusion coefficients of internal energy,  , and turbulent energy,   , also affect the explodability.The smaller diffusion makes the shock revival faster.Our comparison between the two reveals that the diffusion coefficients of internal energy has a greater impact.(see Figure 6).  has relatively smaller impact and    has no impact on the explosion in our setup.We analyze the reason why the lower diffusion of   makes the explosion stronger.It is contrary to the expectations of Yamasaki & Yamada (2006), which demonstrate how the energy diffusion helps the shock revival.The lower diffusion of   , the greater the turbulent energy is generated.The turbulent pressure prevents the accretion and the advection time becomes longer.Since   does not affect turbulent generation, the impact becomes smaller.
Afterward, we address caveats in this study.Since 1D+ approach is a relatively new idea, the method has not been maturely established yet.The method based on Reynolds decomposition should include the turbulent mass flux (Müller 2019), which is ignored in this study (it is same as Mabanta et al. 2019;Couch et al. 2020;Boccioli et al. 2021).On the other hand, the method based on Favre decomposition causes a different problem (see Appendix A).More sophisticated methods should be developed.In any case, a detailed comparison would be necessary to understand the difference.
While this study focus on only the convection in the gain region, the convection in proto-neutron star (PNS) would enhance the neutrino luminosity and hence affect the dynamics of the shock (Keil et al. 1996;Buras et al. 2006;Nagakura et al. 2020).Phenomenological treatments on this should be investigated (e.g., Hüdepohl 2014).
Though    does not change the dynamics in our study since we switch off the turbulent effect in PNS by hand, this parameter would be important in the study of PNS convection.Though we have investigated the effect of the shock using two extreme cases, compressive and non-compressive.The treatment is not perfect.If the convection at the early phase is significant, that may change the position of the stalled shock (Nagakura et al. 2018;Harada et al. 2020).The strength of the convection may depend on the detail of the initial perturbation seed (Kazeroni et al. 2018;Kazeroni & Abdikamalov 2020).We need more careful argument on this issue (Takahashi & Yamada 2014;Takahashi et al. 2016;Huete et al. 2018Huete et al. , 2020)).
One essential question is "How accurately does the phenomenological model reproduce the dynamics of the multi-dimensional simulation?".First we need to evaluate the parameters in multi-dimensional simulations.One difficulty with this issue is the numerical resolution of the multi-dimensional simulations.As Nagakura et al. (2019b) reported, the numerical resolution significantly changes the turbulent properties.We need high-resolution simulation to evaluate the parameters (e.g., Radice et al. 2016;Melson et al. 2020).
Perhaps we have to take a significantly different approach to treat convective turbulence.As discussed in Murphy et al. (2013), the turbulent convection in the gain region is driven by the neutrino heating, not directly by negative entropy gradients as in Eq. ( 14), which is valid only in the linear phase.To precisely model the production term, we may have to employ global modeling (e.g., Mabanta et al. 2019) instead of the local modeling, such as MLT (see also Yokoi et al. 2022, for taking non-local effect in the local modeling).In the global modeling, the turbulent effect is not given by the local variables and needs a global structure of the flow, which is characterized by spatial integration of the variables, analytical closure, and empirical relations (Murphy & Meakin 2011;Mabanta et al. 2019).Though we have taken local modeling in this study just due to its simpleness, we need further effort to go to more detailed modeling in the future.
The next important point is to examine whether the flux due to turbulence is consistant with 3D.As long as the gradient-diffusion approximation used in this study is correct, our results should be justified, but more detailed studies are needed.In addition to flux, turbulent pressure is another physical quantity that is important when evaluating turbulence.
1D+ simulations are robust tools for parametric studies of the corecollapse supernovae, though there are still many problems.A variety of applications are possible, for example; progenitor dependence (Couch et al. 2020); neutrino, gravitational wave, optical light curves (Warren et al. 2020;Barker et al. 2022); detailed neutrino emission (Suwa et al. 2019;Mori et al. 2021); EOS dependence (Yasin et al. 2020;Nakazato et al. 2021);BH formation (da Silva Schneider et al. 2020); quark star phase transition (Fischer et al. 2020); emission of axion-like particle (Fischer et al. 2016) and so on.We hope that this study will take us one step closer to understanding a vast and untapped area of the research field.

Figure 1 .
Figure 1.Comparison of the shock evolution in the compressive and the noncompressive cases.The vertical axis is the shock radius, and the horizontal axis represents the time after the bounce.The red solid curve corresponds to the non-compressive case, and the red dashed curve does the compressive case.The black line denotes the time evolution of the average shock radius in the 3D simulation using the same progenitor.

Figure 2 .
Figure 2. Time evolution of source terms of turbulent energy with the noncompressive and compressive cases.The source terms in the turbulent energy equation is integrate over the gain region, see Eqs. (23)-(26).The terms are shown as a function of the postbounce time in ms.The top panel and bottom panels represent the non-compressive case and the compressive case, respectively.Orange, blue, and green indicate turbulent energy source rate, turbulent energy dissipation rate, and turbulent energy increase rate due to turbulence compression in the gain region, respectively.Red is the neutrino heating rate.

Figure 3 .
Figure 3. Radial profile of turbulent velocity at 30 ms postbouce.The horizontal axis is the radius, and the vertical axis is the turbulent velocity.The solid red line corresponds to the non-compressive case, and the dashed red line is the compressive case.The black line is the result of the same time in the 3D simulation.

Figure 4 .
Figure 4. Shock evolution of  Λ = 0.0, 1.1-1.3.The color of the curve is green for 1.1, orange for 1.2, red for 1.3, and gray for 0.0.The position of the average shock radius of 3D is plotted as black curve.

Figure 5 .
Figure 5. Radial profile of turbulent velocity for 120 ms after core bounce for three models  Λ = 1.1-1.3.The color of the plot is as same as Figure 4.The black line is the turbulent velocity distribution in 3D.

Figure 6 .
Figure 6.Shock evolution for different values of diffusion parameters.The horizontal axis is the time after core bounce and the vertical axis is radius.In the top panel, the green line is the fiducial model (   ,   ,   ) = (1/6, 1/6, 1/6).We change   (solid line) and   (dashed line) and the red line corresponds   = 1/12, and the blue line is   = 1/3.Green dashed line and solid line are overlapped.

Figure 7 .
Figure 7. Radial profiles of turbulent velocity, entropy and radial velocity in different   .From top to bottom panel, we show the radial profiles of turbulent velocity, entropy and radial velocity at 100 ms core bounce.We change the parameter   = 1/3, 1/6, and 1/12 and the colors of the curves are blue for 1/3, green for 1/6, and red for 1/12 as same as Figure6.

Figure 8 .
Figure 8.Time evolution of gain mas and advection timescale for   = 1/3, 1/6, and 1/12.The top and bottom panel shows the time evolution of gain mass [  ⊙ ] and  adv , respectively.Horizontal axis is time after core bounce in the top and bottom panel.The colors of the curves are blue for 1/3, green for 1/6, and red for 1/12 as same as Figure6.

Figure 9 .
Figure 9. Radial profiles of turbulent velocity, entropy and radial velocity in different   .Same as Figure 7 but the dependence on   is shown.The turbulent velocity, entropy and radial velocity at 100 ms core bounce are shown in the top, middle and bottom panels, respectively.The colors of the curves are blue for   = 1/3, green for 1/6, and red for 1/12 as same as Figure6.

Figure 10 .
Figure10.Time evolution of gain mas and advection timescale for   = 1/3, 1/6, and 1/12.Same as Figure8but the dependence on   is shown.he top and bottom panel shows the time evolution of gain mass [  ⊙ ] and  adv , respectively.The colors of the curves are blue for   = 1/3, green for 1/6, and red for 1/12 as same as Figure6.

Figure A1 .
Figure A1.Schematic diagram for global energy balance in turbulent convection.The diagrams of the energy flow for Favre and Reynolds decomposition are depicted in Left and Right panels, respectively.