Response of nuclear-dissociating shocks to vorticity perturbations

TL;DR

This paper analyzes how nuclear dissociation influences the response of standing accretion shocks to vorticity perturbations in core-collapse supernovae, revealing stabilization effects and increased downstream density perturbations.

Contribution

It provides a linear theoretical analysis of shock-vorticity interactions considering nuclear dissociation effects, focusing on non-isotropic perturbations relevant to collapsing stellar shells.

Findings

Nuclear dissociation stabilizes shock oscillations.

Dissociation increases downstream density perturbations.

Exact distributions of rotational and acoustic perturbations are derived.

Abstract

In the context of core-collapse supernova explosions (CCSNe), the interaction of standing accretion shocks with upstream vorticity perturbations is investigated by linear theory analysis. The endothermic effect associated to the nuclear dissociation, which takes place right behind the shock wave, affects the amplitude of the perturbations amplified/generated across the front. For upstream disturbances whose characteristic size is much larger than the post-shock dissociation-layer thickness, the effect of nuclear dissociation can be reduced to that of considering the global endothermic effect that scales with the inflow energy flux. The present study focuses on perturbation fields that are not isotropic, which mimic the perturbations in collapsing convective shells of massive stars. The linear interaction of the shock with bidimensional mono-frequency vorticity perturbations is theoretically addressed, with the limit of highly-stretched vortices being analyzed in detail. The exact spatial distribution of the rotational and acoustic perturbations generated in the post-shock flow are provided along with the transient evolution of the shock front. It is found that nuclear dissociation contributes to stabilize the shock oscillations, but increases the amplitude of the density perturbations downstream. An extension of this work that addresses the interaction with tridimensional isotropic turbulent flows can be found in reference Huete, C., et al. 2018, MNRAS, 475, 3305-3323, which analyzes the effect of the post-shock flow on the critical conditions that ultimately trigger explosion.

Figures

Equations

$${\rho }_1\left(u_1^{\prime }+{\dot{R}}_{\mathrm{shock}}\right)={\rho }_2\left(u_2^{\prime }+{\dot{R}}_{\mathrm{shock}}\right),$$

$$p_1+{\rho }_1{\left(u_1^{\prime }+{\dot{R}}_{\mathrm{shock}}\right)}^2=p_2+{\rho }_2{\left(u_2^{\prime }+{\dot{R}}_{\mathrm{shock}}\right)}^2,$$

$$e_1+\frac{p_1}{{\rho }_1}+\frac{1}{2}{\left(u_1^{\prime }+{\dot{R}}_{\mathrm{shock}}\right)}^2=e_2+\frac{p_2}{{\rho }_2}+\frac{1}{2}{\left(u_2^{\prime }+{\dot{R}}_{\mathrm{shock}}\right)}^2.$$

$$e_1-e_2=\frac{1}{\gamma -1}\frac{p_1}{{\rho }_1}-\frac{1}{\gamma -1}\frac{p_2}{{\rho }_2}+\Delta e_{\mathrm{dis}},$$

$$\Delta e_{\mathrm{dis}}=\epsilon \frac{a_1^2}{\gamma -1}\left(1+\frac{\gamma -1}{2}{M}_1^2\right),$$

$$C_2=\frac{{\rho }_2}{{\rho }_1}=\frac{u_1}{u_2}=\frac{\left(\gamma +1\right)M_1^2}{\left(\gamma -\kappa \right)M_1^2+1},$$

$$P_2=\frac{p_2}{{\rho }_1{u}_1^2}=\frac{\gamma M_1^2(1+\kappa )+1}{\gamma (\gamma +1)M_1^2},$$

$$M_2=\frac{u_2}{a_2}={\left(\gamma C_2{P}_2\right)}^{-1/2}={\left[\frac{\left(\gamma -\kappa \right)M_1^2+1}{\gamma M_1^2(1+\kappa )+1}\right]}^{1/2},$$

$$\kappa ={\left[(1-M_1^{-2}{)}^2+\epsilon (\gamma +1)\left(\gamma -1+2M_1^{-2}\right)\right]}^{1/2}$$

$${\bar{u}}_1\left(x_1,y_1\right)=\frac{u_1-⟨u_1⟩}{{\hat{u}}_1⟨a_2⟩}=\cos \left(k_x{x}_1\right)\cos \left(k_y{y}_1\right),$$

$${\bar{v}}_1\left(x_1,y_1\right)=\frac{v_1-⟨v_1⟩}{{\hat{u}}_1⟨a_2⟩}=\frac{k_x}{k_y}\sin \left(k_x{x}_1\right)\sin \left(k_y{y}_1\right),$$

$${\bar{\omega }}_1\left(x_1,y_1\right)=\frac{\partial {\bar{v}}_1}{\partial (k_y{x}_1)}-\frac{\partial {\bar{u}}_1}{\partial (k_y{y}_1)}=\left(1+\frac{k_x^2}{k_y^2}\right)\cos \left(k_x{x}_1\right)\sin \left(k_y{y}_1\right),$$

$$\bar{p}=\frac{p-⟨p_2⟩}{{\hat{u}}_1\gamma ⟨p_2⟩},\bar{\rho }=\frac{\rho -⟨{\rho }_2⟩}{{\hat{u}}_1⟨{\rho }_2⟩},\bar{u}=\frac{u-⟨u_2⟩}{{\hat{u}}_1⟨a_2⟩},\bar{v}=\frac{v-⟨v_2⟩}{{\hat{u}}_1⟨a_2⟩}.$$

$$\frac{\partial \bar{\rho }}{\partial \tau }+\frac{\partial \bar{u}}{\partial x}+\bar{v}=0,\frac{\partial \bar{u}}{\partial \tau }+\frac{\partial \bar{p}}{\partial x}=0,\frac{\partial \bar{v}}{\partial \tau }-\bar{p}=0,\frac{\partial \bar{p}}{\partial \tau }=\frac{\partial \bar{\rho }}{\partial \tau },$$

