{'fraction_non_alphanumeric': 0.03183003380009657, 'fraction_numerical': 0.016668276195074844, 'mean_word_length': 4.460451381565071, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The impact of progenitor asymmetries on the neutrino-driven convection in core-collapse supernovae

The explosion of massive stars in core-collapse supernovae may be aided by the convective instabilities that develop in their innermost nuclear burning shells. The resulting fluctuations support the explosion by generating additional turbulence behind the supernova shock. It was suggested that the buoyant density perturbations arising from the interaction of the pre-collapse asymmetries with the shock may be the primary contributor to the enhancement of the neutrino-driven turbulent convection in the post-shock region. Employing three-dimensional numerical simulations of a toy model, we investigate the impact of such density perturbations on the post-shock turbulence. We consider a wide range of perturbation parameters. The spatial scale and the amplitude of the perturbations are found to be of comparable importance. The turbulence is particularly enhanced when the perturbation frequency is close to that of the convective turnovers in the gain region. Our analysis confirms that the buoyant density perturbations is indeed the main source of the additional turbulence in the gain region, validating the previous order-of-magnitude estimates.


INTRODUCTION
Multi-dimensional hydrodynamical effects are crucial to power the explosion of massive stars, known as core-collapse supernovae (CCSNe). The gravitational contraction of the iron core of the progenitor into a proto-neutron star triggers a bounce of the core which produces a stalled shock wave at a radius of ∼ 150 km. A successful explosion depends on the revival of this shock by neutrino heating during the first second after the core bounce (see Foglizzo et al. 2015;Janka et al. 2016;Müller 2016 for recent reviews). Except for the lightest progenitors, explosions of massive stars fail in spherical symmetry (Liebendörfer et al. 2001;Kitaura et al. 2006). The shock revival relies on the decisive action of multi-dimensional fluid instabilities in the post-shock region, such as neutrino-driven convection (Herant et al. 1992(Herant et al. , 1994Burrows et al. 1995;Janka & Mueller 1996;Foglizzo et al. 2006) and the Standing Accretion Shock Instability (SASI) (Blondin et al. 2003;Foglizzo et al. 2007;Foglizzo 2009;Guilet & Foglizzo 2012). These instabilities generate large-scale turbulent motions which enhance the efficiency of the neutrino heating and ultimately trigger an asymmetric explosion.
The neutrino mechanism is supported by an increasing number of successful three-dimensional (3D) CCSN simulations across a wide range of progenitor masses (Takiwaki et al. 2012(Takiwaki et al. , 2014Melson et al. 2015a,b;Lentz et al. 2015;Roberts et al. 2016;Müller et al. 2017Müller et al. , 2019Ott et al. 2018;Glas et al. 2019;Vartanyan et al. 2019;Burrows et al. 2019). However, it is not yet clear if the 3D explosions are powerful enough to reproduce observables such as the explosion energy or the distribution of neutron star masses (e.g. Müller 2016). Various ingredients leading to more robust explosions have been considered, including stellar rotation (Takiwaki et al. 2016;Kazeroni et al. 2017;Summa et al. 2018), modifications in the opacities (Melson et al. 2015b), the inclusion of muons (Bollig et al. 2017) or finer numerical resolution (Nagakura et al. 2019;. Another ingredient comes from the asymmetries in the progenitors which could arise naturally in the innermost burning shells of massive stars prior to the collapse of the iron core. Linear theories suggest that fluctuations can be amplified during the collapse (Kovalenko & Eremin 1998;Lai & Goldreich 2000;Takahashi & Yamada 2014). The accretion of these pre-collapse fluctuations generates additional turbulent pressure in the post-shock medium which pro-vides better conditions for shock revival (Couch & Ott 2013. Aspherical accretion could lead to strong shock deformations which may enhance post-shock turbulence even further (Müller & Janka 2015). Pre-collapse perturbations may reduce the critical neutrino luminosity required for an explosion by up to ∼ 20% (Müller et al. , 2017Radice et al. 2018). Through an extensive set of axisymmetric simulations, Müller & Janka (2015) showed that the pre-collapse convective motions dominated by large spatial scales are particularly helpful to revive the stalled shock. Such conditions may be encountered in the O-burning shell and to some degree in the Si-burning shell for a large range of supernova progenitors (Collins et al. 2018). 3D simulations of pre-collapse shell convection revealed a wide variety of the spatial scales and strengths of these fluctuations Müller et al. 2016Müller et al. , 2019Yoshida et al. 2019;Yadav et al. 2019).
Depending on their parameters, the impact of the progenitor perturbations on shock revival ranges from limited to significant.  obtained velocity fluctuations dominated by modes = 4 − 5 in a simulation of Si-burning performed in octant symmetry. The 3D initial conditions led to a slightly earlier shock revival in the subsequent CCSN simulation with a minor impact on the heating conditions. Müller et al. (2016) obtained a dominant mode = 2 in a simulation of O-burning whose impact was decisive to revive the shock in the subsequent CCSN, since the spherically-symmetric progenitor failed at producing an explosion (Müller et al. 2017).
The physics of the coupling between pre-collapse perturbations and post-shock instabilities still remains to be fully understood. Using linear analysis, Takahashi et al. (2016) showed that the accreting perturbations may resonantly amplify SASI (or neutrino-driven convection) only if the perturbation frequency matches the oscillation frequency of SASI (or the convective growth rate). However, such a fine matching is unlikely to be common in nature. As the pre-collapse convective perturbations descend towards the center of the star, they generate pressure waves (Müller & Janka 2015;Abdikamalov & Foglizzo 2019). The interaction of these perturbations with the supernova shock creates acoustic, entropy, and vorticity waves in the post-shock region Huete et al. 2018;Abdikamalov et al. 2018). Of these three, initially the vorticity waves lead to the largest increase of the non-radial motion in the post-shock region. However, the entropy waves, which are accompanied by density variations (in order to maintain pressure equilibrium with their surrounding), become buoyant and generate additional turbulence (Müller et al. , 2017. Müller et al. (2016) estimated that this is the dominant source of additional turbulence in the gain region.
In this study, we investigate the impact of the buoyant entropic density waves 1 on neutrino-driven convection in the gain region. We perform numerical simulations of an idealized model in which parameterized density perturbations are injected into a stationary flow that is unstable to neutrinodriven convection (Kazeroni et al. 2018, hereafter Paper I). Our idealized model does not contain a shock wave, which allows us to disentangle the impact of the density perturbations on the neutrino-driven convection. Controlled simulations enable us to quantify the impact of the perturbation properties, such as the amplitude, the spatial scale, and the temporal frequency on turbulence. Finally, we perform additional two-dimensional (2D) simulations and examine their limitations by comparing them to our 3D models.
Our results confirm that the injection of turbulent kinetic energy in the gain layer mostly results from the work of buoyancy on the density perturbations. The amount of turbulent kinetic energy generated by the perturbation can well be captured by an analytical estimate. The spatial scale and the amplitude of the perturbation have a comparable impact on the additional turbulent kinetic energy. We find that turbulence can be significantly enhanced if the perturbation frequency is close to that of convective turnovers. While large scale modes are found to maximize the turbulent kinetic energy both in 2D and 3D, the strength of the turbulence is largely over-predicted in 2D due to a different turbulent energy cascade.
The structure of the paper is the following. The methodology is detailed in Section 2. The results of the numerical simulations are described in Section 3. The origin and the properties of the additional turbulence related to the accretion of density perturbations are investigated in Section 4. The discrepancies between 2D and 3D simulations are analyzed in Section 5. We summarize and discuss our findings in Section 6.

Stationary flow
As in Paper I, we consider a stationary flow of an ideal gas with γ = 4/3 along vertical z-direction in a 3D Cartesian geometry. The computational domain covers the region −150 km ≤ x, y ≤ 150 km and −450 km ≤ z ≤ 450 km. The central −H ≤ z ≤ H region, where H = 50 km, mimics the gain layer of a CCSN, where neutrino absorption exceeds neutrino emission (cf. Fig. 1). The heating and the gravity in this layer are modeled as described Paper I. The heating is proportional to density and is normalized by the coefficient K H = 2.76 × 10 −3 , while the gravitational potential ∇Φ is normalized by the coefficient K G = 3 (cf. Paper I for the definitions of these quantities). Outside the gain layer, heating and gravity vanish and the stationary flow is uniform. The vertical boundaries are placed far enough to minimize the impact of reflected acoustic waves. The upstream Mach number is set to Ma up = 0.3, where the subscript "up" refers to values in the upstream flow. As mentioned earlier, in order to study the coupling between neutrino-driven convection and injected perturbations in its simplest form, our model does not contain a shock wave to avoid any feedback that would result from its deformation. In addition, the cooling layer surrounding the proto-neutron star is absent from our model for the sake of simplicity.
For convection to develop in the gain layer, the buoyancy has to operate faster than advection through the gain layer. This condition can be expressed in terms of the parameter χ (Foglizzo et al. 2006), which measures the ratio of the advection timescale through the gain region t adv to the convective turnover timescale t conv : where · refers to horizontally averaged quantities and v z to the vertical velocity. The advection timescale through the gain layer is while convective turnover timescale can be approximated as t conv ∼ ω −1 BV , where the Brunt-Väisälä frequency ω BV is defined as with P, ρ, and S being respectively the pressure, the density, and the dimensionless entropy of the flow. The latter is computed as In our setup, the stationary flow is linearly unstable to neutrino-driven convection only if χ χ crit ≈ 2.4 (Kazeroni et al. 2018) 2 . For our choice of parameters K G , K H , and Ma up , the stationary flow has χ = 3, i.e. t conv ≈ t adv /3, which is slightly above the instability threshold.

Perturbation modelling
We first generate a flow where neutrino-driven convection has reached a fully turbulent stage. A simulation of neutrinodriven convection is performed starting from a vertical stationary flow and using a density perturbation of 0.1% amplitude. The non-linear stage of neutrino-driven convection is reached after 10 advection times. The simulation is run for a total duration of 18 advection times. This simulation is used as a reference model to investigate the impact of injected perturbations on the convective instability.
We then take an output of the reference simulation at a time t = 12.8 t adv and inject a parameterized density perturbation through the upper boundary of domain, located at z = 450 km. At this time, neutrino-driven convection in the gain layer has reached the non-linear regime for about than 9 convective turnover times. The simulations of perturbed convection are evolved from t = 12.8 t adv to t = 18 t adv which covers about 15 convective turnover times. We find that the simulations reach a steady state regime by this time and a further evolution does not bring any different outcome.
The infalling density perturbation is modeled as a planar sinusoidal structure in the xz-plane ( Fig. 1): where ρ up is the the density of the upstream stationary flow. In order to maintain the pressure equilibrium with the surrounding flow, the associated temperature variations have the opposite phase. The parameter A is the amplitude, k x 2 Note that for a more realistic model including a shock wave, the instability threshold is χ crit ≈ 3.2 (Foglizzo et al. 2006).  Figure 1. Structure of the perturbation, in the plane y = 0 km, shown in entropy contrast δS = S − S . At a time t = 13.22 t adv , the neutrino-driven convection is fully turbulent in the gain layer located between −H ≤ z ≤ H. Outside the two horizontal solid lines, gravity and heating are turned off, while they are at their full intensity inside the horizontal dashed lines. The flow above the gain layer has a uniform velocity v up . The injected perturbation is characterized by a height H p , a vertical spatial wavelength λ z , and horizontal and vertical wavenumbers respectively labeled k x and k z . The angle α corresponds to the tilt angle from the horizontal axis to the direction of the entropy bands. and k z are the horizontal and vertical components of the wavenumber vector ì k, and ω is the temporal frequency. The horizontal wavenumber is computed as where L x = 300 km is the horizontal extent of the gain region. The integer n x corresponds to the number of wavelengths in the x-direction, which plays a role similar to the angular wavenumber of the spherical harmonics representation. The vertical wavenumber is given by where α corresponds to the tilt angle from the x-axis to the direction of the entropy/density bands as shown in Fig. 1. The vertical spatial wavelength can be obtained by Besides, the temporal frequency ω is given by with v up being the vertical velocity of the upstream uniform flow. The density perturbation can thus be rewritten as 3 The perturbation is injected through the upper boundary for a finite duration of time defined as where H p is the vertical extent of the perturbation (Fig. 1). In the following, H p is normalized by the vertical spatial wavelength λ z . In our study, the free parameters of the injected perturbations are chosen to be A, n x , α, and H p . In the rest of the paper, we will consider the following parameter space: A = {5%; 10%; 20%}, α = {30°; 45°; 60°; 75°; 90°}, n x = {1; 2; 4}, and 0.25 ≤ H p /λ z ≤ 4, or correspondingly 1/32 ≤ ∆T inj /t adv ≤ 1. This parameter space covers a wide range of perturbation properties, including those that where observed in numerical simulations Müller et al. 2016). This will allow us to carefully study the impact of each of these parameters on the evolution of the system.

Numerical simulations
The numerical setup employed in this study is similar to the one of Paper I. The RAMSES code (Teyssier 2002;Fromang et al. 2006) is used to simulate the flow dynamics on a Cartesian grid. At the top and bottom boundaries, we use a constant inflow and outflow conditions set by the stationary flow. We use periodic boundary conditions at the lateral edges. The domain is respectively decomposed into (N x × N z ) = (400 × 1200) and (N x × N y × N z ) = (400 × 400 × 1200) uniform cells in 2D and in 3D simulations. Except in Section 5, which focuses on the impact of dimensionality, all of our simulations are performed in 3D.
In Table 1, we present the full list of the simulations discussed in this work. We evolve three distinct sequences of models where we vary one or two perturbation parameters at a time. The first set of models is used to examine the influence of the perturbation height and the simulations are named according to the following convention. In the models labeled HpX, the height of the perturbation is set to H p = Xλ z with X ranging from 0.25 to 8. The corresponding values of injection timescales range from 1/32 t adv to t adv . The other parameters are set to A = 10%, n x = 2, and α = 30°.
Since the impact of the perturbation horizontal wavenumber and amplitude are expected to be of the same order (Müller & Janka 2015;Müller et al. 2017), in our second sequence of models we vary both A and n x . These simulations are named AX nxY , where X corresponds to the perturbation amplitude, ranging from 5% to 20%, and Y to the horizontal wavenumber, ranging from 1 to 4. The tilt angle is set to α = 30°and the perturbation height to H p = 4λ z .
Finally, the impact of the tilt angle α on the flow dynamics is investigated using the simulations labeled as angX where X denotes the value of the angle, ranging from 30°to Table 1. Overview of the simulations discussed in this work. 90°. In these simulations, the remaining parameters are set to A = 10%, n x = 2, and H p = 4λ z .

Non-linear neutrino-driven convection
In the absence of accreting perturbations, the kinetic energy of neutrino-driven convection results from the neutrino heating. The efficiency factors for the conversion of the latter to the former can be expressed in terms of the mass M g , the net neutrino heating L, and the height 2H of the gain layer (Müller & Janka 2015): and respectively for the horizontal and vertical components of the turbulent kinetic energy. The horizontal turbulent kinetic energy is computed as with v x and v y being the velocity components in the x and y-directions, and V gain being the volume of the gain region. The vertical turbulent kinetic energy is similarly computed as The total heating and mass in the gain layer are barely affected by the convective instability. As a consequence, the efficiency factors are directly proportional to the components of the turbulent kinetic energy. Our simulation of non-perturbed neutrino-driven convection shows that both efficiency factors reach asymptotic values in the non-linear stage of the instability, defined as 10 t adv ≤ t ≤ 18 t adv (Fig. 2). We find that η conv, xy ≈ η conv,z ≈ 0.3, which is consistent with the values obtained in realistic 3D simulations of CCSNe (Müller et al. 2017). Figure 3 illustrates the dynamics of gain region when an injected perturbation accretes through this layer for representative model A10 nx2 (see Table 1). As the perturbation enters the gain region, it shrinks in the vertical direction due to the deceleration of the flow. The compression factor equals the velocity reduction ratio (Eq. 9), which amounts to v up /v z,g ≈ 7.5 in our setup, where v z,g is the average vertical velocity in the gain layer. Once an injected perturbation enters the gain layer, at t = 13.2 t adv , it adjusts itself to the new pressure of the heating region (top left panel). The buoyant forces push the higher-entropy parts of the perturbation above the layer of turbulent convection. The lateral expansion of these buoyant entropy bands generates pinched downflows of cold material which then undergo turbulent convection. As the perturbation accretion continues, the higher-entropy material gradually fills the whole upper part of the gain region (top right panel). This process saturates when the first full vertical wavelength of the perturbation has entered the gain layer. The base of the perturbation is then subject to convective mixing that disrupts the entropy bands while the rest of the perturbation experiences the effects of buoyancy and pressure adjustment (bottom left panel). The accumulation of the buoyant material in the upper part of the gain layer tend to lower the altitude at which the perturbation is disrupted by turbulent convection (bottom right panel).

Perturbation height
In order to study the impact of the perturbation height on the turbulence in the gain layer, we focus on the models with varying H p , labeled HpX in Table 1. We calculate the squared turbulent Mach number in the i-direction with c being the sound speed and v x = v y = 0 since the flow is initially purely vertical. Fig. 4 shows the components of Ma 2 in the x and z-direction. The components of Ma 2 begin to increase as the perturbation starts to accrete through the gain region. After a phase of identical increase in all cases, the growth becomes dependent on H p . The maximum of Ma 2 increases with the perturbation height up to H p = 2λ z , which corresponds to an injection time T inj = t adv /4. For larger values of H p , the maxima do not grow with H p anymore. For such perturbations, the growth of Ma 2 stops before the end of the perturbation accretion through the upper edge of the gain layer. This is likely a consequence of the accumulation of buoyant material in the upper part of the gain layer, as discussed above. This phenomenon tends to reduce the size of the region of neutrino-driven convection and thus to weaken turbulence. After reaching their maxima, the components of Ma 2 remain higher than those of the non-perturbed convection for several advection times. During this post-maximum phase, we observe secondary peaks in Ma 2 in some perturbed models and the non-perturbed case. These peaks are a result of the stochastic nature of the neutrino-driven convection in the gain layer. In a more realistic framework, the additional turbulent pressure generated by the accretion of extended perturbations would move the shock outwards. These perturbations may then be entirely disrupted inside the enlarged gain layer before quenching convection. Figure 5 shows the impact of the perturbation amplitude A and its horizontal wavenumber n x on the neutrino-driven convection. The horizontal and the vertical components of the turbulent kinetic energy in the gain layer are normalized by the averaged components measured in the reference simulation during the non-linear regime of convection, i.e. 10 ≤ t/t adv ≤ 18. In all perturbed models, the turbulent kinetic energy reaches an initial peak as the perturbation accretes through the gain region. The early growth of the turbulent kinetic energy is stronger for perturbations with larger amplitude because larger gravitational potential energy is converted into kinetic energy. Similarly, the turbulent kinetic energy grows faster for lower horizontal wavenumbers. This is a consequence of a larger vertical spatial wavelength over which the mixing occurs (Eq. 8).

Perturbation amplitude and horizontal wavenumber
During the accretion phase of the perturbation, doubling A leads to a similar level of turbulent kinetic energy as  Table 1) after 13.48, 13.7, 13.92, and 14.2 advection timescales (left to right, top to bottom). The horizontal lines are defined in Fig. 1. halving n x . For instance, the three perturbations: A5 nx1, A10 nx2, and A20 nx4 produce roughly the same amount of turbulent kinetic energy as shown in Fig. 5. This suggests that the impact of n x on the convective instability is as large as that of A. This implies that the level of turbulence reached in the gain layer is correlated with the ratio A/n x . We defer a more detailed analysis of this conjecture to Section 4. The correlation between the ratio A/n x and the components of the turbulent kinetic energy holds for both the horizontal and the vertical directions (Fig. 5). The turbulent kinetic energy is initially higher in the horizontal directions. As we will see in the next section, this behavior is specific to the particular value of the tilt angle α used here. When the accretion of the perturbation stops, the two components are found to be roughly equal to each other, similarly to non-perturbed convection (Paper I).
Once the turbulent kinetic energy reaches its maximum, we observe that all components remain larger than their counterparts in the simulation of non-perturbed con-vection. In this post-accretion phase, the turbulent kinetic energy is larger for models with larger perturbation amplitudes, especially in the vertical direction (Fig. 5, right  panel). Compared to the reference model, the vertical turbulent kinetic energy is about 30% and 50% higher for A = 10% and A = 20% respectively. The parameter n x seems to play a minor role on the post-peak phase. For example, the turbulent kinetic energy is about 10% larger in simulations with n x = 1 compared to those with n x = 4.

Tilt angle
In this section, we investigate the influence of the tilt angle α of the perturbation using the set of simulations with varying α, labeled angX in Table 1. The other parameters are kept fixed: A = 10% and n x = 2. The time evolution of E kin, xy and E kin,z for these cases are shown in Fig. 6.
The vertical component E kin,z increases with α (right panel of Fig. 6). For example, the peak value of E kin,z for  Table 1). The leftmost vertical dashed line marks the arrival of the injected perturbation inside the gain layer (in black). The other lines depict the end of the perturbation accretion through the upper edge of the gain layer. For readability, the data are smoothed out over a convective turnover timescale. α = 90°is twice larger than that for α = 30°. This behavior can be explained by the structure of the perturbations. As α approaches 90°, the entropy bands become more vertical and thus less prone to turbulent mixing. Instead, upon entering the gain region, the higher entropy regions inflate horizontally, generating narrow vertical downflows of colder material channeled to the bottom part of the gain layer, as shown in Fig. 7. These downflows are sustained throughout the accretion of the perturbation, generating vertical motion.
In contrast, the turbulent kinetic energy in the horizontal directions does not grow monotonically with α (left panel of Fig. 6). The maximum of E kin, xy is reached for α = 60°a nd is about 50% higher than that reached for α = 90°. This behavior can be explained by comparing the temporal frequency ω of the perturbation, defined in Eq. 9, with the convective turnover frequency ω conv in the gain layer. Following Müller et al. (2016), we compute the average frequency ω conv during the non-linear phase of the convective instability as where v conv corresponds to the velocity of the convective turnovers. In this definition, we consider only the central −H/2 ≤ z ≤ H/2 part of the gain region, because the convective turnovers are mostly confined to this region. Fig. 8 shows that perturbation frequency ω matches the convective turnover frequency ω conv at α crit ≈ 53°. This is close to the value α = 60°that yields the largest E kin, xy among our models (Fig. 6). This suggests that the coupling between the perturbation and the convection reaches a maximum efficiency for ω ≈ ω conv . On the contrary, the peak of E kin, xy decreases with increasing |α − α crit |.
After the maximum of turbulent kinetic energy is reached, the horizontal and vertical components of the turbulent kinetic energy become roughly equal. There is no obvious correlation between α and the amount of turbulent kinetic energy at the late times.

Qualitative model
Upon entering the gain region, the injected perturbations amplify the turbulent kinetic that region. In a quest to explain the mechanism for this amplification, Müller et al. (2016) considered several candidates, including the flux of acoustic energy (F ac ), the conversion of gravitational potential energy into kinetic energy (F pot ), and the flux of transverse kinetic energy (F t ). Based on an order-of-magnitude estimate, they concluded that the additional turbulence generated by accreting perturbations mainly results from the work done by buoyancy. The simplicity of our setup allows us to verify this estimate. In our setup, F ac is given by where L x = L y = 300 km are the horizontal dimensions of the gain layer. δP and δv represent the pressure and velocity fluctuations with respect to horizontally averaged quantities. The term F t is computed as where δv h refers to the fluctuations of horizontal velocity. M corresponds to the mass accretion rate, defined as Finally, F pot is given by ∆Φ is defined as the difference of gravitational potential between the upper and the lower edges of the gain layer. We evaluate the three driving terms at the upper edge of the gain layer using our simulations performed with A = 10% and with n x ranging from 1 to 4. Fig. 9 shows that F pot exceeds F ac and F t by at least an order of magnitude, confirming the hypothesis of Müller et al. (2016). Besides, F pot becomes the dominant energy input over neutrino heating during the accretion of the perturbation.
In the absence of infalling perturbations, the turbulent kinetic energy in the gain layer can be obtained from Eqs. 12-13: and where the average efficiency factors arē η conv, xy ≈η conv,z ≈ 0.3 in the non-linear phase of convection (cf. Fig. 2). The increase of turbulent kinetic energy  Table 1 for the rest of the parameters. All quantities are divided by the average values computed in the non-linear regime of neutrino-driven convection (black). For readability, the data are smoothed out over a turnover timescale.  Table 1 for the rest of the parameters.  (Table 1). At the time t = 14.11 t adv , the injected perturbation has been partially accreted into the gain layer. The initially vertical entropy bands are deformed due to buoyancy force and adjustment to the new pressure. As a consequence, the bands of higher-entropy material (in red) expand laterally and pinch the bands of lower-entropy material (in blue) into narrow downflows reaching deeper regions of the gain layer than the rest of the perturbation. The horizontal lines are defined as in Fig. 1. Figure 8. Temporal frequency ω of the perturbation, normalized by the convective turnover frequency ω conv , as a function of the tilt angle α (solid curve). The horizontal dotted line depicts the equality between these two frequencies which is reached for α crit ≈ 53°.
due to injection of perturbations can be accounted for via a correction term (Müller et al. , 2017: and where The parameter Λ p measures the dissipation length of turbulence inside the injected perturbation.It corresponds to the driving scale of turbulent mixing which amounts to average vertical wavelength of the perturbation in the gain layer, i.e. λ z v z,g /v up (see Section 3.2). This is thus different from the regime of neutrino-driven convection where the dissipation length is the height of the gain layer, i.e. 2H. Note that Λ p would correspond to the size of the infalling convective eddies in spherical geometry, i.e. πr sh / .
Since the vertical wavelength of the perturbation grows with α, it becomes larger than the gain layer height above a certain α (Eq. 8). However, the dissipation length cannot be larger than the gain region, so we limit the maximum value of Λ p to H. This choice is justified for two reasons. First, Müller et al. (2016) advocated that the dissipation length associated to convection is limited to half of the gain layer height. Second, our simulations with large values of α show that turbulent mixing occurs only in the lower half of the gain layer, as it can be seen in Fig. 7. The upper half is dominated by large bands of high-entropy material and narrow downflows without signs of turbulent mixing. Hence, we define the dissipation length as The threshold α, above which Λ p = H, is

Comparison to simulations
In order to assess the robustness of the estimates introduced in the previous section, we compute the average values of the right-hand side of Eqs. 24-25 during the accretion phase of the perturbation, i.e. from t = 13.2 t adv to t = 14.2 t adv , and compare these terms to the average quantities E kin, xy /M g and E kin, xy /M g measured in our simulations. In a more realistic setup, the additional turbulent pressure would result in wider gain region and larger neutrino heating, which cannot be reproduced in our stationary flow. Thus, we do not attempt to model the turbulent kinetic energy during the post-accretion phase of the perturbation. The comparison for the set of models with varying A and n x , labeled AX nxY in Table 1, is displayed in Fig. 10. Note that Λ p < H in all these models. We obtain agreement within ∼ 10% and ∼ 20% respectively for the horizontal (left panels) and the vertical (right panels) directions. In all but one case, the estimates are higher than the values measured in the simulations. A shorter and more elaborate dissipation length than Λ p may improve the agreement even further. Figure 10 also confirms that the additional turbulent kinetic energy is mostly set by the ratio A/n x . For a given ratio, the values measured in the corresponding simulations are remarkably close to each other. The only exceptions are the models with A = 20%. In these cases, large buoyant bubbles compress the region of neutrino-driven convection which reduces the turbulent kinetic energy compared to simulations with smaller perturbation amplitudes (see Section 3.2 for a description of this phenomenon). This suggests that A = 20% could be an upper bound on the perturbation amplitudes that can be accurately studied with our idealized model.
We now consider the sequence of models with varying α, labeled angX in Table 1, for which Λ p = H when α ≥ α thresh ≈ 68°. The results of the comparison are shown in Fig. 11. Overall, the agreement is similar to the previous comparison described above. They amount to about 10% and 20% respectively for the horizontal (left panels) and the vertical (right panels) directions. The agreement is less satisfactory for E kin, xy at large angles. In these cases, the estimate becomes almost twice larger than the values measured in the simulations.
In the horizontal direction (Fig. 11, left panels), the turbulent kinetic energy is rather well approximated when α < α thresh ≈ 68°, that is when Λ p = λ z v up /v z,g . The estimate is well suited to capture the trend of an increasing horizontal turbulent kinetic energy with the angle, up to the maximum located close to α crit ≈ 53°. The gap between the simulations and the estimate becomes much wider when α > α thresh and Λ p = H. The estimate results in a constant turbulent kinetic energy, independent from the angle, whereas E kin,xy decreases from 60°to 90°in our simulations. This decrease may be a consequence of a large ω that exceeds the convective turnover frequency. The temporal modulation of the perturbation, at large α, may become too slow to affect neutrino-driven convection. This suggests that a more elaborate estimate could be designed by considering the distance between α and α crit in the definition of Λ p (Eq. 27). Our study shows that for the horizontal direction, the amount of turbulent kinetic energy is maximized when the convective turnover frequency matches ω. Figure 10. Upper panels: The average ratios E kin, x y /M g and E kin, z /M g measured in the simulations with varying A and n x (models AX nxY in Table 1) during the perturbation accretion (open circles) are compared to the average estimates (filled circles) derived from Eqs. 24-25 for the horizontal direction (left panel) and the vertical one (right panel). These data are plotted as a function of the ratio A/n x . Different perturbation amplitudes are represented by different colours (see lower legend). The horizontal lines represent the average ratios E kin, x y /M g and E kin, z /M g measured during the non-linear phase of neutrino-driven convection. These two average values are multiplied by the correction factor (1 +ψ) 2/3 (black triangles) to indicate the expected increase of turbulent kinetic energy based solely on the input parameters. Lower panels: The ratios of average turbulent kinetic energy measured in the simulations to the estimates are shown as a function of the ratio A/n x , respectively for the horizontal direction (left panel) and the vertical one (right panel). Figure 11. Upper panels: The average ratios E kin, x y /M g and E kin, z /M g measured in the simulations with varying α (models angX in Table 1 during the perturbation accretion (open circles) are compared to the average estimates (filled circles) derived from Eqs. 24-25 for the horizontal direction (left panel) and the vertical one (right panel). These data are plotted as a function of the tilt angle α. The horizontal lines represent the average ratios E kin, x y /M g and E kin, z /M g measured during the non-linear phase of neutrino-driven convection. These two average values are multiplied by the correction factor (1 + ψ) 2/3 (black triangles) to indicate the expected increase of turbulent kinetic energy based solely on the input parameters. Lower panels: The ratios of average turbulent kinetic energy measured in the simulations to the estimates are shown as a function of the angle α, respectively for the horizontal direction (left panel) and the vertical one (right panel).
In the vertical direction, both the estimate of the turbulent kinetic energy and the values obtained in the simulations increase with larger α. The estimate is about 20% higher than the results of the simulations throughout the range of angles considered. E kin,z is an increasing function of α up to an asymptotic value reached for α ∼ α thresh . In this regime of large angles, the size of the mixing region is not set by the vertical wavelength but is instead limited to half of the gain layer height. The maximum turbulent kinetic energy in the vertical direction is thus restricted by the size of the gain layer.
Alternatively, E kin, xy /M g and E kin,xy /M g can be measured by just upscaling their values in the simulation of non-perturbed convection by the correction term (1 + ψ) 2/3 , where ψ is computed using the input parameters of the stationary flow and the perturbation (black triangles in Figs. 10  and 11). This approximation can be understood as the maximum turbulent kinetic energy expected for a given accreting perturbation. The latter estimate is significantly larger than the one from Eqs. 24-25. In fact, this alternate estimate neglects the role of the acoustic feedback of neutrino-driven convection which alters the perturbation before it enters the gain layer. As a consequence, the magnitude of the driving term F pot is reduced compared to the theoretical value expected from the stationary flow and the perturbation parameters. This reduction can be witnessed in Fig. 9 where F pot is not constant during the perturbation accretion but is decreasing. Upscaling the turbulent kinetic energy of the unperturbed convective flow leads to an overestimate by at least a factor 2 in all cases. This alternate method is thus much less accurate than the approach described above.

Turbulent spectra
In order to characterize the change of spatial scales of the convective motions during the accretion phase of the perturbation, we examine the turbulent kinetic energy spectra of different simulations. The turbulent kinetic energy density is first decomposed into Fourier coefficientŝ over each horizontal plane Ω h , contained in the region −H ≤ z ≤ H and then averaged over this height. v turb is the turbulent velocity, defined as To reduce the stochastic character of the turbulent spectra, we also average the Fourier coefficient over a time range [t, t + t conv ]. For a given wavenumber k, the associated component of the power spectrum is then defined as First, we focus on the turbulent spectra for the models AX nxY that have varying A and n x (cf. Table 1). Fig. 12 shows the turbulent spectra taken at the time when the maximum turbulent kinetic energy is reached. The spectra of perturbed cases exceed that of the non-perturbed model. The only outlier is the simulation with A = 5% and n x = 4, which exhibits a negligible increase of turbulent kinetic energy (Fig. 5). The scales at which the spectra E(k) reach Figure 12. Spectra of turbulent kinetic energy computed in the models with varying A and n x (models AX nxY in Table 1). The spectra are computed at the time when maximum turbulent kinetic energy is reached. The spectrum obtained in the simulation of non-perturbed convection is shown with a black dotted curve.
their peaks are closely related to the structure of the perturbations, at least when A ≥ 10%. The spectra of the models with n x = 2 show a peak at k = 2, while the ones with n x = 4 display a peak at k = 5. On the contrary, the spectra of all models with A = 5% peak at k = 1. These spectra resemble that of the non-perturbed convection which also peaks at k = 1. This could be due to a faster disruption of such perturbations which rapidly leads to a system dominated by the spatial scales of the (unperturbed) neutrino-driven convection.
Interestingly, spectra for models with a given ratio A/n x seem to overlap at intermediate scales, e.g. between k = 2 and k = 10. This is consistent with the observation that the turbulent kinetic energy is mostly set by the ratio A/n x (see Section 3.3). However, these overlaps are not systematic and not always clearly visible. For example, the case A10 nx2 seems closer to A20 nx2 instead of A20 nx4. That said, these overlaps involve spatial scales that are sub-dominant in the turbulent kinetic energy budget.
In Paper I, it was shown that the spectrum of unperturbed 3D neutrino-driven turbulent convection scales as ∝ k −5/3 for 10 k 30 (see also Radice et al. 2016;Nagakura et al. 2019). In order to identify the inertial range of turbulence in perturbed models, we show the compensated spectra E(k)k 5/3 in Fig. 13. We find that an inertial range, where the function E(k)k 5/3 is flat, can clearly be observed only for n x = 1. For larger n x , the inertial range is smaller. The perturbations with n x = 1 reproduce more closely the properties of the turbulence generated by (unperturbed) neutrino-driven convection, which favours the largest spatial scales (Kazeroni et al. 2018). Figure 14 shows the turbulent spectra of models with varying α, labeled angX in Table 1. All perturbed cases display a peak in their spectra at k = 2, which corresponds to the input parameter n x = 2 (top panel). The peak is maximized for α = 60°, which is consistent with the angledependence of the turbulent kinetic energy that we observed earlier in Fig. 6. The comparison of these spectra shows that the total turbulent kinetic energy is again maximized for angles such that ω ≈ ω conv .
The spectra of the models with large α are characterized by a pronounced "odd-even" deformation. As the angle approaches 90°, the bands of entropy become vertical and this Figure 13. Spectra of turbulent kinetic energy compensated by k 5/3 for models with varying A and n x (models AX nxY in Table 1). The simulations performed respectively with n x = 1, n x = 2, and n x = 4 are compared to the reference simulation (black dotted curve) in the upper, middle, and lower panels.
prevents efficient disruption of the perturbation, as it was illustrated in Fig. 7. Instead, the perturbation is dominated by the initial spatial scale n x = 2 throughout its accretion into the gain layer. We hypothesize that this behavior is the reason for the "odd-even" shape of the spectra. The large angles not only enhance the modulations of the spectra, but also reduce the range of wavenumbers for which an inertial range could be identified (cf. the bottom panel of Fig. 14). Such a range is hardly noticeable for α ≥ 60°.

IMPACT OF THE DIMENSIONALITY
In order to study the influence of the dimensionality on our results, we simulate in 2D both the unperturbed model and the sequence of models with varying A and n x (labeled AX nxY in Table 1). Fig. 15 shows the ratio of the components of the turbulent kinetic energy between 2D and 3D. For the non-perturbed model, the turbulent kinetic energy is twice larger in 3D than in 2D. This is different from the results of Paper I where a stationary flow with χ = 5 was considered and the turbulent kinetic energy was found to be slightly larger in 2D. In our current work, we consider a flow with χ = 3 which is closer to the marginal stability.  Table 1). The spectra are computed at the time when maximum turbulent kinetic energy is reached. The spectrum obtained in the simulation of nonperturbed convection is shown with a black dotted curve. Lower panel: Same as above but for spectra compensated by k 5/3 . For such a marginally unstable flow, convection is triggered more rapidly in 3D and buoyant bubbles reach higher altitudes and velocities (Kazeroni et al. 2018).
The two-dimensional nature of the flow has a profound impact on the turbulence generated by an accreting perturbation. The turbulent kinetic energy measured in 2D simulations is between 2 to 5 times larger than that in 3D runs, as shown in Fig. 15. The gap between 2D and 3D increases the perturbation amplitude. On the other hand, the horizontal wavenumber does not significantly affect the gap. In 2D, the disruption of the injected perturbation leads to the formation of large scale vortices that tend to merge until reaching the gain layer height. Once formed, such a flow pattern may remain in the gain layer for several advection times before being expelled in the downstream region (Kazeroni et al. 2018). For this reason, large peaks of turbulent kinetic energy can be observed long after the perturbation has fully accreted into the gain region (Fig. 15).
Contrary to the 3D simulations, the turbulent kinetic energy in 2D simulations does not correlate with the ratio A/n x (Fig. 16, left panels). For a given A/n x , the mean values of E kin, x /M g and E kin,z /M g are spread over a wide range that overlaps with the distribution of values obtained for other values of A/n x . Both E kin, x /M g and E kin,z /M g increase with lowering n x (Fig. 16, right panels). At the same time, there is no clear correlation between the perturbation amplitude A and the ratios E kin, x /M g and E kin,z /M g . This is related to the inverse turbulent cascade that takes place in 2D. When a perturbation with n x = 1 accretes into the gain region, its geometry already corresponds to that of the largest possible scale of the convection. This results in a more rapid growth  Table 1 for the rest of the parameters. For readability, the data are smoothed out over a turnover timescale. Figure 16. The average ratios E kin, x /M gain (upper panels) and E kin, z /M gain (lower panels) measured in the 2D simulations with varying A and n x are plotted against the ratio A/n x (left panels) and the perturbation amplitude A (right panels). See the models AX nxY in Table 1 for the rest of the parameters. of the turbulence across the gain region. For this reason, the amplitude is less important than n x .

CONCLUSIONS AND DISCUSSION
We have studied the coupling between injected density perturbations and neutrino-driven convection in the gain layer of a CCSN, employing parameterized 3D simulations of an idealized model. Our study considered a stationary flow where neutrino-driven convection has reached the non-linear stage (Kazeroni et al. 2018) and into which density perturbations are injected. The latter originates from the interaction between pre-collapse inhomogeneities and the stalled shock of a CCSN. The amplitude, the dominant spatial scale, and the temporal frequency of the injected density perturbations were varied in a systematic manner to characterize the additional turbulence generated by the accretion of these fluctuations. A comparison between 2D and 3D simulations was conducted to clarify the impact of the dimensionality on the additional turbulence. Our main findings can be summarized as follows.
1. The turbulent kinetic energy in the gain layer scales as A/n x , where A is perturbation amplitude and n x is the horizontal wavenumber. For example, doubling A or halving n x result in a similar increase of the turbulent kinetic energy. As a consequence, the impact of small-scale perturbations is weaker. This supports the hypothesis that only asymmetries at the largest spatial scales, e.g. = 1 − 2, could be beneficial to revive the shock wave (e.g. Müller & Janka 2015;Müller et al. 2016).
2. The influence of the temporal frequency of the perturbation on neutrino-driven convection was investigated by varying the tilt angle α between the horizontal axis and the direction of the entropy bands of the perturbation. In the horizontal direction, the turbulent kinetic energy is maximized when the frequency of the perturbation is close to that of the convective turnovers. The vertical component increases monotonically with α. When α ≈ 90°, the vertical entropy bands become less prone to turbulent mixing and their disruption leads to greater vertical displacements of matter which increases the turbulent kinetic energy in that direction. Our work suggests that density perturbations may non-linearly increase the turbulent kinetic energy of the convective instability by a few tens of percents for a wide range of perturbation frequencies.
3. We assessed the accuracy of the qualitative model of the turbulent kinetic energy generated by an accreting perturbation by Müller et al. (2016Müller et al. ( , 2017). An agreement within 20% was obtained between the simulations and the estimate for most of our parameter space. The only exception arises at large tilt angles, for which low perturbation frequencies barely affect neutrino-driven convection. Our study unambiguously demonstrates that the additional turbulent kinetic energy originates from the work of buoyancy onto the density perturbations, at least when the neutrino-driven convection dominates the post-shock dynamics, in agreement with Müller et al. (2016Müller et al. ( , 2017. 4. The turbulent kinetic energy in the gain layer increases with the scale of the perturbation both in 2D and in 3D. However, 2D simulations overestimate the additional turbulent kinetic energy by a factor of ∼ 2 − 5 compared to their 3D counterparts. This difference stems from the direction of the turbulent energy cascade (Murphy et al. 2013;Kazeroni et al. 2018). Unlike 3D models, the fluctuations are not disrupted into smaller scales in 2D. This leads to a rapid growth of large-scale turbulent motions, resulting in unphysically large turbulent kinetic energies (cf. Section 5 for more details).
Our study represents an important step towards a comprehensive description of the interaction of pre-collapse fluctuations with non-linear instabilities in the gain layer of CCSNe. We have clarified the origin and the properties of the additional turbulence related to the coupling between density perturbations and neutrino-driven convection. Our results are in line with the few 3D CCSN simulations performed from the final minutes of the progenitor before collapse to the shock revival. Simplified models remain essential to analyze in depth the turbulence resulting from the interplay between different types of pre-collapse fluctuations and the variety of post-shock instabilities while covering a broad parameter space.
In our study, we neglected the role played by the stationary shock. The inclusion of the latter could lead to a subtle interplay between SASI and neutrino-driven convection as witnessed by Müller et al. (2017). In that study, the accreting perturbation, coupled to the shock oscillations, excited turbulent motions primarily in the radial direction whereas turbulence was enhanced equivalently in all directions for our neutrino-driven convection dominated flow. Adding a shock wave to our model could reveal how neutrino-driven convection and SASI are enhanced or damped by an accreting perturbation. Our setup including a stationary shock would also be insightful to quantify the reduction of the critical luminosity for any type of accreting perturbation interacting with a given post-shock instability. This will be the subject of the future studies.