# The non-linear onset of neutrino-driven convection in two and three-dimensional core-collapse supernovae

Rémi. Kazeroni, Brendan.K. Krueger, Jérôme Guilet, Thierry Foglizzo, Daniel Pomarède

► To cite this version:

Rémi. Kazeroni, Brendan.K. Krueger, Jérôme Guilet, Thierry Foglizzo, Daniel Pomarède. The non-linear onset of neutrino-driven convection in two and three-dimensional core-collapse supernovae. Monthly Notices of the Royal Astronomical Society, 2018, 480 (1), pp.261-280. 10.1093/mn-ras/sty1742. hal-01724901

HAL Id: hal-01724901

https://hal.science/hal-01724901v1

Submitted on 15 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The non-linear onset of neutrino-driven convection in two- and three-dimensional core-collapse supernovae

Rémi Kazeroni,<sup>1,2★</sup> Brendan K. Krueger,<sup>2,3</sup> Jérôme Guilet,<sup>2,1★</sup> Thierry Foglizzo<sup>2★</sup> and Daniel Pomarède<sup>4</sup>

Accepted 2018 June 28. Received 2018 June 15; in original form 2018 February 22

## ABSTRACT

A toy model of the post-shock region of core-collapse supernovae is used to study the non-linear development of turbulent motions driven by convection in the presence of advection. Our numerical simulations indicate that buoyant perturbations of density are able to trigger self-sustained convection only when the instability is not linearly stabilized by advection. Large amplitude perturbations produced by strong shock oscillations or combustion inhomogeneities before the collapse of the progenitor are efficiently shredded through phase mixing and generate a turbulent cascade. Our model enables us to investigate several physical arguments that had been proposed to explain the impact of the dimensionality on the onset of explosions in global simulations of core-collapse supernovae. Three-dimensional (3D) simulations are found to lead to higher entropy values than two-dimensional (2D) ones. We attribute this to greater turbulent mixing and dissipation of the kinetic energy into heat in 3D. Our results show that the increase of entropy is enhanced with finer numerical resolution and larger perturbation amplitude.

**Key words:** accretion, accretion discs – hydrodynamics – instabilities – shock waves – turbulence – supernovae: general.

## 1 INTRODUCTION

During the first second after the core of a massive star collapses into a proto-neutron star, the multidimensional hydrodynamics of the innermost region of few hundreds of kilometres plays a crucial role in driving an explosion of the star, known as a core-collapse supernova (CCSN). The delayed neutrino mechanism (Bethe & Wilson 1985) represents the most favoured mechanism to explain a large majority of CCSNe (see Foglizzo et al. 2015; Janka, Melson & Summa 2016; Müller 2016 for recent reviews). The shock wave formed at core bounce turns into an accretion shock and stalls in the iron core at a radius of about 150 km due to energy loss induced by neutrino emission and dissociation of heavy nuclei. Neutrinos carry away most of the energy (~10⁵<fnref n="3" /> erg) released in the gravitational contraction of the core to a proto-neutron star. A small fraction of the neutrino flux is absorbed in the gain layer, a region just below the stalled shock wave where neutrino absorption dominates neutrino emission. A successful explosion is powered and a neutron star is

