Implicit large eddy simulations of anisotropic weakly compressible turbulence with application to core-collapse supernovae

(Abridged) In the implicit large eddy simulation (ILES) paradigm, the dissipative nature of high-resolution shock-capturing schemes is exploited to provide an implicit model of turbulence. Recent 3D simulations suggest that turbulence might play a crucial role in core-collapse supernova explosions, however the fidelity with which turbulence is simulated in these studies is unclear. Especially considering that the accuracy of ILES for the regime of interest in CCSN, weakly compressible and strongly anisotropic, has not been systematically assessed before. In this paper we assess the accuracy of ILES using numerical methods most commonly employed in computational astrophysics by means of a number of local simulations of driven, weakly compressible, anisotropic turbulence. We report a detailed analysis of the way in which the turbulent cascade is influenced by the numerics. Our results suggest that anisotropy and compressibility in CCSN turbulence have little effect on the turbulent kinetic energy spectrum and a Kolmogorov $k^{-5/3}$ scaling is obtained in the inertial range. We find that, on the one hand, the kinetic energy dissipation rate at large scales is correctly captured even at relatively low resolutions, suggesting that very high effective Reynolds number can be achieved at the largest scales of the simulation. On the other hand, the dynamics at intermediate scales appears to be completely dominated by the so-called bottleneck effect, \ie the pile up of kinetic energy close to the dissipation range due to the partial suppression of the energy cascade by numerical viscosity. An inertial range is not recovered until the point where relatively high resolution $\sim 512^3$, which would be difficult to realize in global simulations, is reached. We discuss the consequences for CCSN simulations.


I. INTRODUCTION
Despite decades of studies and compelling evidence that a significant fraction 1 of stars with initial masses in excess of ∼8 solar masses explode as core-collapse supernovae (CCSN) at the end of their evolution, the exact details of the explosion mechanism are still uncertain [2][3][4][5] . Current state-of-the art 3D simulations either fail to explode or have explosion energies that fall short of the observed energies by factors of a few for most of the progenitor mass range [4][5][6] .
The dynamics at the center of a star undergoing core collapse is shaped by a delicate balance between competing effects where all of the known forces: gravity, electromagnetism, weak and strong interactions, are important. The task of modeling these systems is made particularly challenging by the fact that the generation of the asymptotic explosion energies, although enormous (∼ 10 44 J), requires a rather subtle, percent-level imbalance between nonlinear processes over many dynamical times.
The flow of plasma in the core of a star going supernova is known to be unstable to convection [7][8][9][10] and/or to another large scale instability known as standing accretion shock instability 11,12 . In any case, given the very large Reynolds numbers, as large as ∼ 10 17 in the region of interest 13 (the so-called gain region, where neutrino heating dominates over neutrino cooling), it is expected that the resulting flow will be fully turbulent. It has been suggested 14,15 recently that turbulence and, in particular, turbulent pressure could tip the balance of the forces in favor of explosion. In this respect, anisotropy is of key importance, because it results in an effective radial pressure support with adiabatic index γ turb = 2, much larger than that of thermal (radiation) pressure (γ th ≃ 4/3). This means that turbulent kinetic energy is a much more valuable source of radial pressure support than thermal energy (see Appendix A).
All of the current numerical simulations employ the implicit large eddy simulation (ILES) paradigm 16,17 (also known as monotone integrated LES (MILES)) of exploiting the dissipative nature of high resolution shock capturing (HRSC) methods as an implicit turbulence model. However, the combination of the use of rather dissipative schemes a) Electronic mail: dradice@caltech.edu and the relatively low spatial resolution that can be achieved in global simulations is such that the fidelity with which turbulence is captured is questionable 13 .
To be useful in the context of CCSN simulations, an ILES should, at the very least, account for the right rate of decay of the kinetic energy at the largest scales while avoiding unphysical pile up of energy at smaller scales. Unfortunately, all of the current simulations seem to be strongly dominated by the so-called bottleneck effect 13 , which corresponds to an inefficient energy transfer across intermediate scales due to the viscous suppression of non-linear interaction with smaller scales [18][19][20][21][22] . Current global simulations achieve resolutions, in the turbulent region, comparable to those of 30 3 − 70 3 lattices in periodic domains 13,15,23 . At these resolutions, almost all of the dynamical range of the simulations can be expected to be directly affected by numerical viscosity 24 . The fidelity with which turbulence is captured in these simulations will then depend on the degree with which the numerical truncation error approximates an LES closure.
In this respect, it has been shown by Garnier et al. 25 and Johnsen et al. 26 that many HRSC methods can be too dissipative to yield a faithful description of turbulence at low resolutions. These studies, however, considered a different regime, decaying isotropic turbulence, while turbulence in a core-collapse supernova, as well as in many other astrophysical settings, is often strongly anisotropic 14,15,27 as rotational invariance is broken by gravity. Garnier et al. 25 and Johnsen et al. 26 also considered different numerical schemes with respect to those used in supernova simulations. Both of these aspects can, in principle, be important. First of all, strong anisotropies could potentially influence the turbulence dynamics at the level of the energy cascade and of the dissipation 28 . Secondly, some of the schemes used in computational astrophysics, such as the piecewise parabolic method (PPM) 29 as well as some of the MUSCL 30 schemes, have been shown, differently from some of the methods considered by Garnier et al. 25 and Johnsen et al. 26 , to be well suited for ILES 31,32 .
The aim of this work is to fill the gap between existing theoretical studies and the particular applications of our interest. To this end we use a publicly available code, FLASH, which is widely used in the computational astrophysics community, and perform a series of simulations of turbulence in a regime relevant for core-collapse supernovae: driven at large scale, with large anisotropies and mildly compressible. We use five different numerical setups and, for each, several resolutions in the range from 64 3 to 512 3 in a periodic domain. We study in detail the way in which the energy cascade across different scales is represented by our ILES and we discuss the use of local or lower dimensional diagnostics that can be used to assess the quality of a global simulation in a complex geometry where 3D spectra are not readily available.
The rest of this paper is organized as follows. First, in Section II, we discuss the exact setup of our simulations and the diagnostic quantities used in our analysis. Then, in Section III, we discuss the basic characteristics of the flow realized in our simulations. In Section IV, we present a detailed analysis of the way in which the energy cascade is captured by the different schemes at different scales. In particular, we quantify the accuracy with which different methods capture the decay rate of energy from the largest scales and the way in which energy is distributed across scales. We discuss the role of anisotropies in the context of the 4/5−law, a fundamental exact relation for isotropic and incompressible turbulence relating the statistics of velocity fluctuations with the energy dissipation rate (see Section II C), in Section V. We explore the use of the 2D, transverse, energy spectrum as a diagnostic for 3D simulations in Section VI. Finally, we present a brief summary of our main findings, as well as a discussion of their implications for CCSN simulations in Section VII. Appendix A, contains some supplemental background material on the role of turbulence in the explosion mechanism of CCSN.