$$\frac{{\partial }^2\bar{p}}{\partial {\tau }^2}=\frac{{\partial }^2\bar{p}}{\partial x^2}-\bar{p}.$$

$$\left(C_2-1\right){\dot{\xi }}_s=C_2{\bar{u}}_s-M_2{C}_2{\bar{\rho }}_s-{\bar{u}}_1,\frac{}{}$$

$${\bar{p}}_s=2M_2\left({\bar{u}}_s-{\bar{u}}_1\right)-M_2^2{\bar{\rho }}_s,\frac{}{}$$

$$M_1^2{M}_2^2{\bar{\rho }}_s={\Pi }_s{\bar{p}}_s-{\Delta }_s\left({\dot{\xi }}_s-{\bar{u}}_1\right),\frac{}{}$$

$${\bar{v}}_s=M_2\left(C_2-1\right){\xi }_s+{\bar{v}}_1,$$

$${\Pi }_s=\frac{M_1^2{\left[1+M_1^2\left(1-\kappa \right)\right]}^2}{{\left(M_1^2+1\right)}^2-M_1^4{\kappa }^2}$$

$${\Delta }_s=\epsilon \frac{2M_2{M}_1^4\left(\gamma -1\right)\left[1+M_1^2\left(1-\kappa \right)\right]}{{\left(M_1^2+1\right)}^2-M_1^4{\kappa }^2},$$

$$\frac{d{\xi }_s}{d\tau }={\sigma }_a{\bar{p}}_s+\cos \left(\frac{k_x}{k_y}{C}_2{M}_2\tau \right),$$

$$\left({\sigma }_b+M_2\right)\frac{\partial {\bar{p}}_s}{\partial \tau }$$

$${\sigma }_a=\frac{C_2\left(M_1^2-{\Pi }_s\right)}{2M_2{M}_1^2\left(C_2-1\right)+C_2{\Delta }_s},{\sigma }_b=\frac{M_1^2+{\Pi }_s+{\Delta }_s{\sigma }_a}{2M_2{M}_1^2}$$

$$x=r\sinh \chi ,\tau =r\cosh \chi$$

$$r\frac{{\partial }^2\bar{p}}{\partial r^2}+\frac{\partial \bar{p}}{\partial r}+r\bar{p}=\frac{1}{r}\frac{{\partial }^2\bar{p}}{\partial {\chi }^2}.$$

$$\frac{d{\xi }_s(r)}{dr}=\frac{{\sigma }_a}{\sqrt{1-M_2^2}}{\bar{p}}_s(r)+\frac{1}{\sqrt{1-M_2^2}}\cos \left(\zeta r\right)$$

$${\frac{1}{r}\frac{\partial {\bar{p}}_s}{\partial \chi }|}_s=-{\sigma }_b\frac{\partial {\bar{p}}_s(r)}{\partial r}-\frac{M_2^2\left(C_2-1\right)}{\sqrt{1-M_2^2}}{\xi }_s(r)+\zeta \frac{C_2-1}{C_2}\sin \left(\zeta r\right),$$

$$\zeta =\frac{k_x}{k_y}\frac{M_2{C}_2}{\sqrt{1-M_2^2}}=\frac{{\omega }_s}{\sqrt{1-M_2^2}}=\frac{1}{\tan \theta }\frac{M_2{C}_2}{\sqrt{1-M_2^2}}$$

$$\mathcal{F}(s,\chi )={\int }_0^{\infty }f(r,\chi )e^{-sr}\mathrm{d}\mathrm{r}$$

$${\mathcal{P}}_s(s)=\frac{s\left(1+{\sigma }_b\right){\bar{p}}_{s0}}{s\sqrt{s^2+1}+{\sigma }_b{s}^2+{\sigma }_c}+\frac{s\sigma }{\left(s\sqrt{s^2+1}+{\sigma }_b{s}^2+{\sigma }_c\right)\left(s^2+{\zeta }^2\right)},$$

$$\sigma =\frac{C_2-1}{C_2}\left({\zeta }^2-\frac{M_1^2}{M_1^2-1}\right)$$

$${\sigma }_c=\frac{M_2^2\left(C_2-1\right)}{1-M_2^2}{\sigma }_a.$$

$${\bar{p}}_s(r)=-\frac{2}{\pi }{\int }_0^1\cos (zr)f(z)\mathrm{d}z+\frac{2\sigma }{\pi }{\int }_0^1\frac{\cos (zr)-\cos (\zeta r)}{{\zeta }^2-z^2}f(z)\mathrm{d}z$$

Full text

††: Phys. Scr.

Keywords: Shock Wave, Supernovae, Turbulence

1 Introduction

A distinguishing feature of massive stars is that the nuclear fusion occurring in their cores continues even after the hydrogen fuel is exhausted. The high temperatures induced by the strong gravity allow heavier elements such as helium and carbon to fuse sequentially. The energy released due to the mass defect between products and reactants keep the star stable, extending the lifetime of the star. Nonetheless, the characteristic time of these advanced fusion stages decreases very rapidly with the nuclear mass, with heavier nuclei burning on a timescale orders of magnitude shorter than the hydrogen sequence, which counts in million of years [1]. In stars with initial masses $\sim(8-100)M_{\odot}$, the sequential nuclear fusion lasts until the formation of iron nuclei, a point beyond which nuclear fusion is no longer exothermic. As a result, the iron nuclei accumulates in the center, forming a core supported by the pressure of degenerate electrons. When the core reaches its maximum mass of $\simeq 1.4M_{\odot}$, pressure begins to decay and the hydrodynamical stability breaks down, triggering a collapse of the iron core (e.g., [2] for a recent review).

