Skip to main content

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

Abstract

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. The ILES approach has been applied to different contexts, with varying degrees of success. It is the de-facto standard in many astrophysical simulations and in particular in studies of core-collapse supernovae (CCSN). 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. Anisotropy, in particular, could impact the dissipative properties of the flow and enhance the turbulent pressure in the radial direction, favouring the explosion. 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. Our simulations employ several different methods and span a wide range of resolutions. 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 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, i.e., 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 high resolution 5123, which would be difficult to realize in global simulations, is reached. We discuss the consequences for CCSN simulations.

1 Introduction

Despite decades of studies and compelling evidence that a significant fraction (Clausen et al. 2015) 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 (Woosley and Janka 2005; Janka et al. 2012; Burrows 2013; Foglizzo et al. 2015). 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 (Janka 2012; Burrows 2013; Foglizzo et al. 2015).

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 (\({\sim}10^{44}~\mathrm{J}\)), requires a rather subtle, percent-level imbalance between non-linear processes over many dynamical times.

The flow of plasma in the core of a star going supernova is known to be unstable to convection (Herant 1995; Burrows et al. 1995; Janka and Müller 1996; Foglizzo et al. 2006) and/or to another large scale instability known as standing accretion shock instability (Blondin et al. 2003; Foglizzo et al. 2007). In any case, given the very large Reynolds numbers, as large as 1017 in the region of interest (Abdikamalov et al. 2015) (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 (Murphy et al. 2013; Couch and Ott 2015) 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 \(\gamma_{\mathrm{turb}} = 2\), much larger than that of thermal (radiation) pressure (\(\gamma_{\mathrm{th}} \simeq4/3\)). This means that turbulent kinetic energy is a much more valuable source of radial pressure support than thermal energy (see Appendix).

All of the current numerical simulations employ the implicit large eddy simulation (ILES) paradigm (Garnier et al. 2000; Grinstein et al. 2011) (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 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 (Abdikamalov et al. 2015).

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 (Abdikamalov et al. 2015), which corresponds to an inefficient energy transfer across intermediate scales due to the viscous suppression of non-linear interaction with smaller scales (Yakhot and Zakharov 1993; She and Jackson 1993; Falkovich 1994; Verma and Donzis 2007; Frisch et al. 2008). Current global simulations achieve resolutions, in the turbulent region, comparable to those of 303-703 lattices in periodic domains (Couch and O’Connor 2014; Couch and Ott 2015; Abdikamalov et al. 2015). At these resolutions, almost all of the dynamical range of the simulations can be expected to be directly affected by numerical viscosity (Sytine et al. 2000). 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. (1999) and Johnsen et al. (2010) 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 (Arnett et al. 2009; Murphy et al. 2013; Couch and Ott 2015) as rotational invariance is broken by gravity. Garnier et al. (1999) and Johnsen et al. (2010) 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 (Casciola et al. 2007). Secondly, some of the schemes used in computational astrophysics, such as the piecewise parabolic method (PPM) (Colella and Woodward 1984) as well as some of the MUSCL (Toro 1999) schemes, have been shown, differently from some of the methods considered by Garnier et al. (1999) and Johnsen et al. (2010), to be well suited for ILES (Schmidt et al. 2006; Thornber et al. 2007).

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 (Fryxell et al. 2000; Dubey et al. 2009; Lee et al. 2014), 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 643 to 5123 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 2, we discuss the exact setup of our simulations and the diagnostic quantities used in our analysis. Then, in Section 3, we discuss the basic characteristics of the flow realized in our simulations. In Section 4, 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 2.3), in Section 5. We explore the use of the 2D, transverse, energy spectrum as a diagnostic for 3D simulations in Section 6. Finally, we present a brief summary of our main findings, as well as a discussion of their implications for CCSN simulations in Section 7. Appendix contains some supplemental background material on the role of turbulence in the explosion mechanism of CCSN.

2 Methods

2.1 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

$$ \partial_{t} \rho+ \nabla\cdot( \rho \mathbf{v}) = 0 $$
(1)

and the momentum equation

$$ \partial_{t} (\rho\mathbf{v}) + \nabla\cdot(\rho \mathbf{v}\otimes\mathbf{v}+ p \mathbf{I}) = \rho \mathbf{a}. $$
(2)

These equations are closed with a simple isentropic equation of state,

$$ p = \rho^{4/3}, $$
(3)

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) and (2) are solved using the directionally-unsplit hydrodynamics solver of the open-source FLASH simulation framework. FLASH implements the corner transport upwind method (Colella 1990) for fully directionally-unsplit evolution of the Euler equations (Lee and Deane 2009; Lee 2013). FLASH includes several options for the order of spatial reconstruction (Lee et al. 2014), including 2nd-order TVD (Toro 1999), 3rd-order PPM (Colella and Woodward 1984), and 5th-order WENOZ (Borges et al. 2008). Fluxes are computed at 2nd-order accuracy using one of a number of approximate Riemann solvers included in FLASH, such as HLLE (Einfeldt 1988) and HLLC (Toro et al. 1994). Second-order accuracy in time is achieved via a characteristic tracing evolution of the Riemann solver input states to the time step midpoint (Colella and Woodward 1984). We remark that, in accordance with the ILES, paradigm, we do not include any additional sub-grid scale model, but relied on the implicit turbulent closure built in the numerical schemes we use for the integration of the hydrodynamics equation.