A. Numerical methods
We consider a compressible fluid with a prescribed acceleration, a, in a unit-box with periodic boundary conditions. The code that we employ for these simulations, FLASH, solves the gas-dynamics equations in conservation form. In particular we evolve the continuity equation and the momentum equation These equations are closed with a simple isentropic equation of state, that can be considered as a rough description of a gas dominated by radiation pressure. Since the equation of state ensures an adiabatic evolution we do not need to solve the energy equation as equations (1), (2) and (3) suffice to fully describe the flow. Equations (1) & (2) are solved using the directionally-unsplit hydrodynamics solver of the open-source FLASH simulation framework [33][34][35] . FLASH implements the corner transport upwind method 36 for fully directionally-unsplit evolution of the Euler equations 37,38 . FLASH includes several options for the order of spatial reconstruction 35 , including 2 nd -order TVD 30 , 3 rd -order PPM 29 , and 5 th -order WENOZ 39 . Fluxes are computed at 2 nd order accuracy using one of a number of approximate Riemann solvers included in FLASH, such as HLLE 40 , HLLC 41 , and the Roe solver 42 . Second-order accuracy in time is achieved via a characteristic tracing evolution of the Riemann solver input states to the time step midpoint 29 .
Turbulence is driven using the stirring module of FLASH. This module uses the Ornstein-Uhlenbeck process 43 to generate stirring modes in Fourier space. This yields an acceleration field which smoothly decorrelates 44 over a timescale T s . The FLASH implementation permits the use of any arbitrary combination of solenoidal and compressive modes 45 . For our runs, we set T s = 0.5, we use only solenoidal forcing and we restrict the accelerating field to be nonzero only in the first four Fourier modes. This forcing is designed to mimic the influence of some larger scale weakly compressible flow. In the CCSN context, simulations show that the turbulence is highly anisotropic, being roughly twice as strong in the radial direction as either tangential direction 14,15,46,47 since it is driven by buoyancy due to a negative radial entropy gradient. In order to emulate this behavior, the accelerating field in the x−direction (which is going to play the role of the radial direction) is scaled by a constant factor (before the solenoidal projection of the acceleration field) such that R xx ≃ 2R yy ≃ 2R zz , where is the Reynolds stress tensor (to simplify the notation we considered a frame in which ρv = 0) and · denotes an ensemble average. Finally, the overall strength of the stirring is tuned to achieve a RMS Mach number of ≃ 0.3, which is typically observed in realistic CCSN simulations 48,49 .

B. Energy transfer equations
In order to study the cascade of the specific kinetic energy (which we will refer to simply as "kinetic energy" or "energy" in the following), |v| 2 /2, we will consider an energy budget equation across different scales, analogous to the one commonly employed in the study of incompressible, isotropic turbulence, e.g., 50 . In particular, we consider the momentum equation (2) in non-conservation form, where V = 1/ρ is the specific volume of the gas. We can use equation (5) to derive an evolution equation for the Fourier transform of the velocitŷ Transforming both sides of equation (5) we obtain where * denotes the convolution operator, i.e., If we multiply both sides of equation (7) byv * and take the real part, we obtain an equation for the 3D energy spectrum where Here E is the energy spectrum (the velocity power spectral density (PSD)) and T is the same transfer term as in the classical incompressible equations and ǫ is the energy injection rate. The C term vanishes in the incompressible limit and represents the interaction between kinetic and acoustic modes. In practice, in our models, C is found to be at least one order of magnitude smaller than T at all scales and it is thus negligible. In any case, we retain C in the analysis below.
For each of the spectral quantities, S, being E, T, C or ǫ, we define the integrated spectrum, S(k), as δ(·) being the Dirac delta function. Integrating equation (9), we obtain the following one-dimensional energy balance equation This can also be written in terms of the energy flux across scales, as Notice that we did not assume isotropy in any of the above. Equation (15) is derived in the inviscid limit. In practice, our evolution method introduces dissipation in the form of "numerical viscosity". This can be quantified in terms of the residual This can be used to define a wave number dependent "effective viscosity": We remark that ν does not, in general, correspond to a classical shear or bulk viscosity, but can nevertheless be interpreted as a relative measure of the dissipation acting at different wave numbers (see, e.g., [51][52][53] for alternative approaches).
In practice, since we will be working in the stationary case, after having taken the appropriate time averages, R(k) reduces to Finally, since we are working in a periodic domain, which we take of size L x = L y = L z = 1, all of the spectra are quantized and non-trivial only for k x , k y and k z integers. Furthermore, all of the integrals in wave number space reduce to summations. Integrals over spherical shells are transformed to weighted sums following 44 : where N k is the number of discrete wave-numbers in the shell k − 1/2 < |k| ≤ k + 1/2.