left behind the CCSN if the energy deposition is efficient enough to revive the stalled shock wave, preventing the formation of a black hole (O'Connor & Ott 2011).

Self-consistent numerical simulations of CCSNe including, among others, detailed neutrino transport and microphysics, three-dimensional (3D) hydrodynamics and general relativity still remain computationally challenging. Approximations in treatment of the physics and the numerics, such as a reduction of the dimensionality, are required to gain insight in the CCSN problem. Numerical calculations performed in 1D demonstrated that spherical symmetry prevents the onset of explosions (Liebendörfer et al. 2001) except for the lowest CCSN progenitor masses (Kitaura, Janka & Hillebrandt 2006). A successful explosion relies on the decisive action of multidimensional hydrodynamical instabilities which break the spherical symmetry of the collapse and generate vigorous non-radial motions to large spatial scales. This increases the advection time of matter through the gain region, enhances the heating efficiency and creates more favourable conditions to initiate a CCSN (e.g. Murphy & Burrows 2008). Axisymmetric simulations (2D) with various approximations established that two distinct instabilities can govern the dynamics of the gain region prior to the shock revival (Müller, Janka & Heger 2012; Fernández et al. 2014). These results were later

<footnote n="1"><sup>1</sup>Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany</footnote>

<footnote n="2"><sup>2</sup>IRFU/Département d'Astrophysique, Laboratoire AIM, CEA/DRF-CNRS-Université Paris Diderot, CEA-Saclay F-91191, France</footnote>

<footnote n="3"><sup>3</sup>Eulerian Codes (XCP-2), Los Alamos National Laboratory, Los Alamos, NM 87545, USA</footnote>

<footnote n="4"><sup>4</sup>Institut de Recherche sur les Lois Fondamentales de l'Univers, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France</footnote>

<footnote>* E-mail: kazeroni@MPA-Garching.MPG.DE (RK); jerome.guilet@cea.fr (JG); thierry.foglizzo@cea.fr (TF)</footnote>
confirmed in 3D (Takiwaki, Kotake & Suwa 2014; Melson, Janka & Marek 2015a; Melson et al. 2015b; Lentz et al. 2015; Roberts et al. 2016). Neutrino-driven convection was seen in some of the earliest multidimensional CCSNe simulations (Herant, Benz & Colgate 1992; Herant et al. 1994; Burrows, Hayes & Fryxell 1995; Janka & Mueller 1996). The negative entropy gradient induced by neutrino heating can be unstable to buoyancy and convective overturns develop in the gain region, characterized by intermediate spatial scales $l \sim 5-6$. The second instability is the standing accretion shock instability (SASI, Blondin, Mezzacappa & DeMarino 2003) mediated by an advective-acoustic cycle (Foglizzo et al. 2007; Foglizzo 2009; Guilet & Foglizzo 2012). It generates global oscillatory motions of the shock wave associated with the largest spatial scales $l \sim 1-2$.

Large-scale asymmetries in the collapsing core are supported by some observational evidence. The development of early asymmetries is suggested by the spectropolarimetric measurement of a type II-P supernova during the nebular phase of the explosion (Leonard et al. 2006). High kicks of young neutron stars (Arzoumanian, Chernoff & Cordes 2002) could result from a large-scale asymmetry $l = 1$ present at the onset of the explosion (Scheck et al. 2004, 2006). The one-sided spatial distribution of $^{44}$Ti observed in Cassiopea A (Grefenstette et al. 2014) seems to confirm numerical calculations that predict large-scale asymmetries seeded by hydrodynamical instabilities (Wongwathanarat, Janka & Müller 2013; Wongwathanarat et al. 2017). Convection and SASI could leave different imprints on the gravitational waves (Müller, Janka & Marek 2013; Kuroda, Takiwaki & Kotake 2016; Andresen et al. 2017) and the neutrino signals (Lund et al. 2012; Tamborra et al. 2013, 2014a) that might be detected for a future nearby CCSN.

A majority of 2D simulations produce explosions that are underenergetic compared to the range of explosion energies expected from observations (Müller 2015). It remains unclear whether releasing the axisymmetric constraint reduces the discrepancies, in part because running self-consistent 3D simulations long enough to infer the explosion energy is barely feasible. Moreover, the impact of the dimensionality on the likelihood to explode is still controversial. Models assuming a light-bulb approximation for the neutrino transport reached different conclusions regarding the critical neutrino luminosity required to obtain an explosion. Burrows, Dolence & Murphy (2012) and Dolence et al. (2013) found slightly more favourable conditions in 3D, while Hanke et al. (2012) saw almost no difference and Couch (2013) obtained delayed explosions in 3D. More accurate neutrino transport schemes provided less ambiguous results since the shock revival either fails in 3D (Hanke et al. 2013; Tamborra et al. 2014b) or is notably delayed compared to the corresponding 2D simulation (Takiwaki et al. 2014; Melson et al. 2015b; Lentz et al. 2015). The first 3D simulations that produced more robust explosions associated with a faster growth of the energy (Melson et al. 2015a; Müller 2015) raised hope that 3D is not necessarily detrimental to CCSNe. Besides, the inclusion of multidimensional initial conditions obtained by simulating shell burning a few minutes prior to collapse (Müller et al. 2016) may lead to an explosion energy in the ballpark of observations (Müller et al. 2017).

Several physical arguments have been put forward to interpret the discrepancies between 2D and 3D simulations. The turbulent energy cascade acts in different directions in 2D and 3D (Hanke et al. 2012). In 2D, the reverse turbulent cascade feeds the largest spatial scales (Kraichnan 1967) and favours the formation of large buoyant bubbles that help to trigger shock revival (Fernández et al. 2014). On the contrary, the forward cascade in 3D transfers energy to small scales and induces more dissipation of turbulent kinetic

energy. Turbulence created by the instabilities in the post-shock region provides additional pressure support to revive the shock wave. On the one hand, it was proposed that turbulent kinetic energy, related to the largest spatial scales, is artificially overpredicted in 2D simulations (Murphy, Dolence & Burrows 2013; Couch & Ott 2015) and this could lead to earlier explosions. On the other hand, turbulent dissipation to small scales could play a positive role in 3D (Mabanta & Murphy 2018), but this has to be tested in multidimensional simulations. It was also suggested that 2D buoyant bubbles undergo a weaker drag force due to their particular shapes and drive earlier explosions (Couch 2013). In their simulations of low-mass progenitors dominated by the convective instability, Melson et al. (2015a) and Müller (2015) observed that positive 3D effects could come into play after shock revival. Greater fragmentation and mixing to small scales at the edge of downflows are able to decelerate them, keeping more matter in the gain region and fostering explosions in 3D. Whether these effects are generic and play a role in more massive progenitors is not yet known.

Using analytical calculations, Foglizzo, Scheck & Janka (2006) showed that convection can be stabilized by a fast advection through the gain region. Therefore, a negative entropy gradient is not a sufficient instability criterion because the exchange of matter induced by buoyancy should take place before the perturbation leaves the gain region. The linear analysis of Foglizzo et al. (2006) is strictly speaking only valid for perturbations of small amplitudes. It has been proposed that strong perturbations could possibly trigger convection in situations where it is linearly stabilized by advection (Scheck et al. 2008; Fernández et al. 2014). SASI acts as a source of perturbations because the instability pushes the shock outward and creates entropy-vorticity waves which can feed secondary convection (Scheck et al. 2008; Guilet, Sato & Foglizzo 2010; Cardall & Budiardja 2015; Summa et al. 2016). Progenitor inhomogeneities produced by convective shell burning create multidimensional initial conditions which can ease the shock revival under certain circumstances (Couch & Ott 2013, 2015; Müller & Janka 2015). It was shown that this mechanism could lead to the formation of large-scale structures ($l \sim 2-4$) which can turn a failed explosion into a successful one (Müller et al. 2016, 2017). Lastly, the noise related to some numerical features such as the grid geometry may also seed convection. A Cartesian grid with several refinement levels towards the centre of the domain generates greater perturbations than a spherical grid. Such perturbations alone could impact the dynamics to the point of selecting the dominant instability during the whole simulation. A $27 M_{\odot}$ progenitor was found to be dominated by convection (Ott et al. 2013) or by SASI after overcoming the grid artefacts (Abdikamalov et al. 2015).

These numerous works have explored different aspects of the multidimensional dynamics and yet several key questions remain uncertain. What is the main physical effect that drives multidimensional simulations closer to explosions? Several effects have been put forward such as a longer residency time-scale in the gain layer (Murphy & Burrows 2008), the turbulent pressure pushing the shock (Couch & Ott 2015) or heating by dissipation of kinetic energy (Mabanta & Murphy 2018), but their relative importance remain to be determined. What are the main differences between 2D and 3D? What is the role of perturbations resulting from the pre-collapse dynamics or seeded by numerical errors? Can they trigger self-sustained convection in the regime where advection linearly stabilizes the instability? Addressing these issues and disentangling the different effects in global CCSN simulations is challenging, such that a simplified setup is an important complementary approach. The simplicity of such models is an asset not only
for the physical interpretation but also to allow an easier exploration of the parameter space and to study the impact of numerical resolution.

In this study, we investigate the competition between advection and convection in its simplest form. An idealized description of the advection of matter through the gain layer during the stalled shock phase of a CCSN is employed. Numerical simulations of the toy model enable us to assess the robustness of linear and non-linear instability criteria. We show that a strong perturbation is not a sufficient condition to lead to fully developed convection and then discuss several regimes of the instability. Various interpretations of the impact of dimensionality on convection are revisited in the light of our model that reproduces the classical properties of turbulence induced by the instability. Our simulations show an excess of heating in 3D due to the turbulent dissipation of kinetic energy. We identify the hydrodynamical processes that could explain this phenomenon based on the particular dynamics at small spatial scales. We discuss the influence of the numerical resolution and show that a greater impact on the dynamics is expected when large amplitude perturbations come into play.

The rest of the paper is organized as follows. The physical model and the numerical methods are detailed in Section 2. Predictions regarding the linear and non-linear onsets of the convective instability are tested in Section 3. An overview of the main discrepancies between the 2D and the 3D dynamics is presented in Section 4. A quantitative analysis of the properties of turbulence is proposed in Section 5 to pinpoint the effects that could explain the main differences between our 2D and 3D simulations. The influence of the numerical resolution is discussed in Section 6. The limitations of our model and the consequences of our results for self-consistent simulations are addressed in Section 7. We summarize our findings and conclude in Section 8.

## 2 PHYSICAL AND NUMERICAL SETUP

Our approach focuses on an idealized description of the accretion flow below the stalled shock wave of a CCSN and includes only the minimal ingredients leading to the convective instability. This phase is modelled by a stationary flow crossing a source layer which represents a fraction of a realistic gain layer and contains gravity and neutrino heating. The shock wave and the cooling layer surrounding the proto-neutron star are absent from our model in order to study convection in its simplest form and preclude any feedback that would interact with the dynamics of the gain region.

### 2.1 Stationary flow

The stationary flow considered in this study is similar to the one studied by Foglizzo et al. (2006). It is planar along the vertical direction (z). The flow is modelled as a perfect gas with an adiabatic index $\gamma = 4/3$. This is a suitable approximation of the equation of state in the post-shock region where the pressure is dominated by relativistic electrons, photons, and electron-positron pairs (Janka 2001).

In our model, the gain region contains the effects of heating and gravity and is located in the central part of the domain. Upstream and downstream from this layer, gravity and heating are turned off and the stationary flow is thus uniform. The two external layers are used to place the boundaries far enough in the vertical direction in order to minimize the impact of numerical reflections on the instability. To avoid any discontinuity at the horizontal edges of the source layer, gravity and heating are smoothed out by a linear ramp

function $\Psi(z)$:

$$ \Psi(z) = \begin{cases} 1 & \text{if } |z| < \frac{H}{2}, \\ 2\left(1 - \frac{z}{H}\right) & \text{if } \frac{H}{2} \le |z| < H, \\ 0 & \text{if } H \le |z|, \end{cases} \qquad (1) $$

where $2H$ corresponds to the height of the gain layer. To ease the comparison with more realistic models, $H$ is set to 50 km. Variations of $H$ due to shock oscillations and the recession of the gain radius over time are neglected for the sake of simplicity. The numerical domain is delimited in the horizontal directions by $-150 \text{ km} \le x, y \le 150 \text{ km}$ and in the vertical direction by $-450 \text{ km} \le z \le 450 \text{ km}$. The horizontal extent of our domain covers about one-third of the angular size of a gain region for a shock radius of 150 km. The gain layer is delimited vertically by $-50 \text{ km} \le z \le 50 \text{ km}$.

The gravitational potential is defined as:

$$ \nabla \Phi = K_G \left( \frac{c_{up}^2}{H} \right) \Psi(z), \qquad (2) $$

where $K_G$ is a dimensionless parameter that quantifies the intensity of gravity, $c$ is the sound speed, and the subscript 'up' refers to quantities upstream from the source layer. The energy deposition by neutrinos per unit of time in the source layer is modelled by a heating function of the form:

$$ \mathcal{L} = K_H \left( \frac{\rho_{up} M_{up} c_{up}^3}{\gamma H} \right) \left( \frac{\rho}{\rho_{up}} \right) \Psi(z), \qquad (3) $$

where $\rho$ corresponds to the density and $M$ to the Mach number defined as: $M = |v_z|/c$ with $v$ being the fluid velocity. The dimensionless numbers $K_G$, $K_H$, and $M_{up}$ represent the free parameters of our model. In the numerical simulations of the model, we consider only variations of $K_H$, the normalization of the heating function. This parameter controls whether convection is linearly unstable or stabilized by advection (see Section 2.2 and Appendix B).

The parameter $M_{up}$ is the equivalent to the post-shock Mach number. Its value depends on the equation of state and the description of the matter. For an ideal gas with an adiabatic index $\gamma = 4/3$, it varies from $M_{up} \approx 0.1$ when nuclei are completely dissociated by the shock wave to $M_{up} \approx 0.3$ for an adiabatic shock (Fernández & Thompson 2009). In this study, we consider the case $M_{up} = 0.3$.

An estimate of the normalization of the gravitational potential can be made using $H$ and $c_{up}$. The latter can be computed as

$$ c_{up}^2 = \frac{v_{up}^2}{M_{up}^2} = \frac{v_{ff}^2}{M_{up}^2 \kappa^2} = \frac{2GM_{NS}}{r_{sh}} \frac{1}{M_{up}^2 \kappa^2}, \qquad (4) $$

where $\kappa$ corresponds to the compression factor across the shock and $v_{ff}$ to the free-fall velocity at the shock radius. The parameter $\kappa$ ranges from about 5 for an adiabatic shock to about 10 for a strong shock that fully dissociates iron nuclei. Using equations (2) and (4), one obtains:

$$ K_G = \nabla \Phi \frac{H}{c_{up}^2} = \left( \frac{r_{sh}}{r} \right)^2 \frac{H}{r_{sh}} \frac{M_{up}^2 \kappa^2}{2}. \qquad (5) $$

During the stalled shock phase, one typically finds that $r_{sh}/r \sim 1.4-3$ and $H/r_{sh} \sim 0.3-0.7$ (e.g. Marek & Janka 2009; Melson et al. 2015a). Thus, the normalization of the gravitational potential is such that $K_G \sim 0.6-5.4$ for an adiabatic shock and $K_G \sim 0.3-2.7$ for a strong shock. In the following, we restrict ourselves to $K_G = 3$.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Vertical profiles of gravity (red), vertical velocity (green), heating (blue), and entropy (orange) across the gain region. The quantities are normalized by their maximum values except the entropy which is shown as $(S - S_0)/S_{max}$, where $S_0$ defines the uniform entropy upstream from the gain region. The upstream flow (positive $z$) is decelerated by the gravity step function and its entropy increases due to heating. The outer edges of the source layer ($\Psi(z) > 0$) are delimited by the vertical solid lines and the region where the source terms are at their full intensity ($\Psi(z) = 1$) is delimited by the vertical dashed lines.</figcaption></figure>

The structure of the stationary flow can be obtained by solving the following equations:

$$ \frac{\partial \rho v_z}{\partial z} = 0, \qquad (6) $$

$$ \frac{\partial}{\partial z} \left( \frac{v^2}{2} + \frac{c^2}{\gamma - 1} + \Phi \right) = \frac{\mathcal{L}}{\rho v_z}, \qquad (7) $$

$$ \frac{\partial S}{\partial z} = \frac{\mathcal{L}}{P v_z}, \qquad (8) $$

where $S$ defines the dimensionless entropy

$$ S \equiv \frac{1}{\gamma - 1} \log \left[ \left( \frac{P}{P_{\text{up}}} \right) \left( \frac{\rho_{\text{up}}}{\rho} \right)^{\gamma} \right], \qquad (9) $$

with $P$ being the pressure. When crossing the gain layer, the subsonic flow is decelerated by the gravity function and its entropy increases due to the heating function which is proportional to the density (Fig. 1).

## 2.2 Linear analysis

The presence of a negative entropy gradient in the gain region is not sufficient to drive convection because the advection time-scale through the gain region is finite (Foglizzo et al. 2006). The flow is stable to convection if the kinetic energy of the accretion flow is high enough to overcome the gravitational potential energy liberated in the vertical exchange of low and high entropy material. Foglizzo et al. (2006) proposed a linear criterion to assess whether convection can develop or if it is stabilized by advection. This criterion relies on a dimensionless parameter that compares the local buoyancy time-scale ($t_{\text{buoy}}$) to the advection time-scale ($t_{\text{adv}}$) through the gain

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Growth rate of the convective instability as a function of the horizontal wavenumber for $K_G = 3$ and $M_{\text{up}} = 0.3$. The growth rates are normalized to the maximum value of the Brunt–Väisälä frequency noted $\omega_{\text{max}}$. The ratio $\chi$ of the advective and convective time-scales is indicated on each curve. The convective instability disappears for $\chi < \chi_{\text{crit}} \sim 2.4$.</figcaption></figure>

region and writes for our setup:

$$ \chi \equiv \int_{z=-50 \text{ km}}^{z=50 \text{ km}} \left| \operatorname{Im} \left( \langle \omega_{\text{BV}}^2 \rangle^{1/2} \right) \frac{dz}{\langle v_z \rangle} \right| \sim \frac{t_{\text{adv}}}{t_{\text{buoy}}}, \qquad (10) $$

where $\omega_{\text{BV}}$ corresponds to the Brunt–Väisälä frequency defined as:

$$ \omega_{\text{BV}} \equiv (\nabla \Phi)^{1/2} \left| \frac{\nabla P}{\gamma P} - \frac{\nabla \rho}{\rho} \right|^{1/2} = \left( \frac{\gamma - 1}{\gamma} \nabla \Phi \nabla S \right)^{1/2} . \qquad (11) $$

and $\langle \cdot \rangle$ to the horizontal average over the gain region. The parameter $\chi$ is computed by using the average entropy gradient $\langle \nabla S \rangle$ in equation (11). The local buoyancy time-scale can be approximated by: $t_{\text{buoy}} \sim \omega_{\text{BV}}^{-1}$, while the advection time-scale through the gain region is defined as:

$$ t_{\text{adv}} \equiv \int_{z=-50 \text{ km}}^{z=50 \text{ km}} \frac{dz}{\langle |v_z| \rangle} . \qquad (12) $$

In the rest of the paper, the time is normalized by the advection time of the unperturbed flow through the gain region. The latter corresponds to about 20 ms and is almost constant over the range of $K_H$ explored in our study.

The linear analysis of the model is detailed in Appendix B. Fig. 2 demonstrates that if the buoyancy time-scale is short enough compared to the advection time-scale ($\chi > \chi_{\text{crit}}$), the flow is unstable for a range of horizontal wavenumbers [$k_{\text{min}}$, $k_{\text{max}}$] in a similar manner as in a shocked flow (fig. 5 in Foglizzo et al. 2006). The absence of a shock in the present setup affects the instability threshold $\chi_{\text{crit}}$ and the asymptotic scaling of $k_{\text{min}}$ as shown in Fig. 3. This figure includes the stability properties of a flow with a strong shock where energy losses by the dissociation of nuclei is chosen such that the post-shock Mach number is the same as in our model. The critical instability threshold is lowered from $\chi_{\text{crit}} \sim 3.5$ with a shock to $\chi_{\text{crit}} \sim 2.4$ without a shock. The effect of a shock on the stability properties is studied analytically in the low mach number limit in Appendix B, which confirms that the instability threshold is set by the parameter $\chi \propto (K_H K_G)^{1/2}/M$ (equation B14) and demonstrates that the asymptotic scaling of $k_{\text{max}} \propto \chi/H$ is not affected by the shock.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Range of horizontal wavenumbers allowing for the convective instability with and without a shock at the upper boundary. The full lines correspond to an upper Mach number $M = 0.3$ with $K_G = 3$. The dashed lines correspond to analytical calculations in the asymptotic limit $K_G \ll 1$, $K_H \ll 1$, and $M \ll 1$ without and with a shock (equations B44 and B61, respectively). The absence of the shock favours the instability at smaller wavenumber and a lower instability threshold $\chi_{crit} \ge 2$.</figcaption></figure>

In the limit $M \ll 1$, the critical instability threshold is lowered from $\chi_{crit} = 3.3$ with a shock to $\chi_{crit} = 2$ without a shock. Altogether, the convergence study illustrated by Fig. B1 suggests that the instability threshold is in the range 2–2.4 without a shock, and in the range 3.3–3.5 with a shock. For comparison, a broader range 2.8–3.5 of critical thresholds was obtained with the combined action of heating, cooling, and dissociation below a shock (fig. 6 in Foglizzo et al. 2006).

## 2.3 Initial perturbation

The multidimensional dynamics is triggered by a density perturbation at pressure equilibrium which is added to the stationary flow upstream from the source layer (Fig. 4). The perturbation is initially located between $z = 75$ and 425 km, thus entirely contained in the upstream region. Its vertical extent is set so that the accretion of the whole perturbation into the gain region corresponds to roughly one advection time of the stationary flow (equation 12). Note that no perturbations are present at the upper boundary condition ($z = 450$ km). In 3D, the perturbation is almost uniform in the transverse direction in order to perform a detailed comparison between the early non-linear regimes in 2D and 3D. A random noise of amplitude 0.1 per cent in density is included in each numerical cell to enable the growth of the instability in the third direction.

The perturbation models an entropy wave produced by the advective–acoustic cycle of SASI (Fig. 4). We aim at studying the coupling between the hydrodynamical instabilities considering a transient excitation of convection by SASI. The absence of a shock wave prevents a possible feedback from the gain region to generate new entropy waves which would continuously feed the instability.

This perturbation could also represent pre-collapse asymmetries originating from combustion inhomogeneities in multidimensional progenitors (e.g. Couch & Ott 2015; Müller et al. 2016) or numeri-

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Structure of the perturbation shown in entropy contrast, where the horizontal average value is subtracted from the flow. Outside of the two black solid lines, gravity, and heating are turned off ($\Psi(z) = 0$), while they are at their full intensity ($\Psi(z) = 1$) inside the black dashed lines.</figcaption></figure>

cal artefacts such as embedded Cartesian grids refined towards the centre of the domain (e.g. Ott et al. 2013).

Our study focuses on entropy perturbations and ignores vorticity perturbations inherent to the advective–acoustic cycle of SASI (Foglizzo et al. 2007) because those would be unstable to the Kelvin–Helmholtz instability in the uniform flow before reaching the gain layer. Entropy perturbations are left unchanged upstream from the gain layer and are simply advected until they possibly become unstable due to the action of buoyancy. To restrict the number of parameters, we set the horizontal wavenumber of the perturbation to $m = 5$. This corresponds to the most unstable mode for the size of the gain region considered in our study.

## 2.4 Numerical simulations

The RAMSES code (Teyssier 2002; Fromang, Hennebelle & Teyssier 2006) is employed to simulate the dynamics on a Cartesian grid. It is a second-order finite-volume code which uses the MUSCL–Hancock scheme. We employ the HLLD Riemann solver (Miyoshi & Kusano 2005) and the monotonized central slope limiter. Periodic boundary conditions are employed at the lateral edges of the domain. A constant inflow determined by the stationary flow is used at the top boundary condition. The bottom one consists of a constant outflow where ghost cells are filled with values from the stationary flow. We find that this choice minimizes the reflections at the bottom of
the box and only marginally perturbs the constant inflow from the upper boundary. The default numerical resolution is such that ($N_x \times N_z$) = (384 x 1152) in 2D and ($N_x \times N_y \times N_z$) = (384 x 384 x 1152) in 3D. The size of the numerical cells is the same in all directions. The impact of the numerical resolution will be addressed in Section 6. The gain region is discretized in 128 vertical cells in the z-direction and this represents about 60 cells per pressure scale height.

Two main parameters are varied to explore the different regimes of the convective instability and the impact of the dimensionality. The parameter $\chi_0$ corresponds to the initial value of $\chi$ (Fig. 3). It is varied from 0 to 5 to study the regime of linear instability ($\chi > \chi_{crit}$) as well as non-linearly triggered convection ($\chi < \chi_{crit}$). We also vary the amplitude of the perturbation $\delta\rho/\rho$ such that 0.1 per cent $\le \delta\rho/\rho \le 30$ per cent. In Section 3, we restrict ourselves to 2D simulations to cover the parameter space ($\chi_0$, $\delta\rho/\rho$) and assess the robustness of linear and non-linear instability criteria. A detailed analysis of the differences between 2D and 3D simulations is presented in Section 4 based on simulations which illustrate the different instability regimes identified in the next section.

# 3 NUMERICAL EXPLORATION OF THE DIFFERENT INSTABILITY REGIMES

## 3.1 Convective instability triggered non-linearly

A buoyant bubble of density $\rho$ is expected to rise against a surrounding flow of density $\rho_0$, regardless of the value of $\chi$, if the density contrast is such that

$$ \delta_{\min} \equiv \frac{|\rho_0 - \rho|}{\rho_0} \approx \frac{\langle |v_z| \rangle}{\langle g \rangle t_{\text{adv}}} \sim \mathcal{O}(1 \text{ per cent}), \quad (13) $$

as proposed by Scheck et al. (2008). This can be understood as a competition between the buoyant acceleration experienced by the bubble and the advection velocity. Fernández et al. (2014) proposed an alternate criterion based on the balance between the buoyant force ($\sim V\delta\rho g$) and the drag force ($\sim 1/2C_D S\rho_0 v^2$) exerted on a bubble:

$$ \delta_{\min} \equiv \frac{|\rho_0 - \rho|}{\rho_0} \approx \frac{C_D \langle |v_z| \rangle^2}{2l_0 \langle g \rangle}, \quad (14) $$

where $C_D$ is the drag coefficient (0.5 for a sphere) and $l_0$ the ratio between the volume V and the cross-section S of the bubble.

To test these criteria, we perform a set of simulations using a linearly stabilized flow with $\chi_0 = 1.5$ and different perturbation amplitudes. The criteria (13) and (14), respectively, give thresholds of 0.5 per cent and 0.8 per cent, which are in good agreement with our numerical simulations (Fig. 5, top panel). However, our results show that perturbation amplitudes slightly above the thresholds (13) and (14) are not strong enough to trigger convective overturns but only temporary buoyant motions. For such perturbation amplitudes, the maximum vertical velocity in the gain region becomes negative after a few advection times showing that ascending motions are completely suppressed. Our simulations indicate that the criteria used to estimate the rise of a buoyant bubble are not sufficient conditions to trigger turbulent convection.

## 3.2 A self-sustained instability?

We now focus on the dynamics in the gain region at later times. In the set of simulations described in the previous section, we observe

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time evolution of the maximum upward velocity in 2D simulations. The velocities are normalized by the upstream velocity. Top panel: set of simulations of linearly stable flows with $\chi_0 = 1.5$ and different perturbation amplitudes. In all cases, the velocity reaches negative values showing that buoyant motions are suppressed. Bottom panel: set of simulations with $\delta\rho/\rho = 30$ per cent and different values of $\chi_0$. Negative values are close to be reached in all cases where $\chi < \chi_{crit} \sim 2.4$ showing that the instability is not self-sustained when triggered in a linearly stable flow.</figcaption></figure>

that even for the largest perturbation amplitude considered, the instability is completely damped after 15 advection time-scales (about 300 ms), since ascending motions no longer exist (Fig. 5, top panel). In each simulation, the accretion of the initial perturbation into the gain layer lasted less than the first advection time-scale. This set of simulations shows that if a large perturbation triggers convection in a linearly stabilized flow, the instability does not reach a permanent turbulent regime and convection is eventually suppressed. The damping time-scale increases with stronger perturbations but the flow always adjusts itself to a linearly stable flow (sub-critical flow) with a value of $\chi$ below the instability threshold.

The damping time-scale also depends somehow on the distance to the linear instability threshold (Fig. 5, bottom panel). In a second set of simulations, performed with $\delta\rho/\rho = 30$ per cent and several values of $\chi_0$, we observe that only the cases with an initial value of $\chi$ above the instability threshold are able to maintain turbulent convective motions over at least 15 advection time-scales. Our results show that a single excitation is not sufficient to trigger a self-sustained instability when the flow is linearly stable. This can only occur in situations where the dynamics is continuously fed by non-linear perturbations such as in a simulation performed with a Cartesian grid which would produce significant noise (Ott et al. 2013) or when SASI reaches large amplitudes (Cardall & Budiardja
2015; Summa et al. 2016). Moreover, our results show that the perturbative analysis of Foglizzo et al. (2006) holds in cases where the perturbed flow is already in the non-linear regime. In more realistic 3D simulations, the $\chi$ parameter seems to be a reliable criterion to assess whether convection ($\chi \gtrsim 3$) or SASI ($\chi \lesssim 3$) dominates the dynamics shortly after the spherical symmetry is broken (Hanke et al. 2013; Takiwaki et al. 2014; Couch & O'Connor 2014; Abdikamalov et al. 2015). The asymptotic value of $\chi$ in the turbulent convective regime will be discussed in Sections 5 and 6 where we investigate the impact of dimensionality and resolution.

The exploration of the parameter space ($\chi_0$, $\delta\rho/\rho$) enables us to distinguish three different regimes:

(i) linear instability if $\chi_0 > \chi_{crit}$,

(ii) linear stability if $\chi_0 < \chi_{crit}$ and $\delta\rho/\rho \lesssim 1$ per cent,

(iii) transient convection if $\chi_0 < \chi_{crit}$ and $\delta\rho/\rho \gtrsim 1$ per cent.

In the latter case, convection is only transitory and cannot be self-sustained. The instability damping time-scale increases with higher perturbation amplitudes and depends in a less clear way on the nearness to the instability threshold (Fig. 5). Nevertheless, the interaction between ascending bubbles and a shock wave might be strong enough to trigger convection in a linearly stable flow due to a non-linear cycle different from SASI. It cannot be excluded that the inclusion of a shock wave allows new ways to sustain convection, at least for values of $\chi$ slightly below the linear instability condition.

# 4 IMPACT OF THE DIMENSIONALITY

## 4.1 Early non-linear phase

The first differences between 2D and 3D dynamics arise when perturbations reach non-linear amplitudes. If the instability is triggered by a low-amplitude perturbation in a flow where $\chi_0 > \chi_{crit}$, buoyant bubbles rise faster in 3D than in 2D (Fig. 6, top panels). In a simulation run with the parameters $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent, the highest 3D bubbles rise almost twice as fast against the flow and the instability saturates at a much larger amplitude (Fig. 7). This is deduced from the time evolution of the location of the highest buoyant bubble. Its position is obtained by tracking the altitude of the uppermost negative radial entropy gradient which is due to the entropy contrast between the background flow (lower entropy) and the highest buoyant bubble (higher entropy). This transition is much sharper than the one related to the background flow.

A crude estimate of the velocity of the bubbles can be made by computing the terminal velocity $v_{ter}$ which corresponds to the balance between the drag force ($S\rho_0 v_{ter}^2$) and buoyancy ($V(\rho_0 - \rho)g$):

$$v_{ter}^2 = \frac{(\rho_0 - \rho) V g}{\rho_0 S} \qquad (15)$$

The terminal velocity increases with higher volume-to-surface ratio ($V/S$). Couch (2013) interpreted earlier and easier explosions in their 2D simulations as a consequence of larger $V/S$ ratios compared to 3D. This explanation is based on different bubble shapes between 2D and 3D. Large-scale axisymmetric bubbles are considered to be more efficient to overcome ram pressure than 3D bubbles that fragment to small scales. However, the ratio $V/S$ depends mostly on the vertical extension of the bubbles. For a bubble moving in the radial direction, the relevant surface is in the angular directions. The ratio $V/S$ is therefore approximately the vertical extension of the bubbles, within a factor of order unity depending on the exact shape. For example, comparing a spherical bubble of a given radius with a

toroidal one of the same radius, one finds that the terminal velocity is only slightly higher in 2D, by a factor of order unity $\sqrt{3\pi}/8$. Interpreting the role of dimensionality on CCSN simulations as a consequence of specific bubble properties is in conflict with our model where bubbles are clearly faster and reach higher altitudes in 3D than in 2D.

The acoustic feedback produced by the convective motions is able to alter the subsonic upstream flow. As a consequence, the ram pressure against the buoyant bubbles decreases with time (Fig. 7). The upstream acoustic feedback is stronger in 3D, reducing the ram pressure more than in 2D simulations. However, we observe that when the ram pressure deviates by more than 1 per cent between 2D and 3D, the buoyant bubbles have already reached a larger altitude in 3D. This shows that the alteration of the upstream flow by convection motions cannot explain the early discrepancies in terms of altitudes reached by buoyant bubbles. It is conceivable that the acoustic feedback plays a role at a later stage to set the saturation amplitude of the instability. In a model with a shock wave and a supersonic upstream flow, a strong acoustic feedback would increase the post-shock pressure and induce a shock expansion. Our study suggests that this effect could be larger in 3D than in 2D.

Besides, our results seem consistent with conclusions drawn from studies focusing on the role of dimensionality in the Rayleigh-Taylor instability. The mixing zone is found to broaden faster in 3D as plume-like structures penetrate much deeper than 2D planar structures in simulations of incompressible flows without advection (Young et al. 2001; Anuchina et al. 2004). Similar conclusions were obtained in the context of the Rayleigh-Taylor mixing in the stellar envelope during the explosion of massive stars (Kane et al. 2000; Hammer, Janka & Müller 2010). The growth of the Rayleigh-Taylor fingers is faster and the velocities of clumps of heavy elements are higher in 3D. Considering a more elaborate model of the competition between drag force and buoyancy, Hammer et al. (2010) showed that the temporal evolution of the $V/S$ ratio favours a faster rise of the bubbles in 3D. Our model underlines that during the early non-linear phase of the instability, the rise of buoyant bubbles and the widening of the mixing zone are faster in 3D, even when advection is considered.

## 4.2 Heating efficiency

Fig. 8 shows the time evolution of the spatial distribution of entropy in 2D and 3D. Each vertical entropy profile is defined such that a given fraction of the cells have a lower entropy value. In the early non-linear phase, the convective instability exhibits a wider mixing region in 3D, while mixing is more localized in 2D. As a consequence, the vertical profiles related to high entropy values show more pronounced peaks in 2D. The rise of the 3D bubbles stops when they reach the upper edge of gain region above which gravity is switched off (Fig. 7, top panel). After several advection times, the vertical profiles of entropy become higher in the whole gain layer in 3D (Fig. 8, bottom panel).

The generation of entropy in the gain region appears to be higher in 3D than in 2D (Fig. 9, top panel). At the end of the simulations performed with $\chi_0 = 5$, the entropy has increased about four times more in 3D than in 2D. We note that the amplitude of the initial perturbation does not affect the asymptotic dynamics because final entropy values are independent of the perturbation both in 2D and 3D. In the regime of transient convection, the entropy increase is related to the amplitude of the perturbation and is higher in 3D (Fig. 9, top panel, purple curves). The damping time-scale of the instability is much longer in 3D suggesting that perturbations
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Snapshots of entropy for models with $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent after 2.91 (top) and 5.33 (bottom) advection time-scales in 2D (left) and 3D in the vertical slice $y = 0$ (right). The horizontal black lines are defined in Fig. 4. An animated version of this figure is available in the online journal.</figcaption></figure>

may have a more profound impact on convection-dominated cases simulated in 3D than in 2D. A direct comparison of the heating rates, proportional to the mass enclosed in the gain region, does not explain the discrepancies in terms of entropy. It turns out that the rates are slightly higher in 2D throughout the simulations (Fig. 9, middle panel). The strong impact of dimensionality on the entropy variations is suggestive of a heating process that differs significantly between 2D and 3D. We defer the investigation on the origin of this additional heating to Section 5 where the turbulence induced by the instability is analysed in detail.

The morphology of the flow is strongly impacted by the dimensionality. In 2D, the flow is stirred by large vortices, whereas 3D models undergo stronger turbulent mixing to small scales. Phase mixing is essential to disrupt large-scale entropy perturbations even in the absence of heating (case $\chi_0 = 0$ of Fig. 5, bottom panel). The dynamics generates large vortices in 2D whose vertical extension can reach the height of the gain layer. This can be seen as a consequence of a greater conversion of gravitational potential energy to kinetic energy by the Rayleigh-Taylor instability (Young et al.

2001; Cabot 2006). In 2D, the interface between ascending motions and downflows is unperturbed. The downflows are able to sustain throughout the gain layer and channel lower entropy material located in the upper part of the gain layer to the downstream region where heating is absent (Fig. 6, bottom left). On the contrary, 3D downflows are disrupted by turbulent mixing before leaving the gain layer. Such a mixing is known to be more efficient in 3D in the case of the Rayleigh-Taylor instability (Cabot 2006; Hammer et al. 2010). As the flow leaves the gain layer ($z = -50$ km), we observe that in 3D a smaller fraction of the numerical cells have an entropy lower than the stationary flow compared to 2D (Fig. 8, bottom panel). This implies that downflows are more efficient to channel cold entropy material outside of the gain region in 2D, while their fragmentation to small scales in 3D keeps buoyant material for longer times and enhances heating.

The combine action of more efficient downflow braking and stronger upstream acoustic feedback in 3D produces smaller average downflow velocities (Fig. 10). The advection time is thus slightly larger in 3D and this should in principle increase the heat-
<figure><img src="image_9.png" /></figure>

**Figure 7.** Altitude of the upper edge of the uppermost buoyant bubble (top panel) and ram pressure, defined as $ρv_2^2$, (bottom panel) as a function of time for a model with $χ_0 = 5$ and $δρ/ρ = 0.1$ per cent simulated in 3D (thick curves) and 2D (thin dashed curves). The horizontal black lines are defined in Fig. 4. The ram pressure is averaged between the planes $z = 30$ and 45 km. The vertical blue dashed line marks the time at which the ram pressure starts to deviate by more than 1 per cent between 2D and 3D. At this time, buoyant bubbles have already reached higher altitudes in 3D and this initial gap is not related to differences in terms ram pressure against the bubbles. The later evolution is characterized by a lower ram pressure in 3D. This results from a more important acoustic feedback produced by convection which alters the subsonic upstream flow more strongly.

ing efficiency compared to 2D. However, this is contradicted by the previous comparison of the heating rates.

## 4.3 Asymptotic regime

The evolution of the $χ$ parameter is related to the way the flow adjusts itself compared to the marginal stability limit (Fig. 9, bottom panel). The $χ$ parameter is computed from horizontally averaged flows (equation 10) as proposed by Fernández et al. (2014).

When convection is triggered non-linearly by a strong perturbation ($χ_0 = 1.5$, $δρ/ρ = 30$ per cent), the parameter initially takes large values during the advection of the perturbation. This is due to the strong entropy gradients contained in the perturbation (Fig. 4) which induce a larger value of $χ$ for a brief period of time. Then, the value of $χ$ drops much faster in 3D than in 2D. About ~0.1 $t_{adv}$ after the perturbation has started to accrete through the gain layer, the 3D flow already reaches sub-critical values ($χ \lesssim 2.4$) because the perturbation is efficiently disrupted by strong turbulent mixing. It takes more than two advection times to reach the same level in 2D. In both cases, the flow eventually returns to a sub-critical state. Besides, this comparison confirms that even in 3D a strong perturbation cannot lead to self-sustained convection when the flow is linearly stabilized (see Section 3.2).

The discrepancies are more remarkable in the linear instability regime (Fig. 9, bottom panel, black and orange curves). In 2D, the flow maintains a $χ$ value well above the stability threshold ($χ \gtrsim 3$) whereas in 3D, the flow adjusts itself to a slightly sub-critical flow ($χ \approx 2.2$). Slightly lower average vertical velocities in 3D (Fig. 10, top panel) should in principle result in a higher value

<figure><img src="image_10.png" /></figure>

**Figure 8.** Vertical profiles of entropy obtained from 3D (thick curves) and 2D (thin dotted curves) simulations of a model with $χ_0 = 5$ and $δρ/ρ = 0.1$ per cent, after 2.91 (top panel) and 4.5 (bottom panel) advection time-scales. Each profile shows the entropy such that a given fraction of the numerical cells at a given altitude $z$ have a lower value. The entropy of the stationary flow is shown in black.

of $χ$ than in 2D. However, turbulent mixing seems to be part of the explanation of the gap between 2D and 3D values of $χ$. After the instability fully develops, a smaller fraction of the gain layer remains unstable to convection in 3D (Fig. 10, middle panel). In 2D, the flow is dominated by large vortices which create a less stable configuration. The contribution of the lower part of the gain layer ($z \lesssim 0$ km) seems responsible for the greater values of $χ$ in 2D (Fig. 10, bottom panel). These results suggest that any analytical description of turbulence induced by convection in CCSNe (e.g. Murphy & Meakin 2011 in 2D) should reflect that difference when applied to our model.

# 5 THE ROLE OF TURBULENCE

## 5.1 Kinetic energy

Convection may facilitate the shock revival in several ways compared to 1D cases. Non-radial motions induced by instabilities in multidimensional simulations increase the advection time-scale and bring it closer to the time needed for a sufficient energy deposition in the gain layer (Murphy & Burrows 2008). Turbulent pressure generated by convective motions represents an additional pressure support which pushes the shock wave outward and lowers the critical luminosity required to achieve shock revival (Burrows et al. 1995; Murphy et al. 2013; Couch & Ott 2015; Müller & Janka 2015). Both effects may combine to enhance neutrino heating and favour multidimensional explosions. Comparing 2D and 3D simulations,
<figure><img src="image_11.png" /><figcaption>Figure 9. Top panel: time evolution of the average entropy variations in the gain layer in 2D (dashed curves) and 3D (solid curves) for three sets of parameters. The variations are computed as: $\Delta S = \langle S \rangle - \langle S \rangle_0$ where $\langle S \rangle$ is the average entropy at a given time and $\langle S \rangle_0$ at the initial time. Middle panel: time evolution of the ratio of heating rates between 3D and 2D simulations. In all cases, the heating rate is slightly higher in 2D than in 3D. Lower panel: time evolution of the $\chi$ parameter in 2D (dashed curves) and 3D (thick curves). In the regime of transient convection, the flow returns below the critical value $\chi \sim 2.4$. In the linear instability regime, only the 3D dynamics is able to bring the flow to a sub-critical state.</figcaption></figure>

Couch & Ott (2015) concluded that the strength of turbulence is overpredicted in axisymmetric simulations.

The turbulent kinetic energy in the gain layer can be defined as:

$$E_{\text{turb}} = \frac{1}{2} \rho [(v_z - \langle v_z \rangle)^2 + v_x^2 + v_y^2], \quad (16)$$

where $\langle v_z \rangle$ corresponds to the horizontally averaged vertical velocity. In our model, we find that the turbulent kinetic energy is larger in 2D by about 15 per cent during most of the simulation (Fig. 11). Besides, the turbulent kinetic energy is roughly in equipartition in 2D between the two directions. This can be seen as a consequence of the domination of large-scale vortices. In 3D, we find that the horizontal contributions are in equipartition ($E_x \approx E_y$) and their sum amounts to the vertical contribution ($E_z \approx E_x + E_y$). This is the case in 3D models where buoyant motions are induced along a preferred direction (Murphy et al. 2013). Nevertheless, higher Reynolds stresses observed in 2D do not seem to support larger entropy values in our 3D simulations (Fig. 9).

The impact of dimensionality can also be investigated by comparing the turbulent energy cascades in 2D and 3D. To connect our study with previous ones (Hanke et al. 2012; Dolence et al. 2013; Couch &

<figure><img src="image_12.png" /><figcaption>Figure 10. Top panel: vertical profiles of vertical velocity in the gain region of simulations performed in 2D (blue) and 3D (green) with the parameters $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent. The quantities are averaged horizontally and over a period of four advection time-scales. The vertical profiles are normalized by the upstream velocity of the stationary flow. Middle panel: vertical profiles of the Brunt–Väisälä frequency. The profiles are normalized to the maximum value of the Brunt–Väisälä frequency in the gain region, noted $\omega_{\text{max}}$. The upper part of the gain region is unstable to convection in 3D, while the whole layer remains unstable in 2D. Bottom panel: vertical profiles of the cumulative integral of $\partial\chi/\partial z$ computed from the upper edge of the gain region, noted $z, \text{up}$. The higher $\chi$ values obtained in 2D seem to result from the contribution of the lower part of the gain layer.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Turbulent kinetic energy integrated over the gain region as a function of time for simulations performed in 2D (dashed) and 3D (solid) with the parameters $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent. The quantities are normalized by the kinetic energy of the stationary flow upstream from the gain layer, noted $E_{\text{kin, up}}$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Spectra of turbulent horizontal kinetic energy density as a function of the wavenumber modulus *k*. The spectra are computed from 2D (blue) and 3D (green) simulations performed with the parameters $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent. The vertical black line represents the injection scale *k* = 5 related to the initial perturbation in 2D and 3D. In 2D, the reverse turbulent cascade follows a power-law $E(k)\propto k^{-3}$ (purple dotted-dashed line), while the decay law in 3D follows $E(k)\propto k^{-1}$ at large spatial scales and $E(k)\propto k^{-5/3}$ at intermediate ones.</figcaption></figure>

O’Connor 2014; Handy, Plewa & Odrzywołek 2014; Couch & Ott 2015; Abdikamalov et al. 2015; Radice, Couch & Ott 2015; Radice et al. 2016), we compute the turbulent kinetic energy density spectra considering only horizontal motions. The kinetic energy density is decomposed into Fourier coefficients:

$$ \hat{E}_{\perp}(k_x, k_y) = \int_{V_g} \exp\left(-2\pi i \frac{k_x x + k_y y}{L}\right) \sqrt{\rho(x, y)} v_{\perp}(x, y) dV, \quad (17) $$

where *L* represents the horizontal extent of the gain region (*L* = 300 km in our simulations), $V_g$ the volume of the gain region in which the source terms are at their full intensity ($\Psi(z) = 1$), and $v_{\perp}$ the horizontal velocity component. The Fourier coefficients are computed at each discrete vertical coordinate of the grid and averaged over the volume $V_g$. Note that $k_y = 0$ in 2D. The total horizontal kinetic energy density is then

$$ E(k) = \sum_{k-1 < \|\|(k_x, k_y)\| \le k} |\hat{E}_{\perp}(k_x, k_y)|^2. \quad (18) $$

The Fourier spectra are normalized by $\sum_k E(k)$ to obtain unity integrals.

In 2D, the reverse turbulent cascade transfers kinetic energy to the largest spatial scales and this could explain why explosions are easier (Hanke et al. 2012; Couch & O’Connor 2014). In the inertial range, entrophy is transferred from large scales to the dissipation regime following $E(k)\propto k^{-3}$ (see Kraichnan 1967; Hanke et al. 2012 for details). Such a trend seems to hold in our 2D simulations for $5 \lesssim k \lesssim 30$ (Fig. 12), even though the identification of the inertial range is not obvious. In 3D, turbulence proceeds differently because a forward cascade transfers kinetic energy from large scales to the dissipation regime, following a decay power-law $E(k)\propto k^{-5/3}$ (Landau & Lifshitz 1959). Our simulations seem to show such a power-law decay (Fig. 12). There is no consensus on the exact shape of the power law in CCSN simulations. Some studies re-

<figure><img src="image_15.png" /><figcaption>**Figure 13.** 3D visualization of the entropy isosurfaces obtained from a simulation with $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent after 4.92 advection time-scales. Large and intermediate spatial scale structures can clearly be identified in the upper part of the gain layer. Animated versions of this figure are available in the online journal.</figcaption></figure>

ported a decay law close to $E(k)\propto k^{-5/3}$ (Hanke et al. 2012; Handy et al. 2014; Radice et al. 2016), while others observed a shallower decay $E(k)\propto k^{-1}$ (Dolence et al. 2013; Couch & O’Connor 2014; Abdikamalov et al. 2015). Differences may result from a misidentification of the inertial range. A transition from a $k^{-1}$ scaling for $k \lesssim 15$ to a $k^{-5/3}$ scaling for $10 \lesssim k \lesssim 40$ appears in our simulations. The nature of the power law in the inertial range reflects the way dissipation proceeds. A limited numerical resolution could affect the shape of the power-law decay and only the finest resolutions, unachievable for current state-of-the-art CCSN simulations, could lead to a behaviour predicted by Kolmogorov’s theory (Radice et al. 2015, 2016). It can also be argued that the assumptions made in the classical theory do not apply to post-shock turbulence in CCSNe (Abdikamalov et al. 2015).

Although the largest spatial scale structures are more pronounced in 2D, the spatial scales just below, between $k = 3$ and 10, seem almost equally favoured in our 2D and 3D simulations. This is consistent with the large entropy structures seen in the upper part of the gain layer (Fig. 6, bottom right). The presence of large-scale structures can also be observed in Fig. 13 which provides a 3D visualization of entropy isosurfaces of the simulation performed with $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent, when the instability is fully developed. Nevertheless, the discrepancies in terms of turbulent kinetic energy do not seem to be directly connected with a stronger heating that occurs in our 3D simulations.

## 5.2 Turbulent dissipation

In this section, we use a different approach to investigate the impact of dimensionality on turbulence. Following the work of Murphy & Meakin (2011) and Murphy et al. (2013), we employ a mean-field decomposition of the entropy equation into background and turbulent flows to identify the driving agent of the discrepancies in terms of heating. To this end, we decompose the variables into mean and fluctuation using two methods depending on the quantities. The Reynolds decomposition is defined according to $\phi = \bar{\phi} + \phi'$ where the mean-field average of the fluctuation $\phi'$ is $\bar{\phi}' = 0$. The Favrian decomposition corresponds to $\psi = \bar{\psi} + \psi''$. Here, the Favrian-averaged quantity is a density-weighted average:
$\tilde{\psi} = \overline{\rho\psi}/\overline{\rho}$. The Frevian fluctuation is such that $\tilde{\psi}'' = 0$. This second decomposition appears often to be more convenient for terms containing density in equations. We refer the reader to Mocák et al. (2014) for a complete overview of the derivation of hydrodynamics equations using the Reynolds-averaged Navier-Stokes method.

The entropy equation is written as:

$$ \frac{\partial (\rho S)}{\partial t} = -\nabla \cdot (\rho S v_z) + \frac{\rho \dot{q}}{T} + \frac{\rho \epsilon}{T}, \quad (19) $$

where $\dot{q}$ is the local heating rate, $T$ the temperature, and $\epsilon$ the heat due to turbulent dissipation of kinetic energy. Applying a mean-field decomposition to equation (19), we obtain:

$$ \frac{\partial \tilde{\rho} \tilde{S}}{\partial t} = -\nabla \cdot (\tilde{\rho} \tilde{v}_z \tilde{S}) - \nabla \cdot (\overline{\rho v_z'' S''}) + \left(\frac{\overline{\rho \dot{q}}}{T}\right) + \left(\frac{\overline{\rho \epsilon}}{T}\right). \quad (20) $$

The derivation is detailed in Appendix A. These terms are labelled from (i) to (v) and represent:

(i) $\partial \tilde{\rho} \tilde{S} / \partial t$, the generation rate of entropy,

(ii) $\nabla \cdot (\tilde{\rho} \tilde{v}_z \tilde{S})$, the entropy flux due to the mean flow across the gain layer,

(iii) $\nabla \cdot (\overline{\rho v_z'' S''})$, the entropy flux due to the turbulent flow across the gain layer,

(iv) $(\overline{\rho \dot{q}} / T)$, the heating rate due to neutrino absorption,

(v) $(\overline{\rho \epsilon} / T)$, the turbulent dissipation.

This last term is evaluated from equation (20) and can also be approximated by $W_b/T$ (Murphy et al. 2013), where $W_b$ is the buoyancy work, defined as:

$$ W_b = \int \overline{\rho' v_z' g} dV. \quad (21) $$

The term $W_b/T$ is labelled as (vi) in the following. It approximates the turbulent dissipation if all the energy injected into turbulence by the buoyancy force is dissipated into heat. This would be expected in a quasi-steady-state situation if the turbulent energy advected out of the gain region can be neglected.

The terms (i) to (vi) are averaged horizontally and over four advection times during the fully non-linear regime. These terms are then integrated over the height of the gain region and compared to pinpoint which ones play a significant role in the entropy production (Fig. 14). Note that all the terms are positive except the entropy flux due to the turbulent flow in 2D. The neutrino heating rate appears to be the main contributor to the average entropy. It is almost balanced by the entropy flux due to the mean flow which accounts negatively to the entropy production (equation 20). These terms are almost independent of the dimensionality and in particular $\dot{q}$ is only slightly higher in 2D, as already shown in Fig. 9.

The second main heating source is the turbulent dissipation of kinetic energy [term (v)]. In 3D, it accounts for almost half of the heating due to neutrinos. The entropy generation rate is around 20 times higher in 3D than in 2D [term (i)], so is the ratio of turbulent dissipation. This suggests that the additional entropy production in our 3D simulations results from the contribution of the turbulent dissipation. In 3D, the term related to the turbulent flux is slightly lower than the turbulent dissipation but accounts negatively to the average entropy. In 2D, this quantity only marginally contributes to increase the average entropy since it is the smallest term.

The buoyancy work is almost balanced by turbulent dissipation in 3D as shown by (Murphy et al. 2013). This is not satisfied in 2D, where the buoyancy work is an order of magnitude larger than the turbulent dissipation. The much lower dissipation in 2D may be interpreted as a consequence of the large-scale vortices being

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Time-averaged values of entropy production terms in a case with $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent. The definitions of the six terms are given in the text. Top panel: the values obtained in the 3D run (symbols +) are compared to the ones from the 2D run (symbols x). All terms are positive except (iii) in 2D. Bottom panel: ratios of 3D to 2D average values for the six terms. The horizontal dashed line denotes the ratio of one.</figcaption></figure>

expelled from the gain layer (Fig. 6), whereas those are dissipated in 3D. The latter is more efficient to dissipate kinetic energy because of the forward energy cascade to small scales. Turbulent dissipation is less likely in 2D, since large vortices and downflows are not efficiently disrupted.

Using various approaches and setups, it was shown that turbulent motions generated by hydrodynamical instabilities in multidimensional simulations can account for a reduction of the critical neutrino luminosity of 20–30 per cent compared to 1D (Murphy & Burrows 2008; Hanke et al. 2012; Müller & Janka 2015; Fernández 2015). The hierarchy between 2D and 3D is less clear and requires analysis of the properties of turbulence in the gain region in order to pinpoint which ingredients play a crucial role. Applying the turbulence model developed by Murphy & Meakin (2011) and Murphy et al. (2013) on the explodability condition (Burrows & Goshy 1993; Murphy & Dolence 2017), Mabanta & Murphy (2018) conjectured that turbulent dissipation of kinetic energy is the dominant effect in the reduction of the critical luminosity in multidimensional simulations compared to 1D. Our results tend to support their hypothesis and the higher rate of entropy production in 3D seems to be connected to a more efficient dissipation of kinetic energy into heat. Our approach shows that such a physical phenomenon is less effective in 2D due to the particular tendency of favouring large spatial scales.

## 6 DEPENDENCE ON NUMERICAL RESOLUTION

We explore the impact of the numerical resolution on the 2D and 3D dynamics. The default resolution used in the previous sections is such that the uniform grid contains $N_{z,g} = 128$ vertical cells in the gain layer. The numerical resolutions considered in this section are $N_{z,g} = \{32, 64, 128, 256, 512\}$ in 2D and $N_{z,g} = \{32, 64, 128\}$ in 3D. The resolution is changed accordingly in the whole computational domain and in all directions to maintain the numerical cell size identical in each direction. For reference, our highest 3D resolution corresponds to twice the highest resolution considered in the studies of Abdikamalov et al. (2015) and Roberts et al. (2016) and to the
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Results obtained from 2D (thin curves) and 3D (thick curves) simulations performed with the parameters $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent for different resolutions labelled according to the number of vertical cells in the gain layer. Top panel: time evolution of the average entropy in the gain layer in 2D and 3D simulations. Lower panel: time evolution of the $\chi$ parameter in 2D and 3D simulations. The flow adjusts itself below the marginal stability (horizontal dashed line) in 3D for high enough resolutions.</figcaption></figure>

resolution labelled as '2x' in the work of Radice et al. (2016). Figs 15 and 16, respectively, give an overview of the influence of resolution on the dynamics in the linear instability regime and in the regime of transient convection. The main difference between the two regimes is that numerical resolution has a greater impact when the instability is triggered by a perturbation of a large amplitude.

Regarding the linear instability regime, the higher entropy production rate is not dramatically dependent on the resolution either in 2D or in 3D (Fig. 15, top panel). On the one hand, if the number of vertical cells in the gain layer is too low, e.g. $N_{z,g} = 32$, the instability is triggered slightly earlier. This is probably due to the relaxation of the initial flow on the grid which deviates from the stationary solution. Besides an earlier onset of the instability at this low resolution, the entropy values reached may be artificially higher than for other resolutions. This suggests that the dynamics cannot be properly simulated in our model with such a resolution. On the other hand, if $N_{z,g} \ge 64$, the linear growth of the instability and the entropy production are very similar in all cases. This shows that the discrepancies between 2D and 3D concerning the entropy vari-

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Same as in Fig. 15, but for the parameters $\chi_0 = 1.5$ and $\delta\rho/\rho = 30$ per cent.</figcaption></figure>

ations discussed in Section 4.2 are not resolution-dependent. The lowest 3D resolution case produces a value of $\chi$ which is similar to 2D simulations (Fig. 15, bottom panel) and the flow does not adjust itself to a slightly super-critical state. Other resolutions are not shown in 2D because their influence on $\chi$ is very minor. In 3D, a resolution such that $N_{z,g} \ge 64$ seems sufficient to bring the flow to a sub-critical state.

In Section 4.2, we pointed out that 3D favours a larger entropy production when convection is triggered by a large amplitude perturbation. The same conclusion seems to hold with a higher resolution both in 2D and 3D (Fig. 16, top panel). The maximum entropy value reached rises with increasing dimensionality and resolution. Moreover, the damping time-scale is also longer in 3D and with increasing resolution. In 2D, clear signs of convergence are witnessed between the three highest resolutions considered. It is unclear if a similar trend would be obtained in 3D. Only a too coarse resolution, e.g. $N_{z,g} = 32$, prevents the flow from reaching a sub-critical state that is expected for a simulation where $\chi_0 = 1.5$ (Fig. 16, bottom panel).

We now apply the same mean-field decomposition of the entropy equation (20) on 3D simulations of the linear instability regime with different numerical resolutions. The results indicate that except for turbulent dissipation at the lowest resolution, the different
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Top panel: time-averaged values of entropy production terms in a case with $\chi_0 = 5$ and $\delta\rho/\rho = 0.1$ per cent run in 3D with different resolutions, labelled according to the number of vertical cells in the gain layer. The definitions of the six terms are given in Section 5.2. All the terms are positive. Bottom panel: same as above but normalized by the values obtained in the highest resolution run ($N_{z,g} = 128$). The horizontal dashed line denotes the ratio of one.</figcaption></figure>

terms do not vary by more than 30–40 per cent (Fig. 17). This may be expected from the minor variations seen in entropy production rates with resolution (Fig. 15, top panel). Abdikamalov et al. (2015) showed that even for their highest resolution of 66 radial cells in the gain layer, the inertial range of the turbulent energy cascade is barely resolved, thus favouring greater dissipation. This could explain why turbulent dissipation decreases with finer resolution in our simulations. Nevertheless, our study suggests that simulating turbulent dissipation with a level of ~20 per cent accuracy seems achievable with a resolution of about 100 vertical cells in the gain layer (see Müller et al. 2016 for a discussion). Even if a resolution 10 times finer than ours could be needed to accurately resolve the inertial range, it appears that such an improvement does not affect global quantities such as the shock radius (Radice et al. 2016). The influence of resolution on the turbulent dissipation and the entropy variations in a convection-dominated case seems to be already captured by our set of 3D simulations with limited numerical resolutions.

## 7 DISCUSSION

Our comparative study of the 2D and 3D convective instabilities may be affected by some modelling ingredients used in this work. In the following, we discuss the differences with more realistic models and the implications of our results for self-consistent CCSN simulations. Our model consists of a static gain layer that is not delimited by a shock wave above and a stratified cooling layer

below. The absence of a shock wave lowers the instability threshold and favours the instability for long-wavelength perturbations (Fig. 3). Replacing the constant inflow upper boundary condition by a shock wave would provide a continuous source of perturbations. The acoustic retroaction produced by the dynamics in the central region would deform the stationary shock and generate new entropy-vorticity perturbations that may seed the instability. Our study indicates that the upstream acoustic feedback is stronger in 3D and this may correspond to larger shock expansions than in 2D. In Section 3, we considered an extreme case of a perturbation amplitude of 30 per cent to assess its impact on a linearly stabilized flow. Pre-shock density fluctuations related to combustion inhomogeneities may reach up to 10–15 per cent (Müller & Janka 2015). The ones related to SASI may be lower with an upper estimate of 10 per cent according to the models of Kazeroni, Guilet & Foglizzo (2016). Grid induced perturbations could produce fluctuations of 1–10 per cent (Ott et al. 2013). The lower the amplitude, the shorter the damping time-scale of the instability. According to Fig. 5, the damping time-scale may be of the order of only a few advection times for realistic perturbation amplitudes. This implies that a continuous feeding by perturbations could lead to turbulent convection even if the flow is not linearly unstable to convection. This is the case when secondary convection is seeded by large amplitude SASI oscillations (Cardall & Budiardja 2015; Summa et al. 2016). Such a situation can also occur when buoyant bubbles repeatedly push the shock or when numerical features produce noise continuously.

The absence of a cooling layer and a shock wave could impact the way convection modifies the flow. In their 2D models, Fernández et al. (2014) obtained asymptotic values of $\chi$ below the instability threshold for non-exploding, convection-dominated cases whereas this value remains well above the threshold in our 2D simulations. The dynamics seems more efficient at stabilizing the flow in 2D in their study (Fernández et al. 2014, Fig. 7) than in our case (Fig. 10). The 3D dynamics may be less impacted by the boundary conditions since the flow reaches a sub-critical state (Fig. 15). Our 2D results also differ from more realistic models regarding the balance between buoyancy work and turbulent dissipation. The equilibrium is almost reached in 3D, but the buoyancy work is one order of magnitude larger than dissipation in 2D (Fig. 14). Such a balance was reached in the 2D simulations of Murphy & Meakin (2011). A possibility is that 2D vortices are not dissipated as much as in a case where they are kept inside the gain region due to the stabilizing cooling layer. Our setup may amplify the difference between 2D and 3D in terms of the amount of kinetic energy dissipated into heat.

Similarly, our setup may artificially widen the gap between 2D and 3D in terms of entropy production rates. Hanke et al. (2012) observed only slightly higher average entropy values in the gain layer of their 3D simulations compared to their 2D ones. The difference was more pronounced in the work of Dolence et al. (2013), but of the order of 10 per cent. The latter study, using a Cartesian grid, might have been more prone to develop grid perturbations that could have artificially seeded convection. The work of Hanke et al. (2012) employed a spherical grid that generates less fluctuations. It is difficult to speculate whether a Cartesian grid would necessarily produce more heating via turbulent dissipation due to higher amplitude perturbations because the discrepancies between different groups could equally come from several specific approximations used in each study such as the neutrino transport or the microphysics. It should be noted that the average entropy value does not seem to represent a reliable indicator of the readiness of explosion (Hanke et al. 2012; Couch 2013).
Our comparative study of 2D and 3D neutrino-driven convection employs only Cartesian coordinates, neglecting the effects of geometric convergence, while previous works compared 2D axisymmetric geometry with 3D spherical or Cartesian geometries. Foglizzo et al. (2006) showed that the linear instability threshold is not affected by the geometry in their similar model including a shock wave. The same applies to the degree of the first unstable mode. However, conclusions cannot be drawn beyond the linear phase. Our study sheds light on a higher entropy creation in 3D, possibly related to more efficient turbulent dissipation. It remains unclear to which extent the geometry impacts the amount of kinetic energy that is dissipated. Further studies would be needed to understand whether some particular numerical configurations are more prone to generate heat via turbulent dissipation than others among those employed in CCSN simulations.

## 8 CONCLUSION

A physical and numerical experiment was conducted to investigate the turbulent regime of neutrino-driven convection in the gain region of a CCSN. Using an idealized model where heating and gravity are localized in the central region of the domain, we addressed several aspects of the competition between convection and downward advection. Parametrized numerical simulations were employed to assess the robustness of linear predictions developed by Foglizzo et al. (2006) beyond their validity domain. A comparison between 2D and 3D convection enabled us to identify hydrodynamical effects that may result in enhanced heating in global 3D simulations. This includes a wider mixing region due to a faster rise of buoyant bubbles and stronger viscous heating processes leading to higher entropy values in 3D. In addition, we tested several interpretations of the impact of the dimensionality that are based on hydrodynamical considerations because our setup is suitable to isolate the post-shock hydrodynamics from the complexity of the CCSN modelling. A resolution study has been performed in order to anticipate whether future better resolved global simulations would lead to easier explosions. Our main results and their consequences for CCSN simulations can be summarized as follows.

(1) A set of parametric simulations has been run to characterize the impact of a large amplitude perturbation on a flow in which advection can stabilize convection. We confirm that criteria relying on a balance between buoyancy and gravity or drag force (e.g. Scheck et al. 2008; Fernández et al. 2014) predict rather accurately the minimal density contrast above which a buoyant bubble is able to rise against advection. However, if the flow has an initial $\chi$ parameter below the linear instability threshold, we observe that convection is eventually suppressed and the flow returns to its initial configuration, regardless of the perturbation amplitude. The damping time-scale of the temporary convection increases with larger perturbations and proximity to the linear instability threshold. This indicates that criteria applied to a buoyant bubble do not necessarily predict the development of self-sustained turbulent convection. Our study suggests that the criterion established in the linear regime by Foglizzo et al. (2006) is also valid in the non-linear phase. In situations where convection is linearly stable, non-linear turbulent convection can be triggered by a strong perturbation only if a continuous source of excitation exists, such as SASI or numerical errors.

(2) Our model emphasizes that positive 3D effects may already play an important role in the early non-linear regime of convection. In that phase, 3D bubbles rise faster against the flow than their 2D counterparts. Despite the turbulent drag in 3D, buoyant bubbles

resist better against advection in 3D than in 2D. This result is consistent with conclusions drawn from studies on the Rayleigh-Taylor instability with incompressible flows (Young et al. 2001; Anuchina et al. 2004) or Rayleigh-Taylor mixing during the supernova (Kane et al. 2000; Hammer et al. 2010). These results contradict arguments based on the bubble geometry which are thought to favour 2D explosions. A possible explanation is that the velocity of ascending bubbles is mostly set by their radial extension which does not strongly depend on the dimensionality until the saturation of the instability is reached.

(3) For each of the instability regimes identified in our work, much larger average entropy values are measured in 3D than in 2D (Fig. 9). This is a clear sign of an additional source of heating acting in 3D. Part of the discrepancies between 2D and 3D may be ascribed to the differences in terms of dynamics. Large-scale vortices are a natural outcome of 2D simulations and those tend to stir different entropy phases, which is a less efficient mixing process than small-scale mixing. In 2D, the downflows remain unperturbed throughout the gain layer and can channel some matter outside of this region. Turbulent mixing stabilizes a greater fraction of the gain layer in 3D and brings the flow below the marginal stability limit, which is not the case in 2D. Applying a mean-field decomposition on the entropy equation, we identify that turbulent dissipation of kinetic energy could explain the additional heating seen in our 3D simulations. Our results seem to confirm the conjecture of Mabanta & Murphy (2018) that this viscous heating process could play a positive role in multidimensional simulations compared to 1D. These even provide a first evidence that entropy generation is stronger in 3D where the turbulent energy cascade favours smaller spatial scales in which kinetic energy is dissipated into heat.

(4) The impact of resolution is found to be stronger when convection is triggered by a large amplitude perturbation (Fig. 16). In such situations, the entropy production rate increases with dimensionality and resolution and the temporary convection is damped on a longer time-scale. These results suggest that global 3D simulations with improved numerical resolution might benefit from this mechanism. This could be particularly promising for simulations including pre-collapse inhomogeneities which are still in their infancy (Couch & Ott 2015; Müller et al. 2017). This effect is more uncertain when perturbations are numerical artefacts. While this effect would suggest that the discrepancies between Cartesian and spherical grids could widen with finer resolution, this is expected to be counterbalanced by the fact that the numerical noise decreases with resolution. The dynamics seems barely affected by the resolution in the linear instability regime, except for an artificial earlier onset of the instability at low resolution. The level of turbulent dissipation diminishes with higher resolution but seems to be converged with an accuracy of about 20 per cent for a grid of 64 vertical cells in the gain layer. The parametric study of Radice et al. (2016) explored resolutions up to six times higher than that and found no clear trend since their highest resolution case quantitatively resembles their lowest one. So far, convergence studies tend to show that explosions are delayed when the numerical resolution is improved (Hanke et al. 2012; Abdikamalov et al. 2015). It is thus unclear whether higher resolution in CCSN simulations would ease explosions above a certain threshold. It may be that other modelling ingredients would have a stronger impact than resolution. We suggest that the inclusion of large pre-collapse asymmetries could be one of them.

More realistic studies did not exhibit such a strong dependence of the dimensionality on the average entropy in the gain layer,
reaching at most an increase of about 10 per cent in 3D (Dolence et al. 2013). However, we note that in that study, 3D explosions occur earlier than in 2D and this is not the case in the work of Hanke et al. (2012) for which the gap in average entropy values is narrower. It is likely that our setup overestimates the discrepancies between 2D and 3D. A main difference may result from the absence of a stratified cooling layer, in particular in the 2D simulations. Indeed, the vortices are rather expelled out of the gain layer than dissipated inside (Murphy & Meakin 2011) and the flow cannot adjust itself to a sub-critical one as expected from the study of Fernández et al. (2014). The 3D dynamics might be less affected by the specificities of our setup because the turbulent mixing occurs at smaller spatial scales. Nevertheless, our setup contains the minimal set of ingredients for convection to develop in the presence of advection. Our model correctly mimics the properties of turbulence seen in complex CCSN simulations such as the turbulent energy cascade, the anisotropic turbulence and the higher turbulent pressure in 2D. It enabled us to isolate the hydrodynamics from the complexity of the modelling of CCSNe. We found that several arguments commonly employed in CCSN theory to characterize the greater difficulty to achieve explosions in 3D fail at explaining the discrepancies between 2D in 3D in our simplified model.

SASI was mostly neglected in this study except to design the initial perturbation. We note that for SASI-dominated cases, idealized models indicate that 3D spiral modes seem more promising to lead to an explosion than 2D sloshing modes because of the larger non-radial kinetic energy stored in spiral modes (Fernández 2015; Kazeroni et al. 2016). However, more realistic setups produce delayed explosions in 3D non-rotating SASI-dominated models (Melson et al. 2015b) or even failures (Hanke et al. 2013). A possibility is that such an additional energy reservoir in 3D SASI models only supports an explosion in conjunction with other physical ingredients such as rotation. Indeed, stellar rotation rates can enhance spiral modes, but make explosions harder in 2D (Iwakami, Nagakura & Yamada 2014; Kazeroni, Guilet & Foglizzo 2017; Summa et al. 2018). Newly explored microphysics effects could also create more favourable conditions for the development of SASI (Melson et al. 2015b; Bollig et al. 2017). In convection-dominated cases, we speculate that turbulent mixing and dissipation might play a similar role as spiral modes for SASI in that they may favour an explosion but their impact alone could be too modest. Melson et al. (2015a) and Müller (2015) pointed out that favourable 3D effects could contribute to the explosion of convection-dominated cases after shock revival. These effects rely on similar physical processes as in our study: forward turbulent cascade, downflow braking due to the growth of the Kelvin-Helmholtz instability in 3D, more efficient turbulent mixing in 3D. Whether these effects are generic is not yet known.

Finally, a toy model has been used to identify several effects related to the 3D nature of the hydrodynamics that could enhance the heating at the onset of explosion of massive stars. These effects mostly hold for convection-dominated post-shock dynamics, although strong SASI activity can trigger secondary convection. Our idealized model challenges the understanding of convection and the common interpretations of the impact of dimensionality on CCSNe. Our results provide some insights towards a full description of the 3D turbulence generated by neutrino-driven convection. A definitive confirmation of the relevance of the physical processes highlighted in our study will only come from a large set of self-consistent 3D simulations which will become feasible within the next decade.

## ACKNOWLEDGEMENTS

We thank Marc Joos for his help with the analysis tools. We acknowledge insightful discussions with Bernhard Müller, Thomas Janka, Jeremiah Murphy, and Luc Dessart. We are thankful to the referee for helping us improve the manuscript. This work is part of ANR funded project SN2NS ANR-10-BLAN-0503. JG acknowledges support from the European Research Council (grant no. 715368 – MagBURST). This work was granted access to the High-Performance Computing (HPC) resources of Très Grand Centre de Calcul (TGCC) and (Cntre Informatique National de l'Enseignement Supérieur (CINES) under the allocations t2014047094, x2015047094, and t2016047094 made by Grand Équipement National de Calcul Intensif (GENCI) and to the Hydra cluster at Max Planck Computing and Data Facility.

## REFERENCES

* Abdikamalov E. et al., 2015, ApJ, 808, 70

* Andresen H., Müller B., Müller E., Janka H.-T., 2017, MNRAS, 468, 2032

* Anuchina N. N., Volkov V. I., Gordeychuk V. A., Es'kov N. S., Ilyutina O. S., Kozyrev O. M., 2004, J. Comput. Appl. Math., 168, 11

* Arzoumanian Z., Chernoff D. F., Cordes J. M., 2002, ApJ, 568, 289

* Bethe H. A., Wilson J. R., 1985, ApJ, 295, 14

* Blondin J. M., Mezzacappa A., DeMarino C., 2003, ApJ, 584, 971

* Bollig R., Janka H.-T., Lohs A., Martínez-Pinedo G., Horowitz C. J., Melson T., 2017, Phys. Rev. Lett., 119, 242702

* Burrows A., Goshy J., 1993, ApJ, 416, L75

* Burrows A., Hayes J., Fryxell B. A., 1995, ApJ, 450, 830

* Burrows A., Dolence J. C., Murphy J. W., 2012, ApJ, 759, 5

* Cabot W., 2006, Phys. Fluids (1994–present), 18, 045101

* Cardall C. Y., Budiardja R. D., 2015, ApJ, 813, L6

* Couch S. M., 2013, ApJ, 775, 35

* Couch S. M., O'Connor E. P., 2014, ApJ, 785, 123

* Couch S. M., Ott C. D., 2013, ApJ, 778, L7

* Couch S. M., Ott C. D., 2015, ApJ, 799, 5

* Dolence J. C., Burrows A., Murphy J. W., Nordhaus J., 2013, ApJ, 765, 110

* Fernández R., 2015, MNRAS, 452, 2071

* Fernández R., Thompson C., 2009, ApJ, 697, 1827

* Fernández R., Müller B., Foglizzo T., Janka H.-T., 2014, MNRAS, 440, 2763

* Foglizzo T. et al., 2015, PASA, 32, 9

* Foglizzo T., 2009, ApJ, 694, 820

* Foglizzo T., Scheck L., Janka H.-T., 2006, ApJ, 652, 1436

* Foglizzo T., Galletti P., Scheck L., Janka H.-T., 2007, ApJ, 654, 1006

* Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371

* Grefenstette B. W. et al., 2014, Nature, 506, 339

* Guilet J., Foglizzo T., 2012, MNRAS, 421, 546

* Guilet J., Sato J., Foglizzo T., 2010, ApJ, 713, 1350

* Hammer N. J., Janka H.-T., Müller E., 2010, ApJ, 714, 1371

* Handy T., Plewa T., Odrzywołek A., 2014, ApJ, 783, 125

* Hanke F., Marek A., Müller B., Janka H.-T., 2012, ApJ, 755, 138

* Hanke F., Müller B., Wongwathanarat A., Marek A., Janka H.-T., 2013, ApJ, 770, 66

* Herant M., Benz W., Colgate S., 1992, ApJ, 395, 642

* Herant M., Benz W., Hix W. R., Fryer C. L., Colgate S. A., 1994, ApJ, 435, 339

* Iwakami W., Nagakura H., Yamada S., 2014, ApJ, 793, 5

* Janka H.-T., 2001, A&A, 368, 527

* Janka H.-T., Mueller E., 1996, A&A, 306, 167

* Janka H.-T., Melson T., Summa A., 2016, Annu. Rev. Nucl. Part. Sci., 66, 341

* Kane J., Arnett D., Remington B. A., Glendinning S. G., Bazán G., Müller E., Fryxell B. A., Teyssier R., 2000, ApJ, 528, 989

* Kazeroni R., Guilet J., Foglizzo T., 2016, MNRAS, 456, 126

* Kazeroni R., Guilet J., Foglizzo T., 2017, MNRAS, 471, 914
Kitaura F. S., Janka H., Hillebrandt W., 2006, A&A, 450, 345

Kraichnan R. H., 1967, Phys. Fluids, 10, 1417

Kuroda T., Takiwaki T., Kotake K., 2016, ApJS, 222, 20

Landau L. D., Lifshitz E. M., 1959, Theory of Elasticity. Pergamon Press, Reading, MA

Lentz E. J. et al., 2015, ApJ, 807, L31

Leonard D. C. et al., 2006, Nature, 440, 505

Liebendorfer M., Mezzacappa A., Thielemann F.-K., Messer O. E., Hix W. R., Bruenn S. W., 2001, Phys. Rev. D, 63, 103004

Lund T., Wongwathanarat A., Janka H.-T., Müller E., Raffelt G., 2012, Phys. Rev. D, 86, 105031

Mabanta Q. A., Murphy J. W., 2018, ApJ, 856, 22

Marek A., Janka H.-T., 2009, ApJ, 694, 664

Melson T., Janka H.-T., Marek A., 2015a, ApJ, 801, L24

Melson T., Janka H.-T., Bollig R., Hanke F., Marek A., Müller B., 2015b, ApJ, 808, L42

Miyoshi T., Kusano K., 2005, J. Comput. Phys., 208, 315

Mocák M., Meakin C., Viallet M., Arnett D., 2014, preprint (arXiv:1401.176)

Müller B., 2015, MNRAS, 453, 287

Müller B., 2016, PASA, 33, e048

Müller B., Janka H.-T., 2015, MNRAS, 448, 2141

Müller B., Janka H.-T., Heger A., 2012, ApJ, 761, 72

Müller B., Janka H.-T., Marek A., 2013, ApJ, 766, 43

Müller B., Viallet M., Heger A., Janka H.-T., 2016, ApJ, 833, 124

Müller B., Melson T., Heger A., Janka H.-T., 2017, MNRAS, 472, 491

Murphy J. W., Burrows A., 2008, ApJ, 688, 1159

Murphy J. W., Dolence J. C., 2017, ApJ, 834, 183

Murphy J. W., Meakin C., 2011, ApJ, 742, 74

Murphy J. W., Dolence J. C., Burrows A., 2013, ApJ, 771, 52

O'Connor E., Ott C. D., 2011, ApJ, 730, 70

Ott C. D. et al., 2013, ApJ, 768, 115

Radice D., Couch S. M., Ott C. D., 2015, Comput. Astrophys. Cosmol., 2, 7

Radice D., Ott C. D., Abdikamalov E., Couch S. M., Haas R., Schnetter E., 2016, ApJ, 820, 76

Roberts L. F., Ott C. D., Haas R., O'Connor E. P., Diener P., Schnetter E., 2016, ApJ, 831, 98

Scheck L., Plewa T., Janka H.-T., Kifonidis K., Müller E., 2004, Phys. Rev. Lett., 92, 011103

Scheck L., Kifonidis K., Janka H.-T., Müller E., 2006, A&A, 457, 963

Scheck L., Janka H.-T., Foglizzo T., Kifonidis K., 2008, A&A, 477, 931

Summa A., Hanke F., Janka H.-T., Melson T., Marek A., Müller B., 2016, ApJ, 825, 6

Summa A., Janka H.-T., Melson T., Marek A., 2018, ApJ, 852, 28

Takiwaki T., Kotake K., Suwa Y., 2014, ApJ, 786, 83

Tamborra I., Hanke F., Müller B., Janka H.-T., Raffelt G., 2013, Phys. Rev. Lett., 111, 121104

Tamborra I., Raffelt G., Hanke F., Janka H.-T., Müller B., 2014a, Phys. Rev. D, 90, 045032

Tamborra I., Hanke F., Janka H.-T., Müller B., Raffelt G. G., Marek A., 2014b, ApJ, 792, 96

Teyssier R., 2002, A&A, 385, 337

Wongwathanarat A., Janka H.-T., Müller E., 2013, A&A, 552, A126

Wongwathanarat A., Janka H.-T., Müller E., Pllumbi E., Wanajo S., 2017, ApJ, 842, 13

Young Y.-N., Tufo H., Dubey A., Rosner R., 2001, J. Fluid Mech., 447, 377

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

Fig6.mp4

Fig13_time_simulation.mp4

Fig13_flyby.mp4

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by

the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: MEAN-FIELD DECOMPOSITION OF THE ENTROPY EQUATION

Using a Frevian decomposition of S and $v_z$ for the left-hand side of the entropy equation (19), one obtains:

$$ \frac{\partial (\rho S)}{\partial t} + \nabla \cdot (\rho S v_z) = \frac{\partial (\rho \tilde{S})}{\partial t} + \frac{\partial (\rho S'')}{\partial t} + \nabla \cdot (\rho \tilde{S} v_z) + \nabla \cdot (\rho S'' \tilde{v}_z) + \nabla \cdot (\rho S'' v_z'') \quad (A1) $$

The definition of a Frevian average is such that: $\overline{\rho\tilde{\psi}} = \overline{\rho\tilde{\psi}} = \overline{\rho\psi}$. By definition, the average of a Frevian fluctuation cancels: $\overline{\rho\psi''} = \rho(\psi - \tilde{\psi}) = \rho\tilde{\psi} - \overline{\rho\tilde{\psi}} = 0$. Consequently, the average of the following two terms cancels: $\partial(\rho S'')/\partial t$ and $\nabla \cdot (\rho S'' \tilde{v}_z)$.

Using the mean-field decomposition (A1) to the averaged entropy equation (19), one obtains:

$$ \frac{\partial \overline{\rho \tilde{S}}}{\partial t} + \nabla \cdot (\overline{\rho \tilde{v}_z \tilde{S}}) + \nabla \cdot (\overline{\rho v_z'' S''}) = \left(\frac{\overline{\rho \dot{q}}}{T}\right) + \left(\frac{\overline{\rho \epsilon}}{T}\right) \quad (A2) $$

## APPENDIX B: ANALYTICAL SOLUTION OF THE NEUTRAL STABILITY IN THE ASYMPPTOTIC LIMIT $K_G \ll 1, K_H \ll 1, M_{up} \ll 1$

### B1 Differential system

The evolution of perturbations is described by the same differential system and boundary conditions as in Foglizzo et al. (2006) considering two cases for the upper boundary condition: either a shock condition, or an acoustic leaking condition. This latter condition expresses that acoustic waves are free to propagate upward and that no acoustic flux nor advected perturbations are injected downward from the region above the heating region.

In order to reach an analytical characterization of the heating threshold for a convective instability, we chose a regime where the flow is approximately uniform, thus allowing for an analytical integration of the differential system. This is achieved in the limit of small gravity ($K_G \ll 1$) and small external heating ($K_H \ll 1$). In this regime, having comparable advection and buoyancy timescales (i.e. finite non-zero $\chi$ parameter) requires a small enough mach number ($M_{up}^2 \propto K_H K_G \ll 1$). This asymptotic differential system is obtained by introducing some new functions ($\delta\tilde{F}, \delta\tilde{h}$) defined by

$$ \delta\tilde{F} \equiv \frac{1}{M} \left( \frac{\delta f}{c^2} - \frac{\delta S}{\gamma} \right), \quad (B1) $$

$$ \delta\tilde{h} \equiv M \left( \delta h + \frac{\gamma - 1}{\gamma} \delta S \right). \quad (B2) $$

where $\delta f \equiv v\delta v_z + 2c\delta c(\gamma - 1)$ and $\delta h \equiv \delta v_z/v + \delta\rho/\rho$. Conversely,

$$ \frac{\delta \rho}{\rho} = \frac{1}{1 - M^2} \left( M(\delta \tilde{F} - \delta \tilde{h}) - \frac{\gamma - 1}{\gamma} (1 - M^2) \delta S \right), \quad (B3) $$

$$ \frac{\delta v_z}{v} = \frac{1}{1 - M^2} \frac{1}{M} (\delta \tilde{h} - M^2 \delta \tilde{F}), \quad (B4) $$

$$ \frac{\delta c^2}{c^2} = \frac{\gamma - 1}{1 - M^2} \left( M(\delta \tilde{F} - \delta \tilde{h}) + (1 - M^2) \frac{\delta S}{\gamma} \right). \quad (B5) $$
The differential system satisfied by $(\delta \hat{F}, \delta \hat{h}, \delta S, \delta v_x)$ for neutral stability ($\omega = 0$) is the following

$$ \frac{\partial \delta \hat{F}}{\partial z} = \frac{\gamma - 1}{1 - M^2} \frac{\nabla \Phi}{Mc^2} \left( \frac{\delta S}{\gamma} + M \delta \hat{F} \right) - \frac{1}{M} \frac{\partial M}{\partial z} \delta \hat{F} \\ - \left( \delta \hat{h} - \frac{\gamma - 1}{\gamma} M \delta S - \gamma M^2 \delta \hat{F} \right) \frac{1}{1 - M^2} \frac{\gamma - 1}{\gamma} \nabla S, \quad (B6) $$

$$ \frac{\partial \delta \hat{h}}{\partial z} = k_x \frac{i \delta v_x}{c} + \frac{\gamma - 1}{\gamma} M \delta \left( \frac{L}{pv} \right) + \frac{1}{M} \frac{\partial M}{\partial z} \delta \hat{h}, \quad (B7) $$

$$ \frac{\partial \delta S}{\partial z} = \delta \left( \frac{L}{pv} \right), \quad (B8) $$

$$ \frac{\partial i \delta v_x}{\partial z} = -k_x c \delta \hat{F} + k_x \frac{c}{1 - M^2} (\delta \hat{h} - M^2 \delta \hat{F}), \quad (B9) $$

When $M \ll 1$, $H|\nabla S| \ll 1$, and $H\nabla\Phi/c^2 \ll 1$, the flow is approximately uniform and the leading terms of the differential system are independent of the detailed form of the heating function. For analytical simplicity, we neglect the linear transition ramp used to smoothen the gravity profile and the heating function ($\Psi(|z| < H) = 1$) such that $K_H = -H\nabla S$ and $K_G = H\nabla\Phi/c_{up}^2$. Assuming $K_H \propto M_{up}$ and $K_G \propto M_{up}$ in the limit where $M_{up} \ll 1$, the leading coefficients of the differential system are constants of order unity:

$$ \frac{\partial \delta \hat{F}}{\partial z} = \frac{\gamma - 1}{\gamma} \frac{K_G}{HM_{up}} \delta S, \quad (B10) $$

$$ \frac{\partial \delta \hat{h}}{\partial z} = k_x \frac{i \delta v_x}{c_{up}}, \quad (B11) $$

$$ \frac{\partial \delta S}{\partial z} = \frac{K_H}{HM_{up}} \delta \hat{h}, \quad (B12) $$

$$ \frac{\partial}{\partial z} \frac{i \delta v_x}{c_{up}} = k_x (\delta \hat{h} - \delta \hat{F}). \quad (B13) $$

The eigenfunctions $(\delta \hat{F}, \delta \hat{h}, \delta S, i \delta v_x/c_{up})$ are a linear combination of four functions $\exp(\beta_j \chi k_x z)$ with respective coefficients $(f_j, h_j, s_j, v_j)_{j=1,4}$. An important property of the solution is that the parameters $K_H, K_G$, and $M_{up}$ play a role only through the combination $\chi$ defined by

$$ \chi^2 = 4 \frac{\gamma - 1}{\gamma} \frac{K_H K_G}{M_{up}^2}, \quad (B14) $$

$$ = \left( \frac{\omega_{BV} 2H}{v_{up}} \right)^2. \quad (B15) $$

$\beta_{j=1,4}$ are the four roots of the equation $\beta^4 - \beta^2 + 1/4\alpha^2 = 0$ (deduced from the determinant of the system of equations B10–B13) with

$$ \alpha \equiv \frac{H k_x}{\chi}. \quad (B16) $$

If $0 \le \alpha \le 1$, we define the angle $\phi$ with $0 \le \phi \le \pi/2$ such that $\alpha = \cos \phi$:

$$ \beta_1 \equiv \frac{e^{i \frac{\phi}{2}}}{(2 \cos \phi)^{\frac{1}{2}}}, \quad (B17) $$

$$ \beta_2 \equiv \bar{\beta}_1 = \frac{e^{-i \frac{\phi}{2}}}{(2 \cos \phi)^{\frac{1}{2}}}, \quad (B18) $$

$$ \beta_3 \equiv -\beta_1, \quad (B19) $$

$$ \beta_4 \equiv -\beta_2. \quad (B20) $$

## B2 Boundary conditions

### B2.1 Decomposition into advected and acoustic waves

In the uniform semi-infinite regions above and below the region of heating, the perturbation is decomposed into four components: entropic, vortical, and two acoustic waves propagating with the flow (i.e. downward, subscript +) or against the flow (i.e. upward, subscript -)

$$ \delta \hat{F} = \delta \hat{F}_S + \delta \hat{F}_w + \delta \hat{F}_+ + \delta \hat{F}_-. \quad (B21) $$

The identification of acoustic waves propagating can be adapted from equations (C16)–(C21) in Foglizzo et al. (2006). The pressure perturbation $(\delta \hat{F}_{\pm}, \delta \hat{h}_{\pm})$ with a vertical wavenumber $k_z^{\pm}$ is characterized by the absence of perturbation of entropy $(\delta S = 0)$ and vorticity $(k_x c^2 \delta \hat{F}_{\pm} = 0)$. With $k_z^{\pm} = \mp i k_x$ in the low Mach limit, the pressure profile of a purely growing mode $\exp(ik_z^+ z)$ is evanescent upward while $\exp(ik_z^- z)$ is evanescent downward. Using this formula for the expression of $\delta \hat{h}_{\pm}$,

$$ \delta \hat{F}_{\pm} = 0, \quad (B22) $$

$$ \delta \hat{h}_{\pm} = \pm i \frac{\delta v_x^{\pm}}{c}. \quad (B23) $$

The advected perturbation containing entropy and vorticity perturbations is purely advected ($\omega = k_z v$). From equations (B10)–(B13) with $k_z = \omega/v = 0$ for neutral stability, we deduce:

$$ \delta v_x \text{ adv} = 0, \quad (B24) $$

$$ \delta \hat{F}_{\text{adv}} = \delta \hat{h}_{\text{adv}}. \quad (B25) $$

### B2.2 Leaking lower boundary

The lower boundary condition $\delta \hat{h}_{-}^{d} = 0$ is expressed in the low Mach limit:

$$ \delta \hat{F}_{d} - \delta \hat{h}_{d} + \frac{i \delta v_{x}^{d}}{c_{d}} = 0. \quad (B26) $$

### B2.3 Leaking upper boundary condition without a shock

The upper boundary condition corresponds to a pure acoustic wave propagating against the flow: $\delta S = 0$, $\delta F_{\text{adv}} = 0$, and $\delta F_{+} = 0$. In the low Mach limit,

$$ \delta S_{\text{up}} = 0, \quad (B27) $$

$$ \delta \hat{F}_{\text{up}} = 0, \quad (B28) $$

$$ \delta \hat{h}_{\text{up}} = - \frac{i \delta v_x^{\text{up}}}{c_{\text{up}}}. \quad (B29) $$
<figure><img src="image_20.png" /><figcaption>**Figure B1.** Convergence of the instability threshold $\chi(H_{k_x})$ towards the analytical expression (dashed line) when $K_G \ll 1$, either in a subsonic flow ($M_{up} \ll 1$, left-hand panel) or in a post-shock flow ($M_{sh} \ll 1$, right-hand panel). The curves have been calculated without a transition function $\Psi(|z| < 1) = 1$, except for the thick curve noted $K_G = 3L$ which shows the incidence of a linear ramp $\Psi(|z| < 0.5) = 1$. The analytical solutions correspond to equations (B44) and (B61).</figcaption></figure>

## B2.4 Upper boundary condition with a shock

The boundary condition for a perturbed shock with dissociation established by Foglizzo et al. (2006, equations 28–30 and C14), are written using the functions $\delta\hat{F}$, $\delta\hat{h}$, at marginal stability. In the uniform limit where $M_{sh} \ll 1$, $H\nabla S \ll 1$, and $H\nabla\Phi/c_{sh}^2 \ll 1$,

$$ \delta\hat{h}_{\text{sh}} = 0, \quad (B30) $$

$$ \delta S_{\text{sh}} = 0, \quad (B31) $$

$$ \delta \hat{F}_{\text{sh}} = 0. \quad (B32) $$

# B3 Analytical formulation of the criterion for marginal stability

## B3.1 Marginal stability without a shock

With $z_d = -H$, the differential system is translated into linear relations between the coefficients ($f_j$, $h_j$, $s_j$, $v_j$)$_{j=1,4}$:

$$ f_j = \frac{1}{4\alpha^2 \beta_j^2} h_j, \quad (B33) $$

$$ s_j = \frac{K_H}{M \alpha \beta_j \chi} h_j, \quad (B34) $$

$$ v_j = \beta_j h_j. \quad (B35) $$

The boundary conditions are translated into:

$$ \sum_{j=1}^{4} s_j = 0, \quad (B36) $$

$$ \sum_{j=1}^{4} f_j = 0, \quad (B37) $$

$$ \sum_{j=1}^{4} (v_j + h_j) = 0, \quad (B38) $$

$$ \sum_{j=1}^{4} (f_j - h_j + v_j) e^{-2\chi\alpha\beta_j} = 0. \quad (B39) $$

The value of $\alpha(\chi)$ is set by the vanishing determinant of the linear system satisfied by ($h_1$, $h_2$, $h_3$, $h_4$):

$$ \left| \begin{array}{cccc} \beta_2 & \beta_1 & -\beta_2 & -\beta_1 \\ \beta_2^2 & \beta_1^2 & \beta_2^2 & \beta_1^2 \\ 1 + \beta_1 & 1 + \beta_2 & 1 - \beta_1 & 1 - \beta_2 \\ (\beta_1^2 - \beta_1) e^{-2\chi\alpha\beta_1} & (\beta_2^2 - \beta_2) e^{-2\chi\alpha\beta_2} & (\beta_1^2 + \beta_1) e^{2\chi\alpha\beta_1} & (\beta_2^2 + \beta_2) e^{2\chi\alpha\beta_2} \end{array} \right| = 0, \quad (B40) $$

where we have used $1/\beta_1 = (2\cos\phi)\beta_2$. The determinant $\Delta_{\text{leak}}$ can be factorized as follows

$$ \Delta_{\text{leak}} = \frac{8}{4\alpha^2} (\beta_2^2 - \beta_1^2)(Y_1 - Y_2) \left(1 - \frac{1}{Y_1 Y_2}\right), \quad (B41) $$

$$ Y_1 = 4\alpha^2(\beta_1^2 + \beta_1)^2 e^{2\chi\alpha\beta_1}, \quad (B42) $$

$$ Y_2 = 4\alpha^2(\beta_2^2 + \beta_2)^2 e^{2\chi\alpha\beta_2}. \quad (B43) $$

It can be shown that the only vanishing factor in the factorized expression of $\Delta_{\text{leak}}$ is $Y_1 - Y_2$. For $\alpha(\chi) < 1$, the analytical relation between the value of $\chi$ and $H_{k_x}$ is thus reduced to:

$$ \chi = \frac{2^{\frac{1}{2}}}{\cos^{\frac{1}{2}} \phi \sin^{\frac{\phi}{2}}} \arctan \frac{\cos^{\frac{1}{2}} \phi}{2^{\frac{1}{2}} \sin^{\frac{\phi}{2}}}, \quad (B44) $$

$$ \cos \phi \equiv \frac{H k_x}{\chi}. \quad (B45) $$
The solution $Hk_x = \chi$ when $\chi \gg 1$ corresponds to a stabilization a short wavelength when

$$k_x = \frac{2}{|v|} \omega_{BV}(\chi \gg 1). \qquad (B46)$$

It agrees within a factor 2 with the intuitive expectation based on the evanescent vertical profile of the Rayleigh-Taylor instability (exp($\omega_i - k_x|v|$)t) (Guilet et al. 2010). The limit $\phi = \pi/2$ corresponds to the threshold of stability when $k_x \to 0$ and $\chi \to 2$. The flow is linearly stable for $\chi < 2$ (Fig. B1, left-hand panel).

## B4 Marginal stability with a shock

The boundary conditions are translated into:

$$\sum_{j=1}^{4} f_j = 0, \qquad (B47)$$

$$\sum_{j=1}^{4} s_j = 0, \qquad (B48)$$

$$\sum_{j=1}^{4} h_j = 0, \qquad (B49)$$

$$\sum_{j=1}^{4} (f_j - h_j + v_j) e^{-2\chi\alpha\beta_j} = 0, \qquad (B50)$$

with

$$f_j = \frac{1}{4\alpha^2\beta_j^2} v_j, \qquad (B51)$$

$$s_j = \frac{K_H}{M\alpha\beta_j^2\chi} v_j, \qquad (B52)$$

$$h_j = \frac{1}{\beta_j} v_j. \qquad (B53)$$

The value of $\alpha(\chi)$ is set by the vanishing determinant of the linear system satisfied by $(v_1, v_2, v_3, v_4)$:

$$\left| \begin{array}{cccc} \beta_2 & \beta_1 & -\beta_2 & -\beta_1 \\ \beta_2^2 & \beta_1^2 & \beta_2^2 & \beta_1^2 \\ \beta_3^2 & \beta_1^3 & -\beta_2^3 & -\beta_1^3 \\ (1 - \beta_1)e^{-2\chi\alpha\beta_1} & (1 - \beta_2)e^{-2\chi\alpha\beta_2} & (1 + \beta_1)e^{2\chi\alpha\beta_1} & (1 + \beta_2)e^{2\chi\alpha\beta_2} \end{array} \right| = 0 \qquad (B54)$$

Defining

$$Y_3 = (1 + \beta_1)e^{2\chi\alpha\beta_1} + (1 - \beta_1)e^{-2\chi\alpha\beta_1}, \qquad (B55)$$

$$Y_4 = (1 + \beta_2)e^{2\chi\alpha\beta_2} + (1 - \beta_2)e^{-2\chi\alpha\beta_2}, \qquad (B56)$$

the determinant $\Delta_{\text{shock}}$ is factorized as follows:

$$\Delta_{\text{shock}} = 8\beta_1\beta_2(\beta_1^2 - \beta_2^2)(\beta_1^2 Y_3 - \beta_2^2 Y_4) \qquad (B57)$$

For $\alpha = Hk_x/\chi < 1$, we transform the equation $\beta_1^2 Y_3 - \beta_2^2 Y_4 = 0$ into an implicit equation relating $\chi$ and $Hk_x$:

$$\cos \phi \equiv \frac{Hk_x}{\chi}, \qquad (B58)$$

$$X_c \equiv \chi \times (2 \cos \phi)^{\frac{1}{2}} \cos \frac{\phi}{2}, \qquad (B59)$$

$$X_s \equiv \chi \times (2 \cos \phi)^{\frac{1}{2}} \sin \frac{\phi}{2}, \qquad (B60)$$

$$(2 \cos \phi)^{\frac{1}{2}} (\cos \phi \sin X_s \tanh X_c + \sin \phi \cos X_s) + \cos \frac{3\phi}{2} \sin X_s + \sin \frac{3\phi}{2} \cos X_s \tanh X_c = 0. \qquad (B61)$$

This equation has two solutions $k_{\min}$, $k_{\max}$ when $\chi > \chi_{\text{crit}} \simeq 3.291$. Using a Taylor expansion of $\tanh X_c$ and $\tan X_s$ for $\chi \gg 1$, $Hk_{\min} \ll 1$, $\phi \sim \pi/2$, $\alpha \ll 1$, we obtain

$$Hk_{\min} = \frac{3}{\chi^2}. \qquad (B62)$$

If $\chi \gg 1$ and $Hk_{\max} \gg 1$, $\phi \ll 1$, $\alpha = 1$,

$$Hk_{\max} = \chi. \qquad (B63)$$

The threshold for global instability corresponds to $Hk_{\min} = Hk_{\max} \sim 1$ and $\chi_{\text{crit}} \simeq 3.291$ (Fig. B1, right-hand panel).

This paper has been typeset from a TeX/LA TeX file prepared by the author.