The core collapse accelerates until the central density becomes as high as nuclear density ($\sim 2\times 10^{14}\,\mathrm{g/cm^{3}}$), a point where nuclear matter stiffens. This abruptly halts the collapse of the inner iron core, leading to the formation of a shock wave at the boundary of the inner core. The shock has to expel the stellar envelope and thus power core-collapse supernova (CCSN) explosion, leaving behind a stable neutron star (NS). The propagation of the shock, however, does not progress smoothly. The inherent pressure and temperature rise across the shock produces heavy-nuclei breaking as it propagates outwards, with associated energy consumption. In addition, the hot material behind the shock cools rapidly due to copious neutrino emission. As a result, the shock quickly loses its energy and turns in a stalled accretion shock within milliseconds after formation (see the sketch in Figure 1). Despite the decades of effort, the details of how to revive the shock and power CCSN explosion remain unclear (e.g., [3, 4] for recent reviews).

Along with the non-trivial, yet conventional, gas-dynamics effects, CCSNe is a very rich problem that comprises many different phenomena. The newly-born NS cools and contracts, releasing $\sim 10^{53}\,\mathrm{erg}$ potential binding energy as neutrino radiation. The shock-compressed matter is sufficiently opaque to absorb a small fraction of these neutrinos (e.g., [5, 6]). The deposited neutrino energy plays a key role in powering the $\sim 10^{51}\,\mathrm{erg}$ explosion111Core-collapse supernovae with explosion energies as high as $\sim 10^{52}\,\mathrm{erg}$ have also been observed. Also known as hypernovae, such explosions are relatively rare and are believed to be powered by the rotational kinetic energy of rapidly rotating protoneutron stars (e.g., [7, 8]).. Neutrino heating leads to negative gradient of entropy, driving vigorous neutrino-driven turbulent convection in the postshock region ([9, 10] for recent reviews), exerting additional pressure behind the shock [11]. In addition, due to the standing accretion shock instability (SASI) [12, 13, 14], the shock undergoes large-scale non-radial oscillations [15, 16], which expands the size of the region subject to net neutrino heating, resulting in higher heating efficiencies.

Recently, it was shown that the convective instabilities that develop in the innermost nuclear-burning shells of massive stars can have an important impact on the explosion [17, 18, 19, 20, 21, 23]. The convective motion in oxygen and silicon shells may persist even during iron core collapse. As the core collapses, these shells contract and, due to the conservation of angular momentum, the velocities amplify as $\propto r^{-1}$, resulting in an increase by a factor of several (e.g., [25, 28]). Upon reaching the shock, these perturbations interact with the supernova shock and generate additional turbulence in the post-shock region. This augments the turbulent pressure behind the shock (e.g., [20]), resulting in more favorable conditions for producing explosion.

The details of how these perturbations interact with the shock affects the properties of the resulting supernova explosion. Using the linear interaction analysis (LIA), Abdikamalov et al. [25, 26] studied the effect of acoustic, entropy, and vorticity perturbations, which are the three components of a generic weak hydrodynamic turbulent flow [27]. Huete et al. [28] improved their models by taking into account the perturbation of nuclear dissociation energy at the shock. They employed long-time asymptotic expressions to compute the turbulent amplification ratios across the shock front for incoming vorticity waves. The impact of these modes on the explosion condition can be assessed using the notion of the critical (i.e., the minimum) neutrino luminosity necessary for driving the explosion [20]. The effect of the entropic-buoyant turbulent perturbations generated by incident vorticity waves was found to reduce the critical luminosity by $\sim$17–24 per cent, which approximately agrees with the results of three-dimensional simulations of CCSNe [29]. The present study is an extension of [28] and it focuses on the linear interaction of the shock with bidimensional single-mode vorticity perturbations. Employing the mathematical formalism used in describing Richtmyer-Meshkov-type flows [30, 31, 32] and perturbed non-reactive and reactive shocks [33, 34, 35, 36, 37, 38], the exact spatial distribution of the rotational and acoustic perturbations generated in the post-shock flow are provided along with the transient evolution of the shock front towards the permanent oscillatory mode, which, akin to non-ideal gases [39, 40], may change the character of the decay when nuclear dissociation is sufficiently high. The effects of the nuclear dissociation energy, the shock strength, and the characteristic frequency are analyzed in the perturbed flow. Finally, due to accelerated pace of stellar collapse, the inner parts of the star collapses faster than the outer parts (e,g, [41]). As a result, the convective vortices undergo substantial stretching during this phase. Such highly elongated vortices are also considered in this work.

The paper is structured as follows: the problem formulation is shown in Section 2, where the base-flow equations are presented and the linear-perturbation version are provided. The resulting Euler equations are integrated in Section 3, where the transient and the long-time response of the shock front is computed. The complete distribution of the perturbations in the post-shock flow is also shown. The transient evolution of the shock front in the limit of highly elongated vortices is addressed in Section 4. A summary of the results is provided in Section 5.

2 Problem Formulation

2.1 Base-flow equations

Consider a shock wave located at radius $r=R_{\rm shock}(t)$ and assume that the shock thickness $l$ is much smaller than the shock radius ($l\ll R_{\rm shock}$). In this thin-shock limit, one can relate hydrodynamic quantities in the postshock region with those in the preshock region using the conservation equations for the mass, momentum and energy:

[TABLE]

Here, the flow ahead of the shock ($r>R_{\rm shock}$) is denoted with subscript 1, while the flow behind ($r<R_{\rm shock}$) is marked with subscript 2. The variable $u^{\prime}$ is the bulk velocity measured in a reference frame at rest with respect to the center of the star, while variables $\rho$, $p$ and $e$ represent the density, pressure, and internal energy. We model stellar matter as a perfect gas with the polytropic index $\gamma=4/3$ on both sides of the shock. Since we consider scenario of in-falling flow and expanding shock, it is natural to define velocities positive, $u_{1}^{\prime}>0$, $\dot{R}_{{\rm shock}}>0$, for the gas moving inwards and for the shock expanding outwards.