C. Structure functions
The energy spectrum and its sources/fluxes give a comprehensive picture of the energy cascade and can be used to assess the level of convergence of the simulation. Unfortunately, 3D energy spectra and fluxes are not easily accessible in calculations in complex domains and/or with inhomogeneous turbulence. In these cases, local quantities in the physical domain are more easily extracted and analyzed. Hence, one of the goals of this work is to validate the use of indirect measures of convergence of ILES. Among these quantities, the structure functions of the velocity appear to be natural candidates for study.
We define the velocity increments and study the quantities where, · j=0 denotes an ensemble average as well as a mean over all of the angles between v and r (in other words we are looking at the j = 0 component of the SO(3) decomposition of the structure functions 54 ). In the case of homogeneous turbulence S p does not depend on x and is thus a function of only the separation r.
The most important relation involving the structure functions is the so-called 4/5−law, which relates the third order structure function, S 3 (r), with the mean energy dissipation rate, and states that, for incompressible, homogeneous and isotropic turbulence 50 : Equation (25) can be derived from the Navier-Stokes equation for fully-developed, incompressible, homogeneous and isotropic turbulence and it is one of the few exact relations in the theory of turbulence 50 . In the anisotropic or compressible case, however, equation (25) is not strictly valid and could be violated in the data. As we show in Section V, we find equation (25) to be very well satisfied by our data, suggesting that the 3 rd order structure function can be a very useful diagnostic in global simulations.

D. Transverse energy spectrum
Another alternative to the analysis of 3D spectra, which has been adopted by several authors in the core-collapse supernova context 13,23,47,55 , is the use of 2D spectra computed using a spherical harmonics expansion of the velocity field tangential to one or more spherical shells in the simulation. Analogously, we emulate this by looking at quantities in the y − z plane and we define the 2D spectra where v ⊥ is the projection of the velocity perpendicular to the x−direction and we introduced the partial Fourier transform of v ⊥ :ṽ In the limit of infinite Reynolds number / resolution, the 2D spectrum is expected to have the same asymptotic behavior as the 3D spectrum, however it is a-priori unclear if E ⊥ is a good proxy for E at finite resolution. For this reason we find it useful to investigate this here. As was the case for the 3D spectra, also here the spectrum is non-trivial only for integer k y and k z , when periodicity is taken into account. The integral in equation (26) is treated analogously to the integral in the equation (14) for the 3D case, while the average in the x−direction in equation (27) is converted to an average over the x−extent of the simulation box.

III. BASIC FLOW PROPERTIES
We employ the finite-volume HRSC (Godunov) approach in which physical states are reconstructed at inter-cell boundaries and local Riemann problems are solved to compute the physical inter cell fluxes. In particular, we perform five groups of simulations using different numerical methods. Each group is labeled using the name of the reconstruction algorithm and of the Riemann solver. For instance TVD HLLE, denotes a group of simulations done using TVD reconstruction and HLLE Riemann solver. Single simulations are labeled using their resolution so that, for instance, TVD HLLE N128, denotes the TVD HLLE run done using a 128 3 grid. For all of the runs the timestep is chosen to have a CFL, i.e., λ∆t/∆x, of 0.4, λ being the maximum characteristics speed, with the exception of the PPM HLLC CFL0.8 runs where we set the CFL to 0.8. For the TVD runs we use the monotonized central (MC) slope limiter 30 . The runs with PPM use the original flattening and artificial viscosity prescriptions from Colella and Woodward 29 . The artificial viscosity coefficient is 0.1. We remark that the use of the artificial viscosity for PPM is not really necessary in this regime 56 , however our goal is not to perform a study of the turbulent dynamics, but to assess how each numerical method performs when used under the same condition as in a real CCSN simulation where strong shocks need to be handled in some parts of the domain.
For each group of simulations we run four resolutions: 64 3 , 128 3 , 256 3 and 512 3 . The RMS velocity in all of the runs is v rms ≃ 0.4, giving an eddy turnover time τ = 1/v rms ≃ 2.5. All of the simulations are run until time t = 100 (≃ 40 eddy turnover times). The time evolution of a few relevant diagnostics is shown in Figure 1 for our fiducial group of runs (PPM HLLC) at different resolutions. We can see how the flow is accelerated from rest and quickly reaches a steady, fully turbulent, state. In all cases, steady state is reached after t 3 (∼ 1 turnover time) and the diagnostics are insensitive to the resolution. The results for the other runs (not shown) are very similar to the ones of PPM HLLC as they all achieve very similar RMS Mach numbers and Reynolds stresses. All of the analysis shown in the rest of the paper are performed using 380 3D snapshots (evenly spaced in time) of the data in the interval 5 ≤ t ≤ 100.
A first, qualitative, comparison between the different methods can be done by looking at their visualizations. In particular, in Figure 2, we show a visualization of the magnitude of the vorticity in the x-z plane for four of the five schemes (excluding PPM HLLC CFL0.8) at the highest resolution (512 3 ). The data is taken at the final time (t = 100). As it can be seen from the figure, all of the simulations show the presence of thin, elongated, regions of high vorticity, as typically seen in direct numerical simulations (DNS) of homogeneous turbulent flows 57,58 . However, the width and the intensity of the vorticity at these smaller scales depend crucially on the numerical scheme. Methods with small intrinsic numerical viscosity, such as PPM HLLC and WENOZ HLLC, present smaller structures and more intermittent vorticity fields with respect to more dissipative methods, such as PPM HLLE and TVD HLLE.
Vortex tubes appear to be preferentially oriented along the x−axis, i.e., along the direction of anisotropy in the driving. Since regions of intense vorticity are also regions of intense dissipation for the classical Navier-Stokes equations, this anisotropy in the vorticity might result in an anisotropic dissipation tensor 28 . However, given the numerical nature of the viscosity that we employ and the relatively low coverage of the inertial range achieved in our simulations (see Section IV), we do not think that a quantitative analysis of this phenomenon is possible using our data. Furthermore, we also point out that in a CCSN anisotropy is not only due to the anisotropic driving, but also due to the fact that the rotational invariance of the hydrodynamics equations is broken by gravity. This will probably have an effect on the nature of vorticity structures at very small scales.