All of our simulations start with the fluid at rest \(\rho= 1\), \(\mathbf{v}=0\). Turbulence is driven using the stirring module of FLASH. This module uses the Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein 1930) to generate stirring modes in Fourier space. This yields an acceleration field which smoothly decorrelates (Eswaran and Pope 1988) over a timescale \(T_{s}\). The FLASH implementation permits the use of any arbitrary combination of solenoidal and compressive modes (Federrath et al. 2010). 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 and, for this reason, it does not include any compressible component. This is a reasonable approximation for low Mach number convection which is well described by the anelastic approximation, e.g., Verhoeven et al. (2015). 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 (Murphy and Meakin 2011; Murphy et al. 2013; Handy et al. 2014; Couch and Ott 2015) 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} \simeq2 R_{yy} \simeq2 R_{zz}\), where

$$ R_{ij} = \langle\rho v_{i} v_{j} \rangle, $$
(4)

is the Reynolds stress tensor (to simplify the notation we considered a frame in which \(\langle\rho\mathbf{v}\rangle= 0\)) and \(\langle\cdot\rangle\) denotes an ensemble average. Finally, the overall strength of the stirring is tuned to achieve a RMS Mach number of 0.35, which is typically observed in realistic CCSN simulations (Couch and Ott 2013; Müller and Janka 2015).