In the presence of nuclear dissociation in a thin layer behind the shock front, the change of the internal energy across the shock is

[TABLE]

where $\Delta e_{\mathrm{dis}}$ is the specific energy employed in the nuclear dissociation process. For stalled shock in CCSNe, $\Delta e_{\mathrm{dis}}$ can be parametrized as $\Delta e_{\mathrm{dis}}=\varepsilon\upsilon_{\mathrm{FF}}^{2}/2$, where $\upsilon_{\mathrm{FF}}$ free-fall speed and $\epsilon$ is a dimensionless parameter [42, 43]. In this scenario, $\epsilon$ scales as $\sim 0.67{\cal M}_{1.3}^{-1}\left(R_{\rm shock}/150\ {\rm km}\right)$, which results in $\varepsilon$ typically ranging between $0.2$ and $0.5$ [28]. For flows with vanishing Bernoulli parameter above the shock, one can express $\Delta e_{\mathrm{dis}}$ in terms of the preshock Mach number $M_{1}=u_{1}/a_{1}$ (see [28] for the details of the derivation),

[TABLE]

where $a_{1}=(\gamma_{1}p_{1}/\rho_{1})^{1/2}$ is the sound speed in the preshock region and $u_{1}=u_{1}^{\prime}+\dot{R}_{{\rm shock}}$ is the preshock speed in the reference frame at rest with respect to the shock.

The fluid properties behind the shock can be conveniently expressed as functions of $\varepsilon$ and preshock Mach number $M_{1}$:

[TABLE]

[TABLE]

[TABLE]

where $M_{2}$ is the mean Mach number in the post-shock region and $u_{2}=u_{2}^{\prime}+\dot{R}_{{\rm shock}}$ is the speed of the postshock flow in the reference frame at rest with respect to the shock. The function

[TABLE]

contains the effect of nuclear dissociation. In the limit of vanishing nuclear dissociation, $1-\kappa\sim M_{1}^{-2}$, expressions (1d)-(1f) reduce to the classical Rankine-Hugoniot relations.

The progress of the nuclear dissociation taking place behind the shock can be quantified in terms of variable $z$, such that $z=0$ refers to values right behind the adiabatic shock and $z=1$ refers to the corresponding properties once nuclear dissociation has been accomplished. Normalized with the flow properties right behind the shock, which are directly obtained from (1d)-(1g) with $\varepsilon=0$, the inner properties are computed as a function of $z$ in Figure 2. The effect of nuclear dissociation is to increase the density ratio with the corresponding velocity decrease, and also to slightly increase the downstream pressure. The rate of the heavy-nuclei breaking $\dot{z}$ would provide the characteristic reaction time, whose combination with the characteristic velocity $a_{1}$ can be used to scale dissociation layer thickness $\ell\sim a_{1}/\dot{z}$.

2.2 Linear perturbation analysis

The upstream flow is seldom uniform. The in-falling matter is immersed in a strong gravity field that triggers convective instabilities. Vortex cells are then formed upstream and they perturb the shock front. The postshock flow is correspondingly perturbed. In order to study how the shock wave reacts to perturbations, a canonical case is selected to be studied: the interaction of harmonic vorticity perturbations with the shock wave in the planar fast-reaction limit $R_{\rm shock}\gg\lambda_{y}\sim\lambda_{x}\gg\ell$. Such an interaction is sketched in Figure 3, where the disturbed shock induces pressure, density and velocity changes in the flow downstream.

The amplitude of the velocity perturbations is assumed to be much smaller than the background flow properties. The small dimensionless amplitude factor $\hat{u}_{1}\sim\left(u_{1}-\langle u_{1}\rangle\right)/\langle u_{1}\rangle\ll 1$, is used to scale the preshock and postshock perturbation variables. The incident shear wave in the frame $(x_{1},y_{1})$ comoving with the in-falling fluid particles is expressed as a divergence-free velocity field

[TABLE]

where $\bar{u}_{1}$ and $\bar{v}_{1}$ are the order-unity stream-wise and transverse velocity perturbation components. The angle brackets represent the time-averaged value of the flow variables. The dimensionless vorticity function, associated to the rotational velocity perturbation (1ha)-(1hb), is

[TABLE]

where $\vec{k}=(k_{x},k_{y})$ is the perturbation wavenumber in the pre-shock region, which can be expressed in terms of the wavelengths $k_{x}=2\pi/\lambda_{x}$ and $k_{y}=2\pi/\lambda_{x}$ sketched in Figure 3, or in terms of the incident shear angle $\theta=\tan^{-1}(k_{y}/k_{x})$. For the analysis of the postshock flow, it is most natural to use a reference frame comoving with the postshock flow. For this reason, hereafter, the dimensionless coordinates $x=k_{y}x_{2}$ and $y=k_{y}y_{2}$ and the dimensionless time $\tau=a_{2}k_{y}t$ are used to describe the solution in the postshock region.

As a result of the interaction, the shock ripples and the fluid downstream is correspondingly altered with acoustic and entropic-vorticity waves. The former travels at the speed of sound downstream $a_{2}$ and the latter moves with the fluid particles. In the postshock region, the dimensionless pressure, density and velocity perturbations are defined as order-unity functions

[TABLE]

The adiabatic Euler equations governing the postshock flow are written in terms of these variables. Since $\bar{p}$ and $\bar{v}$ are proportional to $\cos(y)$ and $\sin(y)$, the conservation equations for mass, $x$-momentum, $y$-momentum, and energy

[TABLE]

can be combined to yield an equation for pressure perturbation $\bar{p}$:

[TABLE]

This is a periodically-symmetric two-dimensional wave equation, the solution of which yields the perturbation field in the postshock region.