IV. THE ENERGY CASCADE
In this section we focus our analysis on the accuracy with which the energy cascade is captured by our ILES runs. First, we focus on the largest scales of the simulation with the goal of quantifying the accuracy in the decay rate of the energy as a function of the resolution for the different methods. Next, we will look at the energy distribution at smaller scales where, in resolved simulations, the inertial range starts. Finally, we will look at the dynamics in the dissipation region and summarize.  (16), obtained with different numerical methods and resolutions. The energy flux is shown normalized to the average dissipation rate given by equation (24). From left to right and from top to bottom we show the results obtained with PPM HLLC, PPM HLLC CFL0.8, PPM HLLE, TVD HLLE and WENOZ HLLC. The bottom right panel show a comparison of all of the methods at 512 3 . All of the schemes show a good level of accuracy in the energy flux from the largest scales, with errors smaller than a few % already at low resolutions. The differences between the schemes become more marked at large wave numbers where the numerical dissipation starts to interfere with the energy cascade.

A. Energy decay rate
In the limit of very large Reynolds number it is assumed, in standard turbulence phenomenology 50 , that there exists a range of wave numbers (the inertial range) where energy injection and dissipation can be neglected in equation (17). In this range we can write (compressible effects are negligible in our simulations): so that stationarity requires Π(k) ≃ const. In particular, since energy is conserved, one finds Π(k) ≃ ǫ . This means that, in the limit of large Reynolds numbers, the energy decay rate depends only on the macroscopic properties of the flow (and in particular not on the nature of the viscosity), a fact that has also been verified numerically 59 . The significance of this property and its importance for the modeling of turbulence cannot be overstated.
In the context of CCSN simulations this means that the large scale kinetic energy, a crucial quantity for the dynamics of the explosion 15 , can be faithfully captured even with simulations achieving modest Reynolds numbers.
For an ILES, a basic requirement, then, is that a sufficiently high resolution should be achieved to correctly represent the energy cascade at the largest scales. What qualifies as a sufficiently high resolution is of course dependent on the details of the closure built into the scheme (and on the accuracy required for the particular application). To quantify this, we can estimate the level of accuracy that can be reached at any given resolution, using our local simulations. In particular, we can study directly the energy flux across scales, defined by equation (16). This is shown in Figure 3 for all of the different runs,.
As discussed before, we expect that Π(k) ≃ ǫ over an extended region in Fourier space should be a direct indication that a simulation has been able to recover an inertial range. Perhaps not surprisingly, in light of previous results 24 , we find that regions where Π ≃ ǫ as wide as a few wave numbers 4 k 10 only appear at the highest resolutions (we will discuss the inertial range in more detail in Section IV B). However, the amount of energy decaying from the largest scales reaches an asymptotic value much quicker than that implying that the total kinetic energy budget at the largest scales is well resolved even at modest resolutions. We can make a more quantitative statement concerning the energy decay rate by looking at the peak of the energy flux as a function of resolution, as shown in Figure 4. We can see that at 128 3 points all of the simulations have a deviation from the asymptotic energy decay rate of less than 10%. The least dissipative methods already have an error close to the 2% level. A comparison between PPM HLLE and PPM HLLC reveals the profound impact that the choice of the Riemann solver has even at relatively large scale (more on the dissipative properties of the different schemes in Section IV C).