2.2 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), \(\vert \mathbf{v}\vert ^{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., Frisch (1996). In particular, we consider the momentum equation (2) in non-conservation form,

$$ \partial_{t} \mathbf{v}+ (\mathbf{v}\cdot\nabla) \mathbf{v}= - \mathcal {V}\nabla p + \mathbf {a}, $$
(5)

where \(\mathcal {V}= 1/\rho\) is the specific volume of the gas.

We can use equation (5) to derive an evolution equation for the Fourier transform of the velocity

$$ \hat{\mathbf{v}}(\mathbf {k}) = \int_{\mathbb{R}^{3}} e^{-2\pi \mathrm {i}\mathbf {k}\cdot \mathbf {x}} \mathbf{v}(\mathbf {x}) \,\mathrm {d}^{3} \mathbf {x}. $$
(6)

Transforming both sides of equation (5) we obtain

$$ \partial_{t} \hat{\mathbf{v}} + \hat{\mathbf{v}} \ast2 \pi \mathrm {i}\mathbf {k}\otimes\hat{\mathbf{v}} = - \hat{\mathcal {V}} \ast2 \pi \mathrm {i}\mathbf {k}\hat{p} + \hat{ \mathbf {a}}, $$
(7)

where denotes the convolution operator, i.e.,

$$ [f \ast g](\mathbf {k}) = \int_{\mathbb{R}^{3}} f(\mathbf {q}) g(\mathbf {k}-\mathbf {q}) \,\mathrm {d}^{3}\mathbf {q}. $$
(8)

If we multiply both sides of equation (7) by \(\hat{\mathbf{v}}^{\ast}\) and take the real part, we obtain an equation for the 3D energy spectrum

$$ \partial_{t} E(\mathbf {k}) = T(\mathbf {k}) + C(\mathbf {k}) + \epsilon(\mathbf {k}), $$
(9)

where

$$\begin{aligned}& E(\mathbf {k}) = \frac{1}{2} \hat{\mathbf{v}}\cdot\hat{ \mathbf{v}}^{\ast}, \end{aligned}$$
(10)
$$\begin{aligned}& T(\mathbf {k}) = - 2\pi \Re (\hat{\mathbf{v}} \ast \mathrm {i}\mathbf {k}\otimes\hat {\mathbf{v}}) \cdot \hat{\mathbf{v}}^{\ast}, \end{aligned}$$
(11)
$$\begin{aligned}& C(\mathbf {k}) = - 2\pi \Re (\hat{\mathcal {V}}\ast \mathrm {i}\mathbf {k}p) \cdot\hat {\mathbf{v}}^{\ast}, \end{aligned}$$
(12)
$$\begin{aligned}& \epsilon(\mathbf {k}) = \Re \hat{\mathbf {a}}\cdot\hat{\mathbf{v}}^{\ast}. \end{aligned}$$
(13)

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

$$ S(k) = \int_{\mathbb{R}^{3}} S(\mathbf {k}) \delta\bigl(\vert \mathbf {k}\vert - k\bigr) \,\mathrm {d}\mathbf {k}, $$
(14)

\(\delta(\cdot)\) being the Dirac delta function.

Integrating equation (9), we obtain the following one-dimensional energy balance equation

$$ \partial_{t} E(k) = T(k) + C(k) + \epsilon(k). $$
(15)

This can also be written in terms of the energy flux across scales,

$$ \Pi(k) = - \int_{0}^{k} T(\xi) \, \mathrm {d}\xi, $$
(16)

as

$$ \partial_{t} E(k) + \partial_{k} \Pi(k) = C(k) + \epsilon(k). $$
(17)

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

$$ R(k) = \partial_{t} E(k) - T(k) - C(k) - \epsilon(k). $$
(18)

This can be used to define a wave number dependent numerical viscosity:

$$ \nu(k) = - \frac{1}{2} \frac{R(k)}{k^{2} E(k)}. $$
(19)

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., Fureby and Grinstein (1999); Aspden et al. (2009); Zhou et al. (2014) 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

$$ R(k) = - T(k) - C(k) - \epsilon(k). $$
(20)

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 Eswaran and Pope (1988):

$$ E(k) = \frac{4\pi k^{2}}{N_{k}} \sum_{k-1/2 < \vert \mathbf {k}\vert \leq k+1/2} E(\mathbf {k}), $$
(21)

where \(N_{k}\) is the number of discrete wave-numbers in the shell \(k-1/2 < \vert \boldsymbol{k}\vert \leq k+1/2\).

2.3 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

$$ \delta v(\mathbf {x}, \mathbf {r}) = \bigl[\mathbf{v}(\mathbf {x}+ \mathbf {r}) - \mathbf{v}(\mathbf {x}) \bigr]\cdot\frac{\mathbf {r}}{r} $$
(22)

and study the quantities

$$ S_{p}(r) = \bigl\langle \delta v^{p} \bigr\rangle _{j=0}, $$
(23)

where, \(\langle\cdot\rangle_{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 \(\mathrm{SO}(3)\) decomposition of the structure functions (Biferale and Procaccia 2005)). 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,

$$ \langle\epsilon\rangle= \int_{0}^{\infty}\epsilon(k) \,\mathrm {d}k, $$
(24)

and states that, for incompressible, homogeneous and isotropic turbulence (Frisch 1996):

$$ S_{3}(r) = - \frac{4}{5} \langle\epsilon\rangle r. $$
(25)

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 (Frisch 1996). 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 5, we find equation (25) to be very well satisfied by our data, suggesting that the 3rd order structure function can be a very useful diagnostic in global simulations.

2.4 Transverse energy spectrum

Another alternative to the analysis of 3D spectra, which has been adopted by several authors in the core-collapse supernova context (Dolence et al. 2013; Couch and O’Connor 2014; Handy et al. 2014; Abdikamalov et al. 2015), 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

$$ E_{\perp}(k_{\perp}) = \frac{1}{2} \int _{\mathbb{R}^{2}} \tilde{\mathbf {v}}_{\perp}\cdot \tilde{ \mathbf{v}}^{\ast}_{\perp}\delta \Bigl(\sqrt{k_{y}^{2} + k_{z}^{2}} - k_{\perp}\Bigr) \,\mathrm {d}k_{y} \,\mathrm {d}k_{z}, $$
(26)

where \(\mathbf{v}_{\perp}\) is the projection of the velocity perpendicular to the x-direction and we introduced the partial Fourier transform of \(\mathbf{v}_{\perp}\):

$$ \begin{aligned}[b] \tilde{\mathbf{v}}_{\perp}(k_{y}, k_{z}) = {}& \lim_{L_{x}\to+\infty} \frac{1}{L_{x}} \int ^{L_{x}/2}_{-L_{x}/2} \,\mathrm {d}x \\ &{}\times \int_{\mathbb{R}^{2}} e^{-2\pi \mathrm {i}(k_{y} y + k_{z} z)} \mathbf{v}_{\perp}(x,y,z) \,\mathrm {d}y \,\mathrm {d}z. \end{aligned} $$
(27)

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_{\perp}\) 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.

3 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 1283 grid. For all of the runs the timestep is chosen to have a CFL, i.e., \(c \Delta t / \Delta x\), of 0.4, c 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 (Toro 1999). The runs with PPM use the original flattening and artificial viscosity prescriptions from Colella and Woodward (1984). 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 (Porter and Woodward 1994), 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: 643, 1283, 2563 and 5123. The RMS velocity in all of the runs is \(v_{\mathrm{rms}} \simeq0.4\), giving an eddy turnover time \(\tau= 1/v_{\mathrm{rms}} \simeq2.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\gtrsim3\) (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 \leq t \leq100\).

Figure 1
figure 1

Time evolution of the diagnostic quantities for the fiducial set of runs PPM_HLLC with different resolutions. The left panel shows the root mean square (RMS) Mach number, while the middle and right panels show, respectively, the ratios \(R_{xx}/R_{yy}\) and \(R_{xx}/R_{zz}\), R being the Reynolds stress tensor (equation (4)). Since the x-direction is the anisotropic direction (it would play the role of the radial direction in a CCSN) the ratios \(R_{xx}/R_{yy}\) and \(R_{xx}/R_{zz}\), offer a global measure of the anisotropy of the flow at the largest scale. All of the quantities appear to have reached stationarity after time \(t \gtrsim3\) and oscillate around their target values. All resolutions produce the same qualitative behavior.

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 (5123). 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 (Vincent and Meneguzzi 1991; Ishihara et al. 2007). 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.

Figure 2
figure 2

Square root of the magnitude of the vorticity, \(\pmb{\sqrt{\vert \nabla\times\mathbf{v}\vert }}\) , for four of the simulations with 512 \(\pmb{^{3}}\) resolution in a slice through the middle of the x - z plane at the final time of the simulations ( \(\pmb{t=100}\) ). The panels show simulations using PPM_HLLE_N512, PPM_HLLC_N512, WENOZ_HLLC_N512, and TVD_HLLC_N512 clockwise from the top left. The direction of the anisotropic driving is up in these figures. The colorcode goes linearly from 0 (no vorticity; dark colors) to 15 (light colors) and it is the same for all panels.

4 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.

4.1 Energy decay rate

In the limit of very large Reynolds number it is assumed, in standard turbulence phenomenology (Frisch 1996), 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):

$$ \partial_{t} E(k) + \partial_{k} \Pi(k) \simeq0 , $$
(28)

so that stationarity requires \(\Pi(k) \simeq\mathrm{const}\). In particular, since energy is conserved, one finds \(\Pi(k) \simeq\langle\epsilon\rangle \). 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 (Kaneda et al. 2003). 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 (Couch and Ott 2015), 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.

Figure 3
figure 3

Energy flux, as defined by equation ( 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 5123. 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.

As discussed before, we expect that \(\Pi(k) \simeq\langle\epsilon \rangle\) 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 (Sytine et al. 2000), we find that regions where \(\Pi\simeq\langle\epsilon\rangle\) as wide as a few wave numbers \(4 \lesssim k \lesssim10\) only appear at the highest resolutions (we will discuss the inertial range in more detail in Section 4.2). 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 1283 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 4.3).

Figure 4
figure 4

Dissipation rate of the energy at the largest scales due to the turbulent cascade (not including direct dissipation by the numerical viscosity) as a function of resolution and for all of the schemes. The dissipation rate is normalized so as to be 1 in the limit of large Reynolds numbers/resolution. At 1283 points all of the schemes show an error of less than 10%, with the HLLC schemes already close to the 2% level.

4.2 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 \Delta 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 \Delta x\) coincide with the dimensional one for the highest resolution runs. With this choice, \(512 k \Delta x=512\) corresponds to a wavelength of a single grid point, \(512 k \Delta x=256\) corresponds to a wavelength of two grid points and so on.