The problem thus reduces to solving the wave equation (1hl) as a function of $x$ and $\tau$ for $\tau>0$. The spatial domain is bounded by the leading (first) reflected acoustic wave propagating backwards, $x=-\tau$, and the shock wave front traveling upwards, $x=M_{2}\tau$. One of the boundary conditions stems from the assumption of isolated shock, according to which no reflected sonic waves can reach the shock from behind. This is consistent with the large-radius approximation, $R_{\rm shock}k_{y}\gg 1$, and the linear perturbation analysis that neglects the second-order interaction of these sonic waves with the entropic or rotational perturbations downstream. The upstream flow distorting the shock front is what ultimately determines the amplitude of the shock oscillations and the post-shock inhomogeneities.

The boundary conditions at the shock are readily obtained from the linearized Rankine-Hugoniot relations (1a)-(1c) along with the conservation of momentum transverse to the shock front, namely

[TABLE]

where $\dot{\xi}_{s}$ is the time derivative of the dimensionless shock front deformation. As sketched in Figure 3, the transverse velocity induced behind the shock, included in (1hmd), generates a counter-pressure effect to balance the shock ripple. The characteristic shock reaction time would determine the shock evolution, which in turn is affected by the energy absorbed in the nuclear dissociation process. This effect is quantified in the energy equation (1hmc) via functions

[TABLE]

and

[TABLE]

that allow us to differentiate adiabatic shock waves from reacting shock waves such as shocks with nuclear dissociation or detonations. Functions (1hmn) and (1hmo) are equivalent to those provided in [28], with the latter being here corrected by a factor $1/M_{2}^{2}M_{1}^{2}$.

From (1hma)-(1hmd), one can obtain one of the two relations for the shock boundary condition involving $\bar{\xi}_{s}$ and $\bar{p}_{s}$, while the other can be obtained using the material derivative of the longitudinal velocity perturbation behind the shock, namely

[TABLE]

where

[TABLE]

are the factors accompanying the pressure perturbation.

The initial condition for the shock perturbations can be obtained from the requirement that the initial shock is planar, i.e., $\bar{\xi}_{s}=\bar{v}_{s}=0$. Consequently, the initial pressure and streamwise velocity perturbations obey condition $\bar{u}_{s}+\bar{p}_{s}=0$, thus yielding $\bar{p}_{s0}(\sigma_{b}+1)=1$, for the initial pressure perturbation right behind the shock.

3 Results

3.1 Temporal evolution of the shock front

In order to study the transient response of the shock to upstream perturbation, the transformation

[TABLE]

is employed [35, 38]. The initial condition, $\tau=0$, corresponds to the moment when the initial unperturbed shock first meets the incident vorticity perturbations. The $\chi=$const condition represents a planar surface moving in the postshock gas along the $x$ direction, from the weak discontinuity at $x=0$ $(\chi=0)$ to the reacting shock front at $x=M_{2}\tau$ $(\tanh\chi_{s}=M_{2})$. In terms of these variables, equation (1hl) for sound waves reads

[TABLE]

The boundary conditions at the shock front reduces to

[TABLE]

and

[TABLE]

where

[TABLE]

is the characteristic shock oscillation frequency induced by the incident shear wave.

The Laplace transform is conveniently employed to reduce the above system of partial differential equations to an algebraic system. That is, the integral

[TABLE]

applied to the functions defining the shock boundary conditions yields an algebraic system of function $s$. The Laplace transform of $\left.\frac{1}{r}\frac{\partial\bar{p}_{s}}{\partial\chi}\right|_{s}$ can be computed using the isolated boundary condition, namely $\sqrt{s^{2}+1}\mathscr{P}_{s}-\bar{p}_{s0}$. From this, one can obtain the Laplace transform of the pressure perturbation at the shock:

[TABLE]

where

[TABLE]

is the factor accounting for the periodic excitation amplitude, and

[TABLE]

The solution of equation  (1hmps) for the pressure field can be expressed as a combination of the Bessel functions [44], as shown in [28] in this particular context. It is however illustrative to use the inverse of the Laplace transform to the function (1hmpx), which yields

[TABLE]

as the temporal evolution of the shock pressure perturbations, with

[TABLE]

being the auxiliary function. The corresponding $\bar{p}_{s}(\tau)$ is readily given by the variable change $r=\sqrt{1-M_{2}^{2}}\ \tau$.

Self-induced stable oscillations are found to depend on the slope of the Rankine-Hugoniot curve [45, 46, 47, 39, 40, 48] and the corresponding spontaneous acoustic radiation is determined by the condition $\sigma_{c}>\sigma_{b}$. A parametrical study reveals that it would occur only for $\varepsilon<0$, i.e., when the net positive energy release increases with the shock intensity. In this case, as $\sigma_{c}<\sigma_{b}$ the shock will oscillate only with the excitement frequency coming from upstream perturbations, $\omega_{s}=C_{2}M_{2}k_{x}/k_{y}$, thereby yielding an asymptotic response qualitatively similar to the one found for adiabatic shock waves [35]. Nonetheless, the endothermic contribution may have a qualitative impact on the transient evolution towards the long-time dynamics. Although transient evolution always decays in time with $\tau^{-3/2}$ for $\varepsilon>0$, when the function

[TABLE]

is lower than unity (or $\Lambda<0$ in [39]), the initial degree of damping is significantly modified. That is, for $\eta<1$, corresponding highly endothermic shocks, the oscillations associated to the transient response are effectively shorten, while the contrary occurs for $\eta>1$, the latter case corresponding to the regular shock family of solutions. The delimiting curve $\eta=1$ is computed in Figure 4 for $\gamma=4/3$ as a function of the shock strength $M_{1}$ and the endothermic parameter $\varepsilon$.

Irrespective of the transient behavior, the long-time reaction of the shock pressure to mono-frequency perturbations is

[TABLE]