B. Energy spectra
Obviously, not all of the dynamics of turbulence can be reduced to the rate at which kinetic energy decays from the injection scale. The internal dynamics of the energy cascade, far from the injection scale and far from the dissipation range, can also play an important role in many applications. To analyze this aspect we consider in Figure 5 the energy spectrum of the velocity defined by equation (10). The spectra are compensated by k 5/3 to highlight regions with Kolmogorov scaling, which might be expected in the inertial range. Since we want to focus on quantities that do not depend (or depend weakly) on the nature of the energy injection at large scale, we show all of the spectra as a function of a dimensionless wave number, 512 k ∆x. The rationale behind this normalization is that, first of all, we assume the Kolmogorov scale η to be proportional to the grid spacing. Secondly, the 512 factor is introduced to have the dimensionless k, 512 k ∆x coincide with the dimensional one for the highest resolution runs. With this choice, 512 k ∆x = 512 corresponds to a wavelength of a single grid point, 512 k ∆x = 256 corresponds to a wavelength of two grid points and so on.
Looking at any of the groups of runs in Figure 5, one can immediately notice that the spectra obtained at different resolutions do not collapse into a single curve in the dissipation region, as would be required by Kolmogorov's first similarity hypothesis 50 (cf., Gotoh et al. 60 ). This lack of convergence in the dissipation region could be due to the non-linear viscosity of HRSC schemes. This, in turn, could result in an anomalous scaling of η with the grid spacing. Such scaling has been reported in the past for ILES, but it is not very well understood 52 . The good agreement between the three different groups of simulations employing the HLLC Riemann solver seems to support this hypothesis and suggests that the nonlinear viscosity introduced by the Riemann solver is an important ingredient in setting this scaling.
Convergence appears to be recovered at larger scales 8 ∆x (512 k ∆x 64), but the spectra appear to be dominated by the bottleneck effect. This manifests itself as a bump in the compensated spectra extending from the dissipation range until the end of the inertial range, for the simulations that show one (e.g., until 512 k ∆x = 10 for the HLLC runs), or until the energy injection scale (512 k ∆x = 4), for the simulations that show no or little inertial range (TVD HLLE). The bottleneck effect is a viscous phenomenon which is also observed in direct numerical simulations. However, in the present context where viscosity is of numerical origin, it is at the very least questionable if a pronounced bottleneck is a desirable feature of the modeling. In astrophysical flows, where the Reynolds numbers are typically very large, this pile up of energy at large scales is unphysical and could affect the quantitative and qualitative outcome of a simulation 13 . A quantification of the bottleneck effect in terms of the energy budget is discussed in Section IV D. At even larger scales, an inertial range (E ∼ k −5/3 and Π ∼ const, see Figure 3) seems to be recovered by the least dissipative schemes (PPM and WENOZ with HLLC) in the region 4 k 10. PPM HLLE and TVD HLLE have a more limited region, a few wave numbers at most, that could be interpreted as being an inertial range. We note that this resolution is not particularly high in comparison with state of the art DNS 59,61 , but it would already correspond to an extremely high resolution in global CCSN simulations that typically have ∼ 30 − −70 zones across the turbulent region 13 .
The overall behavior of the spectra, as obtained by all schemes, is consistent with Kolmogorov's theory of turbulence. The anisotropic contributions to the angle-integrated spectra are too small to be detected in our data.

C. Effective viscosity
At very small scales (∼ several grid points) the dynamics is dominated by the numerical viscosity. This can be estimated from the residual of the energy equation (17) or, equivalently, by the effective numerical viscosity ν(k) (equation 19). The latter is shown in Figure 6 for all of schemes and resolutions.
The first thing to notice is that the effective viscosity provided by all numerical schemes is not constant, but differs by roughly an order of magnitude between low and high k. Having a wave number dependent viscosity is a desirable feature expected in any LES model (explicit or otherwise). Nevertheless, this makes the definition and calculation of the Reynolds number achieved in a simulation ambiguous. Meaningful ways to estimate the Reynolds number from ILES have been proposed 53 and they can be used to ease the comparison between different simulations and assess their quality. However, one has to be very careful while using any quoted Reynolds number from an ILES, to estimate things like the dynamical range achieved by a simulation. Two other features can be observed in most of the effective viscosity profiles. First, many of them exhibit a sudden reversal at high wave numbers. This is due to the fact that the numerical viscosity does not behave like a shear viscosity so that, although the numerical diffusion is strong at those scales, the effective viscosity appears small because of a partial decoupling between vorticity and dissipation. Second, at high resolution and at the largest scales, the effective viscosity is close to zero or even slightly negative. The reason is that the residual of equation (9) oscillates around zero and it is too small to be reliably extracted from our data: a much longer integration time would be needed to accumulate enough statistics for it.
Finally, a comparison between the effective viscosity reveals two interesting effects. First, by comparing PPM HLLC and PPM HLLE, we see that the choice of the Riemann solver affects the viscosity at basically all scales. Second, if we compare PPM HLLC, PPM HLLC CFL0.8 and WENOZ HLLC, we see that doubling the timestep has an effect comparable to changing the reconstruction method at intermediate scales (40 k 100).

D. The energy distribution
So far we have been concerned with the energy decay rate from the largest scales, which we have shown to be well captured by the ILES (Section IV A), and with the energy transfer in the inertial range, which we have seen to be described accurately only at much higher resolutions (Section IV B). In a turbulent flow both of these aspects are important and a good ILES should display a distribution of energy across vortical structures at different scales that is as close as possible to the asymptotic one. Obviously, there is a limit to the accuracy that any ILES can achieve at a fixed resolution. Here, we make this statement more quantitative by considering the amount of kinetic energy that is well resolved by each simulation at a given resolution. We introduce the cumulative energy spectrum, the integral of the energy spectrum: This quantity is plotted in Figure 7, where it is normalized by to obtain the cumulative distribution function of the kinetic energy. As a reference, we also show the cumulative energy spectrum estimated from Kolmogorov's theory: We find that as the resolution increases, all schemes appear to be converging to the predictions of Kolmogorov's theory. The results at finite resolution, however, are not encouraging: at 64 3 only ∼ 50% or less of the kinetic energy is in well resolved structures, while the other ∼ 50% have piled up at rather large scale, with a cumulative excess of ∼ 20% at the grid scale, mostly because of the bottleneck effect. At higher resolutions, the amount of kinetic energy well captured by the ILES increases, but at 512 3 this is still only about 80% of the energy and there is still a cumulative excess of over ∼ 5% at the grid scale (ℓ ∼ ∆x).