Figure 5
figure 5

Energy spectra (equation ( 10 )) obtained with different numerical methods and resolutions. The energy spectra are compensated by a \(k^{5/3}\) spectrum, so that any region with Kolmogorov scaling should appear roughly flat. Furthermore, the spectra are all plotted as a function of the dimensionless wave number \(512 k \Delta x\) (the 512 factor is introduced to have the dimensionless wave number coincide with the dimensional one for the 5123 runs). The first five panels show the PPM_HLLC (upper left), PPM_HLLC_CFL0.8 (upper center), PPM_HLLE (upper right), TVD_HLLE (lower left) and WENOZ_HLLC (lower center) group of runs. The last panel (lower right) shows a comparison of all of the methods at the highest resolution (5123). An inertial range seems to be recovered only at the highest resolutions (perhaps with the exception of TVD_HLLE where no inertial range is visible). All schemes employing the HLLC Riemann solver are in very good agreement.

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 (Frisch 1996) (cf. Gotoh et al. (2002)). 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 (Aspden et al. 2009). 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 \(\gtrsim8 \Delta x\) (\(512 k \Delta x \lesssim64\)), 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 \Delta x = 10\) for the HLLC runs), or until the energy injection scale (\(512 k \Delta 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 (Abdikamalov et al. 2015). A quantification of the bottleneck effect in terms of the energy budget is discussed in Section 4.4.

At even larger scales, an inertial range (\(E\sim k^{-5/3}\) and \(\Pi \sim \mathrm{const}\), see Figure 3) seems to be recovered by the least dissipative schemes (PPM and WENOZ with HLLC) in the region \(4\lesssim k \lesssim10\). 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 (Kaneda et al. 2003; Federrath 2013), but it would already correspond to an extremely high resolution in global CCSN simulations that typically have 30-70 zones across the turbulent region (Abdikamalov et al. 2015).

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.

4.3 Numerical 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 \(\nu(k)\) (equation (19)). The latter is shown in Figure 6 for all schemes and resolutions.

Figure 6
figure 6

Numerical viscosity as a function of the wave number measured for all schemes and resolutions. The numerical viscosity is estimated using the procedure outlined in Section  2 and it is defined by equation ( 19 ). The different panels are, from left to right and from top to bottom, 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 5123. The numerical viscosity shows large variations across the wave number space. The choice of the Riemann solver plays a role that is at least as important as the choice of the reconstruction method in affecting the numerical viscosity throughout the entire the spectrum.

The first thing to notice is that the numerical 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 effective Reynolds number achieved in a simulation ambiguous. Meaningful ways to estimate it for ILES have been proposed (Zhou et al. 2014) 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, because the dissipative properties of ILES differ considerably from the ones of the true Navier-Stokes equations.

Two other features can be observed in most of the numerical 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 numerical viscosity appears small because of a partial decoupling between vorticity and dissipation. Second, at high resolution and at the largest scales, the numerical 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 numerical 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 appears to have an effect comparable to the difference between the PPM and WENOZ reconstructions at intermediate scales (\(40\lesssim k \lesssim100\)).

4.4 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 4.1), and with the energy transfer in the inertial range, which we have seen to be described accurately only at much higher resolutions (Section 4.2). 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:

$$ \mathcal{E}(k) = \int_{0}^{k} E( \xi) \,\mathrm {d}\xi. $$
(29)

This quantity is plotted in Figure 7, where it is normalized by

$$ \frac{v_{\mathrm{rms}}^{2}}{2} = \int_{0}^{+\infty} E(k) \,\mathrm {d}k $$
(30)

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:

$$\begin{aligned}& \mathcal{E}_{\mathrm{K41}}(k) = \int_{0}^{k} E_{\mathrm{K41}}(\xi) \,\mathrm {d}\xi, \end{aligned}$$
(31)
$$\begin{aligned}& \begin{aligned}[b] &E_{\mathrm{K41}}(k) \\ &\quad = \textstyle\begin{cases} E_{\mathtt{PPM\_HLLC\_N512}}(k), & \text{if } k \leq4 , \\ E_{\mathtt{PPM\_HLLC\_N512}}(4) (\frac{k}{4} )^{-5/3}, & k > 4 . \end{cases}\displaystyle \end{aligned} \end{aligned}$$
(32)

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 643 only \({\sim}50\%\) or less of the kinetic energy is in well resolved structures, while the other \({\sim}50\%\) have piled up at rather large scale, with a cumulative excess of \({\sim}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 5123 this is still only about 80% of the energy and there is still a cumulative excess of over \({\sim}5\%\) at the grid scale (\(\ell\sim \Delta x\)).

Figure 7
figure 7

Cumulated energy distribution (equation ( 29 )) for all methods and resolutions, normalized by a factor \(\pmb{2/v_{\mathrm{rms}}^{2}}\) to be equal to 1 for large k . As a reference for comparison we also plot the asymptotic profile expected from Kolmogorov’s theory (equation (31)). The different panels are, from left to right and from top to bottom, 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 5123. At low resolution all of the schemes show an excess of energy at intermediate scales, due to the bottleneck. Only at the highest resolution at least, roughly, 80% of the energy is correctly resolved.

5 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 (Porter et al. 2002) and driven (Benzi et al. 2008) turbulence. In the anisotropic case, however, anisotropic contributions cannot be excluded (Biferale et al. 2002), although they are known to be subdominant in some important cases (Calzavarini et al. 2002; Biferale et al. 2003; Kaneda et al. 2008). In this section we show that equation (25) is consistent with our data over a wide range of scales.

We compute the 3rd-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 3rd 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 \(-\frac{5}{4} r^{-1} \langle\epsilon\rangle^{-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 \(\eta\sim\Delta x\) and plot the structure functions versus \(r/\Delta x\).

Figure 8
figure 8

Compensated 3rd-order structure functions (equation ( 23 )) for all the numerical methods and resolutions. The structure functions are compensated and scaled so that they should be close to one where the \(4/5\)-law (equation (25)) is verified. The data is plotted as a function of the dimensionless separation \(r/\Delta x\). The first five panels show the PPM_HLLC (upper left), PPM_HLLC_CFL0.8 (upper center), PPM_HLLE (upper right), TVD_HLLE (lower left) and WENOZ_HLLC (lower center) group of runs. The last panel (lower right) shows a comparison of all of the methods at the highest resolution (5123). The \(4/5\)-law is very well verified in our data suggesting that (1) anisotropic corrections are subdominant and (2) all of the simulations behave in a way consistent with large Reynolds numbers turbulence at the largest scales.

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. (2008) and has been known to be true also for Rayleigh-Bénard convection in most regimes (Lohse and Xia 2010). Our results provide an important new example where this appears to hold true. Secondly, for all of our simulations at 5123, we find

$$ \max_{r} \biggl\{ -\frac{5}{4} r^{-1} \langle \epsilon\rangle^{-1} \bigl\langle \delta v^{3} \bigr\rangle _{j=0} \biggr\} $$
(33)

within 5% of 1. This level of accuracy is reached in DNS simulations achieving at least a Taylor micro-scale Reynolds number (Kaneda et al. 2008)

$$ R_{\lambda}= \frac{u' \lambda}{\nu} \sim300 , $$
(34)

where ν is the kinematic viscosity, \(u' = \frac{1}{\sqrt{3}} v_{\mathrm{rms}}\) and \(\lambda= (15 \nu u^{\prime2}/\langle\epsilon\rangle )^{1/2}\) is the Taylor micro-scale. This corresponds to a large-scale Reynold numbers \(R = \frac{u' L}{\nu} \sim R_{\lambda}^{2}\), \(L=1\) being the domain size, in excess of 85,000. Reaching these Reynolds numbers in a DNS requires resolutions between 5123 and 10243 using pseudo-spectral methods (Donzis et al. 2008). This large-scale estimate of the Reynolds number is consistent with previous findings (Zhou et al. 2014), although it is several orders of magnitude larger than the one that could be naively estimated using \(\nu_{\max}\). For instance, for PPM_HLLC at 5123, \(\nu_{\max} \simeq1.5 \times 10^{-3}\) and \(v_{\mathrm{rms}} \simeq0.4\) giving

$$ \frac{u' L}{\nu_{\max}} \simeq150 . $$
(35)

This apparent discrepancy is due to the fact that an ILES is not a DNS. 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 had multiple values of the Reynolds number.

6 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 (Dolence et al. 2013; Couch and O’Connor 2014; Handy et al. 2014; Abdikamalov et al. 2015).

Figure 9 shows a comparison of the 2D transverse spectrum \(E_{\perp}(k_{\perp})\) 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_{\perp}(k_{\perp})\) only shows a very small region that suggests an inertial range, \(3 \lesssim k \lesssim5\) (as opposed to \(5 \lesssim k \lesssim 10\) in \(E(k)\)).

Figure 9
figure 9

3D (equation ( 10 ), blue) and 2D, transverse (equation ( 26 ), red) energy spectra for the PPM_HLLC simulations. The energy spectra are compensated by \(k^{5/3}\) to highlight eventual regions with Kolmogorov scaling. The spectra are plotted as a function of the dimensionless wave number \(512 k \Delta x\), as in Figure 5. Although \(E(k)\) and \(E_{\perp}(k_{\perp})\) have similar trends, the use of the transverse spectrum can overestimate the width of the bottleneck region.

Abdikamalov et al. (2015) 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 3rd-order structure functions, that, as we have shown, can give a more accurate description of the energy cascade.

7 Conclusions

The details of the explosion mechanism of CCSNe have eluded our comprehension in spite of more than 50 years of studies (Woosley and Janka 2005; Janka et al. 2012; Burrows 2013; Foglizzo et al. 2015). Recent numerical advances (Murphy et al. 2013; Couch and Ott 2013; Couch and Ott 2015; Müller and Janka 2015) suggest that turbulence might play a fundamental role in tipping over the balance of the forces and lead to successful explosions (see also Appendix). At the same time, the level of accuracy of current simulations, which employ the ILES methodology, is unclear (Abdikamalov et al. 2015). Turbulence in CCSNe is mildly compressible, but strongly anisotropic (Murphy et al. 2013; Couch and Ott 2015). 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 643 to 5123. 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 3rd-order structure functions). Finally, we measured the numerical 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 1283 (e.g., \({\sim}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 (5123). Even at this resolution, which would be difficult to achieve in global simulations, only roughly \({\sim}80\%\) (the exact number depends on the scheme, see Section 4.4) of the energy is resolved, the remaining \({\sim}20\%\) accumulates as excess energy at intermediate scales (the cumulative energy excess at the grid scale alone is as large as \({\sim}5\%\) of the total energy).

Current CCSN simulations have resolutions of at most of 30-70 points covering the gain region (Abdikamalov et al. 2015) (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 4.1). 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. (2015), 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. (2015), 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 \({\sim}20 \mathrm{~ms}\) in a CCSN (Ott et al. 2013; Couch and O’Connor 2014)). 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) (Couch and Ott 2015). 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 due 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 3rd 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. (2006); Martín et al. (2006) or Thornber et al. (2008). 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.

References

Download references

Acknowledgements

We acknowledge helpful discussions with E Abdikamalov, WD Arnett, A Burrows, R Fisher, C Meakin, P Mösta, J Murphy, M Norman, and L Roberts. This research was partially supported by the National Science Foundation under award nos. AST-1212170 and PHY-1151197 and by the Sherman Fairchild Foundation. The simulations were performed on the Caltech compute cluster Zwicky (NSF MRI-R2 award no. PHY-0960291), on the NSF XSEDE network under allocation TG-PHY100033, and on NSF/NCSA BlueWaters under NSF PRAC award no. ACI-1440083. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to David Radice.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

DR ran the simulations, performed the analysis of the data and wrote the basic draft of this paper. SMC assisted with the use of the FLASH code. CDO had the original idea that started this investigation. All of the authors took part in discussions concerning the results and contributed corrections and improvements on the early draft of the manuscript. All authors read and approved the final manuscript.

Appendix: The role of turbulence in core-collapse supernova explosions

Appendix: The role of turbulence in core-collapse supernova explosions

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. (2012) and Burrows (2013). Next, we will discuss the role of turbulence and, in particular, of turbulent pressure on the explosion mechanism, in light of some recent results (Murphy et al. 2013; Couch and Ott 2015).

1.1 A.1 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 (\({\sim}4\mbox{-}5 \times10^{14} \mbox{ 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 neutrino-semitransparent regions behind the shock (Bethe 1990). 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 \({\sim}100\mbox{-}200 \mathrm{~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., O’Connor and Ott (2011)).

During this time, however, the PNS will release a significant fraction of its binding energy in the form of neutrinos (of order \(10^{46} \mathrm{~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 (Bethe 1990; Pejcha and Thompson 2012). 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 Couch and O’Connor (2014)). 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.

1.2 A.2 Turbulent pressure and the Rankine-Hugoniot conditions

Simulations (Burrows et al. 1995; Murphy et al. 2013; Couch and Ott 2015) 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,

$$ s \bigl[ \rho_{d} v_{d}^{r} - \rho_{u} v_{u}^{r} \bigr] = \rho_{d} \bigl(v_{d}^{r}\bigr)^{2} + p_{d} - \rho_{u} \bigl(v_{u}^{r}\bigr)^{2} - p_{u}, $$
(36)

where s is the shock speed and \(\cdot_{d}\) and \(\cdot_{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:

$$\begin{aligned} p_{u} = 0 \qquad \bigl(v_{u}^{r}\bigr)^{2} = \frac{GM}{r}. \end{aligned}$$
(37)

For the shock to expand we must then have

$$ \rho_{d} \bigl(v_{d}^{r} \bigr)^{2} + p_{d} > \rho_{u} \frac{GM}{r}. $$
(38)

In the presence of turbulence, Murphy et al. (2013) suggested to modify equation (38) in a way akin to a Reynolds decomposition and write it as

$$ \rho_{d} \bigl(\bar{v}_{d}^{r} \bigr)^{2} + \rho_{d} \bigl(\delta v_{d}^{r} \bigr)^{2} + p_{d} > \rho_{u} \frac{GM}{r}, $$
(39)

where v is the average velocity and \(\delta\mathbf{v} = \bar{\mathbf{v}} - \mathbf{v}\) is the turbulent velocity. Although not entirely rigorous, equation (39) has been shown to be well verified in the numerical simulations if angular averages are used to compute the respective quantities (Murphy et al. 2013; Couch and Ott 2015). Couch and Ott (2015) 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 (39).

Going beyond the arguments of Murphy et al. (2013), we can reinterpret equation (39) 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 \(\delta v_{d}^{r}\) in terms of the specific turbulent energy

$$ e_{\mathrm{turb}} = \frac{1}{2} \vert \delta\mathbf{v} \vert ^{2} $$
(40)

noting that

$$ \vert \delta\mathbf{v}\vert ^{2} := \bigl(\delta v^{r} \bigr)^{2} + \bigl(\delta v^{\theta}\bigr)^{2} + \bigl( \delta v^{\phi}\bigr)^{2}, $$
(41)

and using the fact that

$$ \bigl(\delta v^{r}\bigr)^{2} \simeq\bigl( \delta v^{\theta}\bigr)^{2} + \bigl(\delta v^{\phi}\bigr)^{2} $$
(42)

in CCSN turbulence, to obtain

$$ \bigl(\delta v^{r}\bigr)^{2} \simeq \frac{1}{2} \vert \delta\mathbf{v} \vert ^{2} = e_{\mathrm{turb}}. $$
(43)

Assuming the pressure varies like \(p = (\gamma- 1) \rho e\), and substituting (43) into (39), we find

$$ \begin{aligned}[b] &\rho_{d} \bigl(\bar{v}_{d}^{r} \bigr)^{2} + (\gamma_{\mathrm{th}} - 1) \rho_{d} e_{d} + (\gamma_{\mathrm{turb}} - 1) \rho_{d} e_{\mathrm{turb}} \\ &\quad > \rho_{u} \frac{GM}{r}, \end{aligned} $$
(44)

where \(\gamma_{\mathrm{th}} \simeq4/3\) is the thermodynamical adiabatic index, \(e_{d}\) is the downstream thermal energy and \(\gamma_{\mathrm{turb}} = 2\) is the equivalent adiabatic index associated with anisotropic CCSN turbulence. Since \(\gamma_{\mathrm{turb}} > \gamma_{\mathrm{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 (42) is dropped and turbulence is assumed to be isotropic, then \(\gamma_{\mathrm{turb}} = 5/3\), which is still larger than \(\gamma_{\mathrm{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_{\mathrm{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 Müller and Janka (2015) who formulated their equations in terms of the turbulent Mach number, as opposed to the turbulent energy.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Radice, D., Couch, S.M. & Ott, C.D. Implicit large eddy simulations of anisotropic weakly compressible turbulence with application to core-collapse supernovae. Comput. Astrophys. 2, 7 (2015). https://doi.org/10.1186/s40668-015-0011-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40668-015-0011-0

Keywords