where

[TABLE]

for $\zeta<1$, and

[TABLE]

for $\zeta>1$. As in previous studies of the interaction of shocks with vorticity perturbations [49, 50, 51, 35, 28], in the long wavelength regime ($\zeta<1$), the sonic disturbances immediately behind the shock consists of two orthogonal contributions $\mathcal{P}_{lr}$, and $\mathcal{P}_{li}$. In this regime, the sonic waves decay exponentially as they move away from the shock. In contrast to this, in the short wavelength regime ($\zeta>1$), the solution is represented by constant-amplitude sonic waves. The critical $\zeta=1$ value corresponds to the case when stable acoustic waves travel parallel to the shock surface in the shock reference frame.

Likewise, the integration of (1hmpb) with respect to $\tau$ yields the temporal evolution of the shock oscillation amplitude:

[TABLE]

with the variable $r=\sqrt{1-M_{2}^{2}}\ \tau$. The associated long-time function of the shock oscillations is

[TABLE]

where the coefficients $\mathcal{J}$ are obtained from the pressure fluctuations at the shock:

[TABLE]

A direct comparison of the long-time response, provided by (1hmpah), and the exact temporal evolution, given by (1hmpag), is computed in Figure 5 as a function of $\tau$ for a shock strength $M_{1}=5$ with two different dissociation sensitivities $\varepsilon=0$ (a) and $\varepsilon=0.4$ (b). The figure is qualitatively similar to figure 5 in [28], where different input parameters had been selected in this occasion. The upstream shear-wave angle is $\theta=60^{\circ}$, which yields dimensionless oscillation frequencies greater than unity in both cases: $\zeta=1.4$ and $\zeta=1.72$ for the adiabatic and endothermic case, respectively. The transient solution is found to achieve the asymptotic regime in a relatively short period of time, with the panel on the left doing it seemingly faster, in agreement with Figure 4. As expected, the amplitude of the long-time oscillations is found to be smaller in the endothermic case.

Figure 6 shows the amplitude of the oscillations as a function of the shock strength $M_{1}$, shear-wave angle $\theta$, and for three different values of the dissociation degree $\varepsilon=0$ (a), $\varepsilon=0.2$ (b), and $\varepsilon=0.4$ (c). In agreement with Figure 5, the amplitude is found to be generally smaller for endothermic shocks. The limit $\zeta=1$ is also computed in Figure 6 as a function of $M_{1}$ and $\theta$ for different endothermic intensities. The zones on the left of these dashed curves correspond to pressure radiating conditions (high-frequency regime), and the zones on the right refer to non-radiating conditions (low-frequency regime), with the area of the latter being reduced with the increase of nuclear dissociation.

3.2 Post-shock flow variables

The downstream flow perturbations are conveniently split decomposed those of acoustic type, which travel at the speed of sound relative to the fluid particles, and those of entropic-rotational nature, which move with the fluid particles [27].

The exact temporal evolution of the pressure field downstream is readily obtained through Bessel functions, as derived in [28], or by direct integration of the wave equation provided that $\bar{p}(x=M_{2}\tau)=\bar{p}_{s}(\tau)$. The former strategy has been employed in Figure 7 to compute the pressure field in the domain $0\leq x\leq M_{2}\tau$. Computations reveal two well-distinguished regimes: acoustically radiating and non-radiating conditions. The acoustic radiation condition is then determined by $\omega_{s}>(1-M_{2}^{2})^{1/2}$, a condition that depends on the upstream shear wave, since $\zeta=\left[0,\infty\right)$ depends on the relative properties of the perturbation field ahead of the shock. Small values of $\zeta$ represent the interaction with upstream vortices highly stretched in the streamwise direction $\lambda_{x}\gg\lambda_{y}$, while the opposite is true for $\zeta\gg 1$. In the latter low mode-number scenario ($\lambda_{x}\ll\lambda_{y}$), the problem reduces to the one-dimensional interaction of the shock with radial perturbation waves.

The sonic waves traveling in the postshock region are functions of $(\omega_{a}\tau-k_{a}x)$, where the frequency $\omega_{a}$ and the wavenumber $k_{a}$ are obtained from the shock oscillation frequency $\omega_{s}=\omega_{a}-M_{2}k_{a}$ and the adiabatic dispersion relation $\omega_{a}^{2}=k_{a}^{2}+1$:

[TABLE]

It is readily seen that $k_{a}$ can be either positive or negative. The $k_{a}<0$ case represents acoustic waves traveling in the direction of the postshock flow, while the waves moving in the opposite direction towards the shock have $k_{a}>0$. The solution corresponding to shock oscillation frequency $\omega_{s}=1$ represents the standing acoustic waves that separate the solution traveling to the left $\omega_{s}>1$ from the one traveling to the right $(1-M_{2}^{2})^{1/2}<\omega_{s}<1$ in the reference frame comoving with the postshock fluid. At large distances from the shock in the downstream region (far larger than the wavelength of the perturbations), the asymptotic pressure and the isentropic density perturbations are given by

[TABLE]

where $\mathcal{P}_{s}$ is the amplitude of the pressure perturbations. The isentropic temperature variations induced by the acoustic radiation are just $\bar{T}_{a}(x,\tau)=\left(\gamma-1\right)\bar{p}(x,\tau)$.

The acoustic contribution of the velocity perturbations are readily obtained through the Euler momentum equations

[TABLE]

which can be used to write the long-time response

[TABLE]

where

[TABLE]

are the associated amplitudes.

In absence of dissipative effects, vorticity disturbances downstream remain frozen to the fluid particles, with the amplitude being determined by the vorticity produced at the shock, namely

[TABLE]

where

[TABLE]

accounts for the shock-compression of the vortices, a one-dimensional effect, and

[TABLE]

indicates the contribution of shock deformations, which is a two-dimensional effect.