V. THE 4/5−LAW
The 4/5−law (equation 25) is not a-priori valid in the regime of turbulence we are considering. However, the 4/5−law has been numerically verified to hold also in some situations outside the domain of validity of its derivation. For instance, for isotropic mildly compressible decaying 62 and driven 63 turbulence. In the anisotropic case, however, anisotropic contributions cannot be excluded 64 , although they are known to be subdominant in some important cases [65][66][67] . In this section we show that equation (25) is consistent with our data over a wide range of scales. We compute the 3 rd order structure functions of the velocity, defined by equation (23), in a rather simple way using a random sample of 20,000 points in each of the 380 3D data dumps of our simulations. At each time, we compute the 3 rd power of the velocity increments for each pair of points and accumulate and average in time the results in bins of size ∆x. The resulting structure functions are shown in Figure 8, compensated by − 5 4 r −1 ǫ −1 , so that the resulting quantity should be equal to one if the 4/5−law is satisfied in our data. As was the case for the energy spectra, we assume η ∼ ∆x and plot the structure functions versus r/∆x.
The degree with which the 4/5−law is satisfied in our data is very good. We see that anisotropic contributions only play a minor role in the angle-integrated formulation of the 4/5−law. This is in agreement with the incompressible DNS of Kaneda et al. 67 and has been known to be true also for Rayleigh-Bénard convection in most regimes 68 . Our results provide an important new example where this appears to hold true. Secondly, for all of our simulations at 512 3 , we find within 5% of 1. This level of accuracy is reached in DNS simulations achieving at least a Taylor micro-scale Reynolds number 67 where ν is the kinematic viscosity, u ′ = 1 these Reynolds numbers in a DNS requires resolutions between 512 3 and 1024 3 using pseudo-spectral methods 69 . This large-scale estimate of the Reynolds number is consistent with previous findings 53 , although it is several orders of magnitude larger than the one that can be estimated using ν max . For instance, for PPM HLLC at 512 3 , ν max ≃ 1.5 × 10 −3 and v rms ≃ 0.4 giving This apparent discrepancy is due to the fact that an ILES is not a DNS and, in particular, there is no unique Reynolds number for an ILES. As a consequence, different quantities that in a DNS depend on the Reynolds number, such as the dissipation rate or the Kolmogorov scale, behave as though the simulation has multiple values of the Reynolds number.

VI. THE TRANSVERSE SPECTRUM
Finally, we want to comment on the use of 2D transverse spectra in 3D simulations, a practice typically employed in the analysis of turbulence in CCSN simulations 13,23,47,55 . Figure 9 shows a comparison of the 2D transverse spectrum E ⊥ (k ⊥ ) from equation (26) and the 3D energy spectrum from equation (10) for the PPM HLLC simulations. The other runs (not shown) have the same qualitative behavior. We can see that the transverse spectrum follows qualitatively the same trend as the 3D spectrum in terms of convergence. They are both roughly compatible with a Kolmogorov scaling, but the bottleneck appears to be more pronounced in the 2D spectrum than in the 3D spectrum. In particular, E ⊥ (k ⊥ ) only shows a very small region that suggests an inertial range, 3 k 5 (as opposed to 5 k 10 in E(k)).
Abdikamalov et al. 13 concluded, also based on the analysis of 2D spectra, that turbulence in CCSN simulations is dominated by the bottleneck effect. Given the resolutions used in CCSN studies, our work supports their conclusion. However, in the light of Figure 9, we recommend that future studies supplement the analysis of 2D spectra with 3 rd order structure functions, that, as we have shown, can give a more accurate description of the energy cascade.

VII. CONCLUSIONS
The details of the explosion mechanism of CCSNe have eluded our comprehension in spite of more than 50 years of studies [2][3][4][5] . Recent numerical advances 14,15,48,49 suggest that turbulence might play a fundamental role in tipping over the balance of the forces and lead to successful explosions (see also Appendix A). At the same time, the level of accuracy of current simulations, which employ the ILES methodology, is unclear 13 . Turbulence in CCSNe is mildly compressible, but strongly anisotropic 14,15 . Simulations use rather dissipative numerical schemes (because they have to deal with strong shock waves and complex microphysics) and relatively low resolution, a combination (anisotropic turbulence and dissipative schemes) that has not been systematically studied before.
With the goal of assessing the reliability of ILES employed in the study of CCSNe, as well as in other areas of physics and astrophysics, we performed a series of local simulations of driven, anisotropic, weakly compressible turbulence. We compared five commonly employed numerical schemes with different reconstruction methods, Riemann solvers, and time step size. Each was run at 4 different resolutions ranging from 64 3 to 512 3 . Our analysis focused on the fidelity with which the turbulent cascade is represented in each model. In particular, we performed an analysis both in Fourier space (with the velocity power-spectra and the energy flux) and in physical space (with the 3 rd order structure functions). Finally, we measured the effective viscosity of each scheme from the residual of the specific kinetic energy equation.
We found that, on the one hand, all of the numerical setups are able to accurately capture the decay rate of kinetic energy from the injection scale, with errors at the few % level already at 128 3 (e.g., ∼ 2.5% for PPM HLLC N128). On the other hand, a large fraction of the energy is at unresolved scales where it piles up due to the bottleneck effect and an inertial range appears only at the highest resolutions (512 3 ). Even at this resolution, which would be difficult to achieve in global simulations, only roughly ∼ 80% (the exact number depends on the scheme, see Section IV D) of the energy is resolved, the remaining ∼ 20% accumulates as excess energy at intermediate scales (the cumulative energy excess at the grid scale alone is as large as ∼ 5% of the total energy).
Current CCSN simulations have resolutions of at most of 30 − 70 points covering the gain region 13 (the energy injection scale). Based on our analysis we expect that at these resolutions even the energy decay rate from the largest scales will not be completely converged, but will show errors of up to tens of percent, depending on the numerical scheme (see Section IV A). At smaller scales, the dynamics is going to be completely dominated by the bottleneck effect. This is in agreement with the findings of Abdikamalov et al. 13 , based on the use of global simulations reaching a maximum resolution of 66 grid points covering radially the extent of the gain region.
Based on our findings, we expect that, if the resolution in global simulations is increased by a factor ∼ 2 from the one of Abdikamalov et al. 13 , the decay rate will be converged to within a few % of the asymptotic value. This implies that the ratio between thermal and kinetic energy, a crucial quantity for the onset of the explosion, will also be converged to within a few %, at least when the energy injection rate changes slowly compared to the eddy turnover time (which is roughly ∼ 20 ms in a CCSN 23,70 ). Unfortunately, while the lead up to explosion occurs over a larger timescale of a few hundred milliseconds, the transition to explosion can happen over much shorter timescales (one turnover time or less) 15 . This means that the dynamics of the cascade over smaller time and length scales in the gain region also needs to be captured correctly since changes in the energy input rate on such short time scales will yield an inaccurate representation of the energy on large scales do to the bottleneck effect. This could require an increase of resolution by a factor ∼ 4 − 8 with respect to current high-resolution simulations. Additional work using semi-global or global simulations will be required to more firmly establish the resolutions requirements at the transition of the explosion.
Concerning the properties of anisotropic turbulence in our simulations, we found anisotropy contributions to the energy spectrum and to the angle-averaged formulation of the 4/5−law to be subdominant: the accuracy with which the 4/5−law is satisfied is limited only by the employed resolution and the energy spectrum appears to be consistent with Kolmogorov k −5/3 scaling. We also found the transverse energy spectrum with respect to the direction of anisotropy, a quantity typically computed in CCSN simulations, to overestimate the bottleneck with respect to the angle-integrated 3D spectrum. For this reason, we recommend future studies of CCSN to supplement (or replace) the analysis of the transverse spectrum with the analysis of the 3 rd order, angle-integrated, structure function (or, where possible, with the 3D spectrum itself).
Our results are, of course, dependent on the choice of the numerical scheme. In particular, we found significant differences in the dissipative properties of schemes employing the HLLC Riemann solver with respect to schemes using the more dissipative HLLE solver. The reconstruction order of the scheme is also important, although, while significant differences are found between TVD and PPM, the differences between PPM and WENOZ are much more minute (despite WENOZ being significantly more computationally expensive than PPM). In the end, none of the schemes we considered seems to be able to yield an accurate representation of the kinetic energy distribution across different scales at an affordable resolution for global CCSN simulations. A possible way forward would be to adopt low-dissipation numerical schemes especially designed for the use in ILES, such as the methods proposed by Hickel et al. 71 , Martín et al. 72 or Thornber et al. 73 . Implementing and testing these schemes will be subject of future work.
An important limitation of the present work is that we considered a very idealized setup. On the one hand, this allows us to benchmark the behavior of ILES in a controlled environment. On the other hand, our simulations cannot fully capture all features of the turbulent convective flow in a CCSN. Unlike the situation in a CCSN, our local simulations did not include a vertical advective velocity field that is due to the accretion of the stellar mantle. However, the advective velocities are nearly constant in the regions of interest and Galilean invariance ensures that our results are unaffected. More limiting is the local nature of our simulations and the inevitable choice of boundary conditions. Moreover, our simulations could not take into account spatial variations in gravity and the large-scale radial convergence of the flow in globally spherical problems like collapsing stars. Addressing these issues will also be subject of future work.
In this appendix we present a brief discussion of the importance of turbulent pressure in the explosion of massive stars. To set the stage, we will briefly summarize what is known of the dynamics of the most common class of CCSNe that are relevant for our later discussion. This is done for the benefit of readers that are not supernova specialists and it is not meant to be a complete review of the status of the field, for which we refer, instead, to the reviews of Janka et al. 3 and Burrows 4 . Next, we will discuss the role of turbulence and, in particular, of turbulent pressure on the explosion mechanism, in light of some recent results 14,15 .