The linear perturbations in the velocity field satisfy the equation

[TABLE]

with the breakdown of irrotational-acoustic and steady-rotational perturbations obeying, separately,

[TABLE]

With the acoustic field being given by (1hmpam) and (1hmpan), the spatial distribution of the rotational-velocity perturbations is now derived to provide the complete velocity field. The solenoidal part is calculated by tracking the vorticity left behind by the oscillating shock front from $\tau=0$, yielding

[TABLE]

with the particular solution $\bar{u}_{p}$ being provided by

[TABLE]

and with the auxiliary integration function $f(z)$ being defined in (1hmpab).

The function (1hmpauav) is qualitatively similar to that shown in [35] and [38] for adiabatic and reacting shocks waves, respectively. Details of its derivation, omitted here for the sake of conciseness, involve the use of the inverse Laplace transform technique. The asymptotic rotational contribution of the velocity field is also written as a piecewise function of the shock oscillation frequency, with the longitudinal

[TABLE]

and transverse

[TABLE]

contributions being characterized by their amplitudes

[TABLE]

and

[TABLE]

for the streamwise and crosswise components in the long-wavelength regime, respectively, and

[TABLE]

for the corresponding short-wavelength regime amplitudes.

The density variations $\bar{\rho}_{e}$ due to the entropy waves is obtained from Rankine-Hugoniot relations (1hma)-(1hmc) by subtracting the contribution of sonic waves:

[TABLE]

where $\mathcal{D}=\left(2M_{2}\sigma_{b}-1\right)/M_{2}^{2}$ is the density perturbation amplitude relative to the pressure at the shock. The temperature variation corresponding to the entropy waves is given by $\bar{T}_{e}(x)=-\bar{\rho}_{e}(x)=-(\mathcal{D}-1)\bar{p}_{s}(\tau=x/M_{2})$, where the temperature is normalized with the base flow temperature.

Figure 8 shows the spatial distribution of the rotational velocity field and entropic density perturbations from the origin $x=0$ for $M_{1}=2$, $\varepsilon=0.2$ and $\zeta=1.5$. The upper panels (a)-(c) display the functions $\bar{u}_{r}$, $\bar{v}_{r}$, and $\bar{\rho}_{e}$ as a function of $x$, respectively, and the lower panel (d) shows the vector velocity field superposed to the density iso-contours in the plane $(x,y)$. The exact solution for the rotational-velocity contribution (solid line) is found to approach the asymptotic solution (red-dashed line) in a fairly short distance, while the entropic-density function exhibits a longer transient period towards the long-time solution.

The amplitude of the asymptotic rotational and acoustic velocity perturbations, as well as the corresponding entropic and acoustic density perturbations, have been conveniently expressed in terms of the shock pressure perturbation amplitude. Likewise, the asymptotic shock oscillation amplitude has been written as a function of $\mathcal{P}$ in equation (1hmpai). It is then immediate to obtain the value of any perturbation variable with the aid of Figure 6. Further computations, as the long-time amplitudes for the velocity and density perturbations, have been computed in Figures A1-A3 of [28]. As found in Figure 6, there exists a peak in the perturbation amplitude near the critical frequency $\zeta=1$ and the effect of endothermicity, along with the associated amplitude change, is to stretch the peak in the frequency domain. This effect occurs for any perturbation variable: pressure, velocity or density. The acoustic contribution is found to provide a negligible contribution in comparison to the rotational or entropic counterpart.

Density perturbations generated by the interaction of upstream asphericities with the shock are found to play a pivotal role in driving post-shock turbulence by buoyancy effects, which translates into a reduction of the critical neutrino luminosity necessary for producing explosion [9, 52]. In [28], the effect of vorticity waves with isotropic distribution of orientations was considered. This is perhaps not an accurate representation. The convective motion is characterized by a dominant eddy with a specific size that undergo stretching due to the accelerated collapse, which increases the value of the shear angle $\theta$.

Anticipating that the entropic contribution dominates the post-shock density perturbations, the correlation between the upstream turbulent Mach number and the entropic-density perturbations is used to evaluate the contribution of the shock-generated density fluctuations through the factor

[TABLE]

which is computed in Figure 9 as a function of the upstream Mach number and the shear-wave angle $\tan\theta=k_{y}/k_{x}$. The upper region, corresponding to the highly-elongated vorticity limit $\zeta\ll 1$ ($\theta\sim 90$), is found to yield stronger density perturbations than those produced for $\zeta\gg 1$ ($\theta\sim 0$), the lower limit. Modification of the endothermic contribution $\varepsilon$ does not change the qualitative picture significantly, but it increases the figures. That is, despite pressure perturbations and shock-ripple deviations decrease with the endothermicity associated to the nuclear dissociation mechanism, entropic-density perturbations grow with the factor $\varepsilon$.

4 Interaction with highly-stretched vortices

Natural convective cells are typically formed when temperature gradients are counter-aligned with the gravity field, conditions that are met in massive stars. When, in addition, they are advected by a highly-accelerating radially-converging flow, they elongate in the radial direction. In this scenario, the upstream perturbation field is strongly anisotropic, so that the canonical shock-turbulence interaction approach is no longer applicable. It is, however, possible to construct a representative problem setup with the formulation previously presented, as depicted in Figure 10.

In the planar-shock limit, i.e., when the characteristic length of the upstream perturbation field is much smaller than the shock radius, the formulation of the problem shown in previous section can be particularized to the case $k_{x}\ll k_{y}$ or, equivalently, $\zeta\ll 1$. In this slender limit, the acoustically-induced shock oscillation period is much shorter than the characteristic residence time of the convective cells crossing out the shock, thereby providing the following upstream modulation

[TABLE]

as the transverse contribution becomes negligible.

The corresponding boundary conditions at the shock front are appropriately adapted to yield

[TABLE]

and

[TABLE]

which have been conveniently written as a function of the variables $r$ and $\chi$ [34]. After some straightforward manipulation, it is obtained the Laplace Transform of the pressure perturbations, namely

[TABLE]

provided that $\left(1+\sigma_{b}\right)\bar{p}_{s0}=1$, and with the factor in (1hmpy) being now

[TABLE]

The initially planar shape of the shock front is distorted as a result of the interaction with the upstream modulated velocity field. For $y=0$, the shock encounters an always positive velocity perturbation aligned with the shock propagation, which pulls the shock front upwards. The opposite would apply for $k^{\prime}y_{1}=y=\pi$, where an always negative velocity input pushes the front backwards. Then, the tangential velocity perturbation induced by the shock corrugation, with the associated mass flux, tends to restore the shock shape. This two counter-effects are responsible of the shock oscillations in the initial stage. Eventually, the shock approaches a steady-state regime when the two opposed effects balance, something that cannot occur when the upstream non-uniform flow is made of vortices with $k_{x}\sim k_{y}$. The corresponding long-time functions take the forms $\bar{p}_{s}\left(\tau\gg 1,y\right)=\bar{p}_{s}^{\infty}\cos(y)$ and $\xi_{s}\left(\tau\gg 1,y\right)=\xi_{s}^{\infty}\cos(y)$ for the shock pressure and ripple variables, respectively, with the associated amplitudes being

[TABLE]

The evolution of the shock pressure and shock ripple perturbations, when the front travels through a stripe-like pattern of longitudinal perturbations, is computed in Figure 11 for $M_{1}=5$, and for $\varepsilon=0.0$ and $\varepsilon=0.4$. The transient behaviour is readily obtained by just taking the limit $\zeta-1\ll 1$ in (1hmpaa) and (1hmpag), respectively. It is observed that endothermic nuclear dissociation exhibits a stabilizing mechanism, since the long-time amplitudes are smaller and the oscillations decay faster for $\varepsilon=0.4$ than for $\varepsilon=0$. This is in consonance with Figure 4, where the case $M_{1}=5$ and $\varepsilon=0.4$ lies on the region $\eta<1$, while regular adiabatic shocks with $\varepsilon=0$ lie on the zone $\eta>1$.

In regard to the far-field perturbations, it is immediate to see that the acoustic contribution is negligible, as pressure perturbations decay exponentially with the distance from the shock for $\zeta<1$. What remains constant, in absence of diffusive effects, is the entropy-vorticity perturbation generated across the shock, namely

[TABLE]

Likewise, by direct inspection of (1hmpaub) and anticipating that tangential velocity generated behind the shock when $\bar{u}_{r}=$constant is of acoustic type, and then evanescent, it is found that

[TABLE]

for the longitudinal and transverse components of the rotational velocity field.

5 Conclusions

When a shock wave encounters a vorticity wave on its way, their interaction results in a deformation of the shock. The latter induces pressure changes that results in the radiation of acoustic waves downstream. Along with the sonic mode, post-shock perturbations include vortical and entropic disturbances that are convected by the fluid particles downstream. The distinguished feature of nuclear-dissociating shocks is the endothermic contribution, which depends on the shock intensity relative to the flow stream. Then, in a likely non-uniform context, perturbations ahead of the shock as the shear-pattern considered in this work, the amount of nuclei that are dissociated is affected by the shock perturbation, which in turn affects the energy balance across the shock, and ultimately the amplitude of the perturbations downstream. This effect can be easily studied when considering the interaction of the shock with intermediate vortical scales, i.e., those whose characteristic length is sufficiently small for the shock to be considered a planar front, yet sufficiently large for the shock to be a seen as a discontinuity front.

This problem is particularly relevant in the context of core-collapse supernova explosions, where the shock after bounce swallows the convective structures generated upstream within the supersonic inwards-traveling mass. The shock-perturbation interaction modifies and creates additional perturbations downstream, which affects the critical conditions for the supernova explosion. When the upstream flow is assumed to be dominantly isotropic, theory predicts that the injection of non-radial motion and the buoyancy-driven convection triggered by entropy waves reduce the critical neutrino luminosity by$\sim$ 12-24 per cent, for typical problem parameters [9, 25, 28]. When convective cells are not sufficiently turbulized by the inwards acceleration of the fluid particles, the isotropic assumption may be inaccurate. Motivated by this fact, in this work the emphasis is placed in bidimensional structures made of iso-density vortices. The temporal evolution of the oscillating shock and the long-time asymptotic expressions have been derived analytically as a function of the dominant governing parameters: the shock strength $M_{1}$, the nuclear dissociation degree $\varepsilon$, and the incident shear-wave angle $\sim k_{x}/k_{y}$. Likewise, the exact and asymptotic spatial distribution of the perturbations in the shocked gas have been derived explicitly with use made of the Laplace Transform technique. The limit of high elongated vortices, corresponding to $k_{x}\ll k_{y}$, is also evaluated in terms of closed-form expressions.

The effect of the endothermic nuclear dissociation is found to diminish the amplitude of the shock oscillations and to reduce the acoustic radiation. Contrarily, the entropic density disturbances grow with the factor $\varepsilon$. In some conditions, for sufficiently endothermic shock, the transient evolution is significantly shorten. The distribution of the downstream entropic and vortical perturbations have been provided analytically. These perturbations serve as initial conditions of the post-shock regime, as density variances become buoyant in a temporal scale that exceeds the shock influence. This phenomenon will be studied in more depth in a future work.

Acknowledgements

This work is supported by the Ministry of Science, MEC (ENE2015-65852-C2-1-R) and Fundación Iberdrola España (BINV-hBbhOeJQ), Spain (for C. Huete), by MES RK state-targeted program BR05236454, MES RK grant No. 346, NU ORAU grant SST 2015021 and NU grant No. 090118FD5348 (for E. Abdikamalov).