The neutrino mechanism
Towards the end of their evolution, massive stars form massive (∼ 1.5 solar mass) iron cores at their center. Since the iron nucleus has the largest binding energy per nucleon, no energy can be extracted from nuclear fusion beyond iron. The iron core is essentially inert and supported against gravity only by the degeneracy pressure of relativistic electrons. The mass of the iron core increases with time as more iron-group material is added by silicon shell burning. Electron capture on protons, which becomes energetically favorable at high densities, depletes the core of electrons and thus reduces the pressure supporting it against gravity. Eventually, the core becomes dynamically unstable and collapses.
During the collapse, the subsonically collapsing inner core (∼ 0.5 solar masses) contracts until it reaches densities comparable to that in atomic nuclei (∼ 4 − 5 × 10 14 g cm −3 ). At this point, the nuclear equation of state stiffens (due to the strong nuclear force). This halts the collapse of the inner core. It stops, bounces back and a proto-neutron star (PNS) is formed. The outer core, however, is still collapsing supersonically and a strong shock wave is launched at the interface between the inner and outer core.
It was once thought that this shock wave would travel outwards dynamically, depositing its energy in the outer layers of the star, causing the explosion. However, multiple numerical simulations performed over multiple decades have consistently shown that the initial shock fails to explode the star. Instead, it stalls due to energy losses to the dissociation of heavy nuclei into free nucleons and to the emission of neutrinos that stream away from the neutrinosemitransparent regions behind the shock 74 . The shock generally stalls within only a few tens of milliseconds of core bounce and turns into an accretion shock standing at a radius of ∼ 100 − 200 km. The accretion rate through the shock is so high (a fraction of a solar mass per second) that, if nothing revitalizes the shock within ∼ 1 − 2 seconds, the gravitational force would overwhelm the nuclear repulsion force, collapsing the core of the supernova to a black hole, precluding explosion (e.g., 75 ).
During this time, however, the PNS will release a significant fraction of its binding energy in the form of neutrinos (of order 10 46 J). Converting a few percent of that energy into kinetic energy would be enough to unbind the stellar envelope and power the supernova explosion. In the standard neutrino mechanism it is theorized that a small fraction neutrinos emitted from the edge of the protoneutron star is re-absorbed in the region right behind the stalled shock. The deposition of neutrino energy leads to higher thermal pressure so that the shock can eventually overcome the ram pressure of accretion and accelerates in a run-away process 74,76 . Turbulence in the heating region behind the shock increases the time a fluid parcel spends in that region and, importantly, turbulent pressure helps in overcoming the ram pressure of accretion (see next Section and 23 ). It is, however, presently unclear if neutrino heating (even if aided by turbulence in launching the explosion) is able to provide enough energy to power the explosions to the energies inferred from astronomical observations.

Turbulent pressure and the Rankine-Hugoniot conditions
Simulations 8,14,15 have shown that turbulence and, in particular, turbulent pressure behind the shock, could play an important role in aiding the explosion. To see why this is the case, we consider the Rankine-Hugoniot momentum condition for a standing accretion shock in a supernova core, where s is the shock speed and · d and · u denote the downstream and upstream values respectively. For the purpose of our discussion, we can assume the upstream gas to be cold and free-falling: For the shock to expand we must then have In the presence of turbulence, Murphy et al. 14 suggested to modify equation (A3) in a way akin to a Reynolds decomposition and write it as where v is the average velocity and δv =v−v is the turbulent velocity. Although not entirely rigorous, equation (A4) has been shown to be well verified in the numerical simulations if angular averages are used to compute the respective quantities 14,15 . Couch and Ott 15 have shown that the turbulent pressure expressed in this fashion can exceed 50% of the thermal pressure, making a very significant contribution to the momentum balance in (A4). Going beyond the arguments of Murphy et al. 14 , we can reinterpret equation (A4) as being the Rankine-Hugoniot condition for a fluid with a modified equation of state, which has two separate internal degrees of freedom: thermodynamical and turbulent. To this aim, we express δv r d in terms of the specific turbulent energy e turb = 1 2 |δv| 2 (A5) noting that and using the fact that (δv r ) 2 ≃ (δv θ ) 2 + (δv φ ) 2 (A7) in CCSN turbulence, to obtain (δv r ) 2 ≃ 1 2 |δv| 2 = e turb .
Assuming the pressure varies like p = (γ − 1)ρe, and substituting (A8) into (A4), we find ρ d (v r d ) 2 + (γ th − 1)ρ d e d + (γ turb − 1)ρ d e turb > ρ u GM r , where γ th ≃ 4/3 is the thermodynamical adiabatic index, e d is the downstream thermal energy and γ turb = 2 is the equivalent adiabatic index associated with anisotropic CCSN turbulence. Since γ turb > γ th , we see that turbulent energy is more efficient, per unit specific internal energy, at pushing the shock than thermal energy.
We point out that, if equation (A7) is dropped and turbulence is assumed to be isotropic, then γ turb = 5/3, which is still larger than γ th , but not as large as for the anisotropic case. This is a simple consequence of the fact that anisotropic turbulence has an anisotropic pressure, which is stronger in the radial direction.
In both cases, since the total energy is conserved, the relevant quantity is the ratio e turb /e. From standard turbulent phenomenology we expect that this ratio will only depend on macroscopic parameters, such as the net heating rate, the accretion rate and so on, and not on the details of the viscosity. For this reason, we expect this ratio to be correctly captured in ILES achieving a sufficiently high resolution.
As a final remark, we point out that a similar argument has been recently proposed by Mueller and Janka 49 who formulated their equations in terms of the turbulent Mach number, as opposed to the turbulent energy.