LINKS BETWEEN THE SHOCK INSTABILITY IN CORE-COLLAPSE SUPERNOVAE AND ASYMMETRIC ACCRETIONS OF ENVELOPES

The explosion mechanism of core-collapse supernovae (CCSNe) is one of the biggest interests in high-energy astrophysics. The challenge is how the stagnant shock wave revives and goes through an iron core. The most anticipated scenario is based on the so-called neutrino-heating mechanism, in which fluid behind the shock gains energy from neutrinos that diffuse out of the proto-neutron star (PNS), and as a result the shock is pushed outward. In this scenario, multidimensional fluid instabilities such as the standing accretion shock instability (SASI) and/or neutrino-driven convections are expected to play an important role in enhancing the efficiency of the neutrino heating and producing shock revival. In fact, it is almost a consensus that one-dimensional, spherically symmetric simulations do not lead to shock revival even with the most detailed physics implemented.

The multidimensional fluid instabilities in CCSNe have been systematically investigated by numerical simulations (Burrows et al. 2012; Bruenn et al. 2013; Hanke et al. 2013; Iwakami et al. 2014a, 2014b; Takiwaki et al. 2014; Abdikamalov et al. 2015; Fernández 2015) and also by linear analysis (Foglizzo et al. 2007; Yamasaki & Yamada 2007; Yamasaki & Foglizzo 2008; Foglizzo 2009; Sato et al. 2009; Guilet & Foglizzo 2012; Foglizzo et al. 2015). The neutrino-driven convection is characterized by small-scale multiple buoyant bubbles and is driven by negative entropy gradients caused by neutrino heating. SASI, on the other hand, is characterized by large-scale deformations of the shock surface induced by sloshing and spiral motions. It is a global instability caused by the so-called advective-acoustic cycle, in which the spherical shock surface is deformed by repetitious round trips of upcoming sonic waves and down-going vorticity waves between the shock and PNS surface. Which of these instabilities, convection or SASI, dominates the shock dynamics is dependent on the mass accretion rate, neutrino luminosity, and hence the progenitors.

Multidimensionality is recently discussed not only in fluid instabilities but also in structures of progenitors, which was revealed by a series of papers by the groups of Arnett (Bazan & Arnett 1998; Asida & Arnett 2000; Meakin & Arnett 2006, 2007; Arnett & Meakin 2011) and Chatzopoulos et al. (2014). They showed numerically that the silicon and oxygen shells outside the iron core in progenitors substantially deviate from spherically symmetric configurations at the onset of collapse due to nuclear burnings and associated violent convections. They have found that the power spectrum of the turbulence peaks at l = 4 or 5. Another group reported, however, the development of a larger-scale convection with l = 2 in a different progenitor just before the collapse (Müller et al. 2016). Their results imply that the patterns of the fluctuation may be dependent on the progenitor.

Such fluctuations in the outer shells will accrete onto the stagnant shock and may affect the subsequent shock dynamics. In fact, it is demonstrated by linear analysis (Lai & Goldreich 2000; Takahashi & Yamada 2014) that the fluctuations can grow rather than damp during the supersonic infall. Furthermore, Couch & Ott (2013, 2015), Müller & Janka (2014), and Couch et al. (2015) showed that the interaction of such nonspherical outer layers with the stagnant shock can produce shock revival for the progenitors that failed to explode without the fluctuations. The systematic study by Müller & Janka (2014) reported that large-scale velocity perturbations efficiently contribute to the shock revival. Fernández et al. (2014) also found that upstream perturbations with l = 1, 2 do not turn the SASI-dominant flows into convectively dominated flows.

Thanks to the above studies, the multidimensional structures of progenitors have been focused as a new key ingredient to successful explosions of CCSNe in addition to the fluid instabilities and other multidimensional players such as stellar rotation and magnetic fields (Guilet & Foglizzo 2010; Iwakami et al. 2014b; Sawai & Yamada 2014; Mösta et al. 2015; Takiwaki et al. 2016). However, it has not been fully understood yet how the asymmetry in the envelopes interacts with the shock instabilities and eventually helps the stalled shock revive.

In this paper, we focus on the links between the shock instabilities and the asymmetric accretion flows that occur as consequences of the multidimensional structures of progenitors. We perform a linear stability analysis of the standing accretion shock by taking into account such asymmetric, unsteady accretion flows in front of the shock as the outer boundary condition, which has never been considered in previous liner analyses. We systematically study the dependence of the results on the typical frequency of the accreting fluctuations for several types of background flows, which are designed to crudely mimic the collapse of fluctuating progenitors. We also obtain some general relations analytically.

This paper is organized as follows. We explain the method and setups in the next section, where basic equations and Laplace transform are introduced together with our models. In Section 3, we present the results of the global linear stability analysis of the standing accretion shock. Discussions are given in Section 4. We summarize our investigation in Section 5. Following the main body, some appendices are added, which include analytic treatments that give the general relations between shock dynamics and perturbations.

As stated in the introduction, we investigate the linear stability of the standing accretion shock against the multidimensional upstream perturbations. We ignore the time variation of the (spherically symmetric, unperturbed) flows, since the typical timescales of the changes in the neutrino luminosity and/or mass accretion rate are much longer than the timescale of the instability.

In this section, we introduce the system of basic equations and the method of Laplace transform, which is the principal tool for the investigation of the linear analysis in this paper.

2.1. Basic Equations

The basic equations that govern accretion flows in the SN core are given as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(\rho {\boldsymbol{v}})+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{vv}}+p{\boldsymbol{I}})=-\rho \displaystyle \frac{{GM}}{{r}^{2}}\displaystyle \frac{{\boldsymbol{r}}}{r},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\varepsilon }{{dt}}+p\displaystyle \frac{d}{{dt}}\left(\displaystyle \frac{1}{\rho }\right)=q,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}({{nY}}_{e})+{\boldsymbol{\nabla }}\cdot ({{nY}}_{e}{\boldsymbol{v}})=\lambda ,\end{eqnarray} \tag{ 4 }$$

in addition to an equation of state (EOS). In the above expression, ρ, p, n, Y_e, ε, and ${\boldsymbol{v}}$ are density, pressure, number density, electron fraction, specific internal energy, and velocity, respectively. The mass of the PNS is denoted by M, which is assumed to be constant, and G is the gravitational constant. The self-gravity is neglected. We incorporate the reactions between the electron-type neutrinos and matter, which are symbolically denoted by q and λ: the former is the net heating rate, and the latter is the net reaction rate of electrons and positrons.

We use the lightbulb approximation instead of solving the neutrino transport (Ohnishi et al. 2006; Scheck et al. 2006). In this prescription, neutrino luminosities (${L}_{{\nu }_{e}}$ and ${L}_{{\bar{\nu }}_{e}}$) and temperatures (${T}_{{\nu }_{e}}$ and ${T}_{{\bar{\nu }}_{e}}$) are arbitrary parameters. Then the radius of the neutrino sphere, ${r}_{{\nu }_{e}}$, is given by the following relation:

$$\begin{eqnarray}&&{L}_{{\nu }_{e}}=\displaystyle \frac{7}{16}4\pi {r}_{{\nu }_{e}}^{2}\sigma {T}_{{\nu }_{e}}^{4},\end{eqnarray} \tag{ 5 }$$

where σ is the Stefan–Boltzmann constant. The radius of the antielectron neutrino sphere is also given in the same manner. The unperturbed flows are given as spherically symmetric, time-independent solutions to these equations with appropriate boundary conditions.

Following Lai & Goldreich (2000) and our previous paper, we linearize these equations for perturbed quantities written as

$$\begin{eqnarray}&&\delta X({\boldsymbol{r}},t)=\sum _{l,m}\delta {X}^{(l,m)}(r,t){Y}_{{lm}}(\theta ,\phi ),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\delta {\boldsymbol{v}}({\boldsymbol{r}},t) & = & \displaystyle \sum _{l,m}\delta {v}_{r}^{(l,m)}(r,t){Y}_{{lm}}(\theta ,\phi )\hat{{\boldsymbol{r}}}\\ & & +\delta {v}_{\perp }^{(l,m)}(r,t)\left[\hat{{\boldsymbol{\theta }}}\displaystyle \frac{\partial {Y}_{{lm}}}{\partial \theta }+\displaystyle \frac{\hat{{\boldsymbol{\phi }}}}{\sin \theta }\displaystyle \frac{\partial {Y}_{{lm}}}{\partial \phi }\right]\\ & & +\delta {v}_{\mathrm{rot}}^{(l,m)}(r,t)\left[-\hat{{\boldsymbol{\phi }}}\displaystyle \frac{\partial {Y}_{{lm}}}{\partial \theta }+\displaystyle \frac{\hat{{\boldsymbol{\theta }}}}{\sin \theta }\displaystyle \frac{\partial {Y}_{{lm}}}{\partial \phi }\right],\end{eqnarray} \tag{ 7 }$$

where X denotes scalar variables. ${Y}_{{lm}}(\theta ,\phi )$ is the spherical harmonics with the indices $l,m$. $\hat{{\boldsymbol{r}}},\hat{{\boldsymbol{\theta }}}$, and $\hat{{\boldsymbol{\phi }}}$ are the unit vectors in spherical coordinates. For spherically symmetric background flows, the system of linearized equations for each combination of indices, (l, m), is decoupled from each other and is given symbolically as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{y}}}{\partial r}=A\displaystyle \frac{\partial {\boldsymbol{y}}}{\partial t}+B{\boldsymbol{y}},\end{eqnarray} \tag{ 8 }$$

where ${\boldsymbol{y}}={\boldsymbol{y}}(r,t)$ denotes the vector of perturbed quantities given as

$$\begin{eqnarray}&&{\boldsymbol{y}}(r,t)={\left(\displaystyle \frac{\delta \rho }{{\rho }_{0}},\displaystyle \frac{\delta {v}_{r}}{{v}_{r0}},\displaystyle \frac{\delta {v}_{\perp }}{{v}_{r0}},\displaystyle \frac{\delta \varepsilon }{{\varepsilon }_{0}},\displaystyle \frac{\delta {Y}_{e}}{{Y}_{e0}},\displaystyle \frac{\delta {v}_{\mathrm{rot}}}{{v}_{r0}}\right)}^{T},\end{eqnarray} \tag{ 9 }$$

where ${(\ldots )}^{T}$ means a transposed vector and the subscript 0 denotes unperturbed states. Here and hereafter the superscript of indices is omitted for notational simplicity. The coefficient matrices, $A=A(r)$ and $B=B(r)$, whose components are made of the background quantities, are given in Appendix A. While the matrices are dependent on l, they are independent of m thanks to the assumption of spherically symmetric background.

We solve the linearized Equation (8) in the region between the standing shock and the PNS surface, which is assumed to coincide with the electron-type neutrino sphere. The initial perturbed state at t = 0 in the shocked region is given as

$$\begin{eqnarray}&&{\boldsymbol{y}}(r,t=0)={{\boldsymbol{y}}}_{0}(r)\ ({r}_{{\nu }_{e}}\lt r\lt {r}_{\mathrm{sh}}),\end{eqnarray} \tag{ 10 }$$

where r_sh is the unperturbed shock radius.

The outer boundary condition imposed at the shock radius is given by the perturbed quantities in front of the shock through the linearized Rankine–Hugoniot relations, which are schematically given by

$$\begin{eqnarray}&&{\boldsymbol{y}}({r}_{\mathrm{sh}},t)=R{\boldsymbol{z}}(t)+\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}{\boldsymbol{c}}+\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}{\boldsymbol{d}},\end{eqnarray} \tag{ 11 }$$

where $\delta {r}_{\mathrm{sh}}$ is the time-dependent perturbed shock radius. R is another matrix, whose components are described only by background quantities. The components of the vectors, ${\boldsymbol{c}}$ and ${\boldsymbol{d}}$, are also given by the unperturbed states. Their explicit forms are found in Appendix A. A given fluctuation in front of the shock is denoted by ${\boldsymbol{z}}(t)$.

The inner boundary is set at the PNS surface. We note here that we can impose only one condition, which determines a remaining degree of freedom at the outer boundary, i.e., the perturbed shock radius, once the upstream flow is given. Otherwise, there is generally no solution that satisfies the inner and outer boundary conditions at the same time. Hence, the inner boundary condition is symbolically represented as

$$\begin{eqnarray}&&f({\boldsymbol{y}}({r}_{{\nu }_{e}},t),t)=0.\end{eqnarray} \tag{ 12 }$$

To summarize, we solve the initial boundary value problem that is described by Equations (8) and (10)–(12) to find the time-dependent shock radius, $\delta {r}_{\mathrm{sh}}/{r}_{\mathrm{sh}}$.

2.2. Mode Analysis in Laplace Transform

To solve this initial boundary value problem, we use the Laplace transform with respect to time (for properties of Laplace transform see Schiff 1999; Takahashi & Yamada 2014). Laplace-transformed equations are given as

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{y}}}^{* }}{{dr}}(r,s)=({sA}+B){{\boldsymbol{y}}}^{* }-A{{\boldsymbol{y}}}_{0}(r),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\boldsymbol{y}}}^{* }({r}_{\mathrm{sh}},s)=(s{\boldsymbol{c}}+{\boldsymbol{d}})\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }(s)}{{r}_{\mathrm{sh}}}+R{{\boldsymbol{z}}}^{* }(s),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{f}^{* }\,({{\boldsymbol{y}}}^{* }({r}_{{\nu }_{e}},s),s)=0,\end{eqnarray} \tag{ 15 }$$

where s is the conjugate variable of t and is a complex number. The superscript ∗ means Laplace-transformed variables, which are complex functions of s in general.

We here emphasize that this is a system of ordinary differential equations with s being a parameter. We find $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$ for each s so that both the inner and outer boundary conditions would be satisfied. We note that such a value of $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$ is easily found for a given s by integrating Equation (13) twice as explained in Appendix B.

The eigenmodes, both stable or unstable, are found in the following way. We first assume the functional form of the perturbed shock radius as

$$\begin{eqnarray}&&\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}=\sum _{j}{a}_{j}{e}^{{{\rm{\Omega }}}_{j}t}\sin ({\omega }_{j}t+{\phi }_{j}),\end{eqnarray} \tag{ 16 }$$

where ${{\rm{\Omega }}}_{j}$ and ${\omega }_{j}$ are the growth/damping and oscillation frequency of the jth mode ($j=1,2,3,\ldots$), respectively, and a_j is the amplitude. To ensure the uniqueness of this expansion, we assume that ${\omega }_{j}\geqslant 0$ and $-\pi /2\leqslant {\phi }_{j}\lt \pi /2$ for any j. Performing Laplace transform, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}=\sum _{j}{a}_{j}\displaystyle \frac{(s-{{\rm{\Omega }}}_{j})\sin {\phi }_{j}+{\omega }_{j}\cos {\phi }_{j}}{{(s-{{\rm{\Omega }}}_{j})}^{2}+{\omega }_{j}^{2}},\end{eqnarray} \tag{ 17 }$$

which has single poles at $s={{\rm{\Omega }}}_{j}\pm i{\omega }_{j}$. Here we used the relation

$$\begin{eqnarray}&&{ \mathcal L }[{e}^{{\rm{\Omega }}t}\sin (\omega t+\phi )]=\displaystyle \frac{(s-{\rm{\Omega }})\sin \phi +\omega \cos \phi }{{(s-{\rm{\Omega }})}^{2}+{\omega }^{2}}.\end{eqnarray} \tag{ 18 }$$

Hence, we know the growth/damping rate and eigenfrequency of a mode from the position of the corresponding pole. The poles of $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$ are found, on the other hand, by solving the initial boundary value problem (13)–(15) for a series of s. Since the functional value changes near a pole, as shown in Figure 1, for example, we can easily find poles by observing the behavior of the function as long as the grid on the complex plane is sufficiently fine.⁶

Zoom In Zoom Out Reset image size

Incidentally, we obtain the time evolution of the shock radius by performing inverse Laplace transform for $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$ numerically. Note that the inverse transform should be done by an integral along a path that is parallel to the imaginary axis and that all singularities lie to the left of the line (Schiff 1999; see also Appendix C). Furthermore, the obtained function is necessarily 0 for $t\lt 0$. If an inappropriate path is chosen for the integral, i.e., if some poles are missed, the time evolution will not remain stationary at $t\lt 0$. We would hence never miss the pole with the maximal growth rate if we look into the time evolution at $t\lt 0$, which is one of the advantages of the Laplace-transform method.

2.3. Models and Parameters

The initial boundary value problem (13)–(15) is solved for various background flows and perturbations that accrete onto the shock surface. A background state is characterized by the following parameters: Throughout the paper, a constant mass of the PNS, M = 1.4 M_⊙, is used. The mass accretion rate and neutrino luminosities are fixed as $\dot{M}=0.6$M_⊙ s⁻¹ and ${T}_{{\nu }_{e}}={T}_{{\bar{\nu }}_{e}}=4.5\,{\rm{MeV}}$. The accretion flow ahead of the shock is assumed to be a free fall of irons whose entropy and electron fraction are given as $S=3{k}_{B}$ and Y_e = 0.5, where k_B is the Boltzmann constant. We employ Shen's EOS (Shen et al. 1998), which takes into account the contributions from nucleons, nuclei, photons, electron, and positrons. For neutrino heating and cooling, the reaction rates of Bruenn (1985) are adopted. The inner boundary is assumed to coincide with the neutrino sphere, where the density is fixed to 10¹¹ g cm⁻³. Neutrino luminosities are also model parameters, and some value is given to ${L}_{\nu }:= {L}_{{\nu }_{e}}={L}_{{\bar{\nu }}_{e}}$ for each.

We list in Table 1 the background-flow models that are representatively used in the study. In the table, ${\omega }_{\mathrm{aac}}$ and ${\omega }_{\mathrm{pac}}$ are characteristic frequencies of the advective-acoustic cycle and purely acoustic cycle (Blondin et al. 2003), respectively, which are given as

$$\begin{eqnarray}&&{\omega }_{\mathrm{aac}}=2\pi {\left[{\int }_{{r}_{\nu }}^{{r}_{\mathrm{sh}}}\left(\displaystyle \frac{1}{| {v}_{r}| }+\displaystyle \frac{1}{{c}_{s}-| {v}_{r}| }\right){dr}\right]}^{-1},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\omega }_{\mathrm{pac}}=2\pi {\left[{\int }_{{r}_{\nu }}^{{r}_{\mathrm{sh}}}\left(\displaystyle \frac{1}{{c}_{s}+| {v}_{r}| }+\displaystyle \frac{1}{{c}_{s}-| {v}_{r}| }\right){dr}\right]}^{-1},\end{eqnarray} \tag{ 20 }$$

where c_s is the sound speed and ${r}_{\nu }={r}_{{\nu }_{e}}={r}_{{\bar{\nu }}_{e}}$ is the neutrino sphere. We also gave the χ-parameter (Foglizzo et al. 2006) in the list, which is defined as

$$\begin{eqnarray}&&\chi ={\int }_{{r}_{\mathrm{gain}}}^{{r}_{\mathrm{sh}}}\left|\displaystyle \frac{{\omega }_{\mathrm{BW}}}{{v}_{r}}\right|{dr},\end{eqnarray} \tag{ 21 }$$

where r_gain is the gain radius, i.e., the bottom boundary of the region with negative entropy gradients (${r}_{\mathrm{gain}}\leqslant r\leqslant {r}_{\mathrm{sh}}$), and ${\omega }_{\mathrm{BW}}$ is the Brunt–Väisälä frequency:

$$\begin{eqnarray}&&{\omega }_{\mathrm{BV}}=\sqrt{\displaystyle \frac{{GM}}{{r}^{2}}\left|\displaystyle \frac{1}{{{\rm{\Gamma }}}_{1}p}\displaystyle \frac{{dp}}{{dr}}-\displaystyle \frac{1}{\rho }\displaystyle \frac{d\rho }{{dr}}\right|},\end{eqnarray} \tag{ 22 }$$

with

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1}:= {\left(\displaystyle \frac{\partial \mathrm{ln}p}{\partial \mathrm{ln}\rho }\right)}_{S,{Y}_{e}}.\end{eqnarray} \tag{ 23 }$$

According to Foglizzo et al. (2006), $\chi \gtrsim 3$ is the condition for the flow being convectively unstable. We also list the mean Brunt–Väisälä frequency for models with a nonvanishing gain region, which is defined by

$$\begin{eqnarray}&&{\bar{\omega }}_{\mathrm{BV}}=\displaystyle \frac{1}{{r}_{\mathrm{sh}}-{r}_{\mathrm{gain}}}{\int }_{{r}_{\mathrm{gain}}}^{{r}_{\mathrm{sh}}}{\omega }_{\mathrm{BV}}{dr}.\end{eqnarray} \tag{ 24 }$$

Table 1.  Background-flow Models

Model	Lν	rsh	rν	${\omega }_{\mathrm{aac}}$	${\omega }_{\mathrm{pac}}$	rgain	${\bar{\omega }}_{\mathrm{BV}}$	χ
	(1052 erg s−1)	(106 cm)	(106 cm)	(ms−1)	(ms−1)	(106 cm)	(ms−1)
M06L20	2.0	4.75	2.94	1.56	5.90	4.75	⋯	0
M06L25	2.5	5.79	3.28	0.973	3.99	5.79	⋯	0
M06L30	3.0	6.93	3.60	0.644	2.84	6.45	0.270	0.143
M06L35	3.5	8.22	3.89	0.444	2.07	6.80	0.320	0.603
M06L40	4.0	9.78	4.15	0.313	1.52	7.18	0.340	1.41
M06L45	4.5	11.8	4.41	0.221	1.10	7.73	0.332	2.56
M06L50	5.0	14.6	4.64	0.152	0.754	8.44	0.302	4.31
M06L55	5.5	19.3	4.87	0.0974	0.470	9.42	0.238	6.77
M06L60	6.0	29.5	5.09	0.0510	0.232	10.8	0.153	11.1

Download table as:  ASCIITypeset image

The linearized equations are solved for an inner boundary condition: $\delta {v}_{r}=0$. As for the outer boundary condition, based on the numerical result of Takahashi & Yamada (2014), we set the perturbed flow that accretes onto the shock surface approximately as

$$\begin{eqnarray}&&\displaystyle \frac{\delta \rho }{{\rho }_{0}}=\sin ({\omega }_{\mathrm{up}}t+\varphi ),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\delta {v}_{r}}{{v}_{r0}}=-0.5\sin ({\omega }_{\mathrm{up}}t+\varphi ),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\delta \varepsilon }{{\varepsilon }_{0}}=\sin ({\omega }_{\mathrm{up}}t+\varphi ).\end{eqnarray} \tag{ 27 }$$

According to their results, the time evolutions of perturbed quantities at a fixed radius are not simple harmonic oscillations. If we focus, however, on the timescale of a few hundred milliseconds, which is relevant for the development of SASI and shock revival itself, the approximation by the sinusoidal function is reasonable, since the temporal variation of the perturbation at the shock front is dictated mostly by the structure of the unperturbed accretion flow and is rather insensitive to the initial fluctuation in the envelope. See Takahashi & Yamada (2014) for more details. It was found in the same paper that the typical frequency is in the range of 10–100 s⁻¹ and tends to become larger with l. In this study we regard the oscillation frequency as a free parameter and vary it in this range.

In this paper, we mainly investigate unstable eigenmodes and the effects of upstream fluctuations on them. We start with the intrinsic instabilities, by which we mean the unstable eigenmodes in the absence of the upstream perturbations. We then show the latter effects.

3.1. Intrinsic Instabilities

As shown in the next section, the growth/damping rate and oscillation frequency of eigenmodes are not changed by the upstream perturbations to the linear order. We can hence extract this information of the intrinsic modes by neglecting the upstream perturbation (the proof is given in the next section). It is then equivalent to the ordinary normal mode analysis of SASI and/or convection with Fourier transform, and it turns out that the results are qualitatively consistent with those of Yamasaki & Yamada (2007) although they used different background flows with a higher mass accretion rate.

We first mention the spherically symmetric mode (l = 0), which is shown in Figure 2 for various neutrino luminosities. Different modes are distinguished with colors. They are normally classified according to the number of radial nodes in the corresponding eigenfunctions: the mode with the smallest number of nodes is called the fundamental mode, while the others are referred to as the first, second, third, ..., overtones as the number of nodes increases (Yamasaki & Yamada 2007). Since we do not obtain eigenfunctions directly in our method, we refer to the mode with the lowest oscillation frequency as the fundamental mode and call other eigenmodes the first, second, and third overtones, and so on, in the ascending order of the oscillation frequency.

Zoom In Zoom Out Reset image size

As is apparent in the left panel of Figure 2, the background flow is stable for radial perturbations as long as the neutrino luminosity is low. It becomes unstable, however, once the luminosity reaches a threshold value, ${L}_{\nu }\sim 4.9\times {10}^{52}$ erg s⁻¹, where the first overtone becomes unstable. For much higher luminosities, the first overtone bifurcates into two branches and turns out to be nonoscillatory. Bifurcations are also seen in higher-l modes, and we discuss the physical interpretation later. The fundamental mode of l = 0 is always nonoscillatory irrespective of neutrino luminosities and was identified as a thermal mode by Yamasaki & Yamada (2007), which has no counterpart for higher-l modes and will disappear if we turn off the perturbation of the heating rate.⁷

Next, we consider nonspherically symmetric modes in Figures 3 and 4, where $l=1,\ldots ,6$ modes are shown. We note that the magnitudes of oscillation frequencies for different modes can interchange with one another as the neutrino luminosity increases. This was not observed in Yamasaki & Yamada (2007).⁸ As demonstrated clearly, the modes with $l=1,2$ become unstable first as the luminosity increases, and the growth rate is highest for the l = 1 mode as long as the neutrino luminosity is less than ${L}_{\nu }\sim 5\times {10}^{52}$ erg s⁻¹. These are SASI modes, which are likely to be driven by the advective-acoustic cycle, since the oscillation frequencies of these unstable modes follow more closely those of the advective-acoustic cycle, ${\omega }_{\mathrm{aac}}$, rather than those of the purely acoustic cycle, ${\omega }_{\mathrm{pac}}$, or of the mean Brunt–Väisälä frequencies, ${\bar{\omega }}_{\mathrm{BV}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the neutrino luminosity becomes higher, the fundamental mode and the first overtone in some cases bifurcate into two branches and become nonoscillatory, as does the first overtone with l = 0. The growth rate in one of the two branches increases rapidly thereafter, and higher-l modes become dominant as a result. To compare these bifurcated modes, we compared them for $0\leqslant l\leqslant 8$ in Figure 5. These bifurcated nonspherically symmetric modes may be interpreted as convective modes (Yamasaki & Yamada 2007). It is found that the value of the χ parameter at the point when one of the bifurcated modes first becomes unstable is ∼4 in our models, which seems to support the claim by Foglizzo et al. (2006).⁹ It is also interesting to point out that the bifurcation seems to occur when the oscillation frequency falls below the frequency of the advective-acoustic cycle, ${\omega }_{\mathrm{aac}}$, as seen in Figures 2–5.¹⁰

Zoom In Zoom Out Reset image size

Since the spherically symmetric mode becomes unstable at ${L}_{\nu }\sim 4.9\times {10}^{52}$ erg s⁻¹ and its growth rate is smaller than the convective modes, nonspherical modes are more important practically in the SN explosion. SASI is dominant in the low-luminosity regime, while convection overwhelms SASI in the high-L_ν regime. These results are in agreement qualitatively with what we have observed in many realistic simulations (e.g., Burrows et al. 2012; Hanke et al. 2013; Iwakami et al. 2014a, 2014b), as well as with the previous linear analysis (Yamasaki & Yamada 2007).

3.2. Upstream Perturbations

Now we proceed to the case in which we impose nonspherical perturbations also in the flow ahead of the shock surface. As proved in the next section, the upstream perturbations excite the intrinsic instabilities while not changing the eigenfrequencies. Put more specifically, the amplitude of an intrinsically unstable mode turns out to be given as a linear combination of the contributions from the upstream perturbations, initial fluctuations in the shocked region, and fluctuations at the inner boundary (see Appendices D and E).

In this paper, we pay attention to the first contribution and evaluate numerically the resultant amplitudes of the modes excited by these upstream fluctuations. They are given by Equation (108) or equivalently Equation (109) as a product of the residue at the corresponding pole of the Laplace-transformed shock radius for the impulsive perturbation and the functional value of the Laplace-transformed upstream perturbation at the pole position.

We consider both SASI and convective instabilities with $0\leqslant l\leqslant 8$, which have been obtained in the previous subsection, and discuss them separately below. We assume the upstream perturbation to have the functional forms given in Equations (25)–(27). We also change the frequency ${\omega }_{\mathrm{up}}$ and phase φ included in them.

3.2.1. Effects on SASI

We first study the excitation of intrinsically unstable SASI modes, which have ${\rm{\Omega }}\gt 0$ and $\omega \gt 0$, by sinusoidal upstream perturbations with $\varphi =0,\pi$/2 and ${\omega }_{\mathrm{up}}=1\mbox{--}{10}^{4}$ s⁻¹.

Figure 6 shows the results for ${L}_{\nu }=3\times {10}^{52}$ erg s⁻¹, where the absolute values of the excited amplitudes are plotted. The amplitudes are normalized by that of the upstream perturbation, and we note that the vertical scale differs from panel to panel. There are two unstable SASI modes with $l=1,2$ found for the neutrino luminosity. Similarly, Figure 7 presents the results for ${L}_{\nu }=4\times {10}^{52}$ erg s⁻¹, where another l = 2 mode becomes unstable, as well as an l = 3 mode. Figures 8 and 9 show the results for still higher luminosities, $5\times {10}^{52}$ erg s⁻¹ and $6\times {10}^{52}$ erg s⁻¹, respectively. We note that some convective modes also become unstable for these high luminosities, which are discussed separately later and not shown in the figures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a general trend of these results, we point out that the excited amplitude is larger for the cosine-type perturbation than for the sine-type if the frequency of the upstream perturbation is much smaller than the oscillation frequency of the eigenmode. As ${\omega }_{\mathrm{up}}$ gets larger, on the other hand, the amplitude of excitation for sine-type perturbations increases rapidly, and when ${\omega }_{\mathrm{up}}$ is comparable to the oscillation frequency of the eigenmode, the amplitudes for both types of perturbations are also similar to each other. For even larger ${\omega }_{\mathrm{up}}$ the order is reversed, with the excitation by the sine-type perturbations being dominant.

The magnitude of the excited amplitudes is in the range from 10⁻⁶ to 10. As shown in the figures, it peaks around ${\omega }_{\mathrm{up}}\sim {\omega }_{j}$ (the oscillation frequency of the jth mode) and becomes smaller as ${\omega }_{\mathrm{up}}$ goes away from ${\omega }_{j}$. At high luminosities, as the upstream frequency becomes higher, however, the amplitudes of all modes become comparable to one another, which is most clearly discernible in Figure 9.

3.2.2. Effects on Convection

The excitation of unstable convective modes by the upstream fluctuations is studied in the same way as for the SASI modes. Since their oscillation frequency vanishes by definition, we show in Figures 10 and 11 the results in the two-dimensional plane of the excited amplitude and growth rate. The former figure displays the results for ${L}_{\nu }=5\times {10}^{52}$ erg s⁻¹, while the latter is for ${L}_{\nu }=6\times {10}^{52}$ erg s⁻¹. We note that the convective modes with $1\leqslant l\leqslant 8$ are unstable for ${L}_{\nu }\gtrsim 5\times {10}^{52}$ erg s⁻¹, as indicated in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As in the SASI modes, we find that the amplitude of excitation is larger for the cosine-type perturbations if ${\omega }_{\mathrm{up}}$ is much smaller, while the excitation is dominated by the sine-type perturbations if ${\omega }_{\mathrm{up}}$ is large enough. The order change appears to occur at ${\omega }_{\mathrm{up}}\sim {{\rm{\Omega }}}_{j}$ (the growth rate of the jth mode).

The convective modes with smaller Ω are likely to be more strongly excited especially when ${\omega }_{\mathrm{up}}$ is very small. As the upstream frequency increases, the difference in the excited amplitudes among various modes tends to become smaller, as seen in the figures, with the modes with small Ω being suppressed more strongly.

3.2.3. Forced-oscillation Modes

When the sinusoidal perturbation is added ahead of the shock, there naturally occurs a mode that oscillates stably with the same frequency as shown analytically in the next section. This is analogous to the forced oscillation of a spring.

The amplitude of this mode is also obtained from Equation (104). The residue at the corresponding pole is evaluated by integrating the Laplace-transformed function of the perturbed shock radius around the pole. We show the results in Figure 12.

Zoom In Zoom Out Reset image size

For the cosine-type perturbations, the amplitude tends to increase as the upstream frequency decreases. In the case of the sine-type fluctuations, there appears a broad peak at ω_up ∼ 10–100 s⁻¹. As for the l-dependence, the amplitude becomes larger for smaller l when the luminosity is small, and hence SASI dominates convections (${L}_{\nu }=3,4\times {10}^{52}$ erg s⁻¹), while it is reversed for larger luminosities, where convections become dominant (${L}_{\nu }=5,6\times {10}^{52}$ erg s⁻¹), although there are some outliers.

Since the growth rate of the mode vanishes, it will eventually be dominated by unstable SASI or convections. For upstream fluctuations of finite amplitudes, which may not be small indeed (Takahashi & Yamada 2014), the forced-oscillation mode may not be negligible, either.

To understand the numerical results presented above, we give some analytical considerations, which can be summarized as follows: (1) the growth/damping rates and oscillation frequencies of intrinsic modes are unaffected by the presence of the upstream perturbations; (2) these intrinsic modes are excited with different amplitudes by the upstream perturbations. In the following subsections, we discuss them in turn.

4.1. Effects of Upstream Perturbations on the Eigenfrequencies of Unstable Modes

We here show that neither the growth rates nor the oscillation frequencies of intrinsic modes are altered by upstream fluctuations except for a special case.

The Laplace-transformed shock evolution is described by Equation (82), which is divided into three parts: the contribution from the upstream perturbations, Equation (86), the one from the initial perturbations in the shocked region, Equation (95), and the one from the inner boundary, Equation (96). Recalling that the growth rates and oscillation frequencies of eigenmodes are obtained from the location of the corresponding poles, we conclude that zero points of the denominator in these equations give the eigenvalues unless cancellation occurs between the denominator and numerator, which is not expected in general. It is also possible that the numerators have their own poles, which is discussed later. Since the denominator is dependent only on the background quantities and the inner boundary condition, the poles that originated from the denominator are unaffected by the presence of the upstream perturbation. This is of course true up to the linear order.

The poles of the numerator correspond to the forced-oscillation mode (see Section 3.2.3). For example, if a sinusoidally oscillating density perturbation $\delta \rho ({r}_{\mathrm{sh}},t)\propto \sin ({\omega }_{\mathrm{up}}t)$ exists just ahead of the shock front, then the Laplace-transformed shock radius, Equation (86), has a pole at $s=\pm {\omega }_{\mathrm{up}}$ (see Equation (18)) This corresponds to an "eigen" mode that oscillates stably at the same frequency as the original upstream perturbation. It is stressed that this mode is stable and has nothing to do with other intrinsic modes in general.

An exceptional case is the resonance, however. It occurs only when the upstream perturbation has coincidentally the same growth rate and oscillation frequency as one of the eigenvalues, i.e., when the upstream perturbation is given by

$$\begin{eqnarray}&&{\boldsymbol{z}}(t)\propto {e}^{{{\rm{\Omega }}}_{n}t}\sin ({\omega }_{n}t+\varphi ),\end{eqnarray} \tag{ 28 }$$

where ${{\rm{\Omega }}}_{n}$ and ${\omega }_{n}$ are the growth rate and oscillation frequency of the nth intrinsic eigenmode, respectively. Then the denominator of Equation (86) is factorized by ${(s-{{\rm{\Omega }}}_{n}\mp {\omega }_{n})}^{2}$, that is, $s={{\rm{\Omega }}}_{n}\pm i{\omega }_{n}$ becomes a double pole. Note that the inverse Laplace transform of a multipole is

$$\begin{eqnarray}&&{{ \mathcal L }}^{-1}\left[\displaystyle \frac{1}{{(s-{s}_{0})}^{\nu }}\right]=\displaystyle \frac{{t}^{\nu -1}{e}^{{s}_{0}t}}{{\rm{\Gamma }}(\nu )},\end{eqnarray} \tag{ 29 }$$

where $\nu \gt 0$ is a real number, s₀ a complex constant, and ${\rm{\Gamma }}(\nu )$ the gamma function. Applying the relation to the double pole, we find that the evolution of the shock radius for the eigenmode becomes proportional to $t\exp ({{\rm{\Omega }}}_{n}t)\sin ({\omega }_{n}t)$, i.e., it acquires an extra power of t. This is particularly important for ${{\rm{\Omega }}}_{n}=0$. We emphasize again, however, that the resonance occurs only when both the oscillation frequencies (and the growth rates if any) are identical for the upstream perturbation and one of the intrinsic modes, which may not be expected in general.

4.2. Amplitudes of Intrinsic Modes Excited by Upstream Perturbations

More important is the fact that the upstream perturbation, upon hitting the shock, excites the intrinsic modes even without initial fluctuations in the downstream or the shocked region.

As discussed in Appendix E, the excited amplitude of the jth mode is found to be roughly proportional to ${z}_{i}^{* }(s={s}_{j})(i=1,\ldots ,n)$, as shown in Equation (107), where ${s}_{j}={{\rm{\Omega }}}_{j}+i{\omega }_{j}$ denotes the position of the corresponding pole. More precisely speaking, the amplitude is given by the absolute value of a complex number, which is a sum of the products of ${z}_{i}^{* }(s={s}_{j})$ and the residue at $s={s}_{j}$ of ${{ \mathcal J }}_{i}^{* }$, where ${{ \mathcal J }}_{i}$ is the response of the shock radius to the impulsive perturbation of the ith component.

It will be informative to give $\mathrm{Res}{{ \mathcal J }}_{j}^{* }$, since it is one of the elements to determine the excited amplitudes as mentioned above. We present their absolute values for the impulsive perturbations to density, radial velocity, and internal energy separately in Figures 13–15. As seen in the figures, the response to the perturbation in the internal energy is a few orders of magnitude smaller than those for the other cases, which are comparable to each other. Hence, the density and/or velocity perturbations play more important roles in exciting intrinsic modes than the internal energy perturbation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Assuming that the upstream perturbation is a sinusoidal function in time, i.e., ${z}_{i}(t)\propto \sin ({\omega }_{\mathrm{up}}t+\varphi )$ and hence ${z}_{i}^{* }(s)\propto (s\sin \varphi +{\omega }_{\mathrm{up}}\cos \varphi )/({s}^{2}+{\omega }_{\mathrm{up}}^{2})$, we obtain

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})=\displaystyle \frac{({{\rm{\Omega }}}_{j}+i{\omega }_{j})\sin \varphi +{\omega }_{\mathrm{up}}\cos \varphi }{{({{\rm{\Omega }}}_{j}+i{\omega }_{j})}^{2}+{\omega }_{\mathrm{up}}^{2}}.\end{eqnarray} \tag{ 30 }$$

We show the values of Equation (30) in the plane of ${\omega }_{\mathrm{up}}$ and φ as a color contour in Figures 16–18. The former two figures illuminate the dependence on ${{\rm{\Omega }}}_{j}$ and ${\omega }_{j}$. Figure 18, on the other hand, compares SASI and convection, choosing the typical values of ${{\rm{\Omega }}}_{j}$ and ${\omega }_{j}$ for these cases: ${{\rm{\Omega }}}_{j}=50$ s⁻¹ and ${\omega }_{j}=300$ s⁻¹ for SASI in the upper panel and ${{\rm{\Omega }}}_{j}=150$ s⁻¹ and ${\omega }_{j}=0$ s⁻¹ for convection in the lower panel. As seen in the figure, ${z}_{i}^{* }$ takes maximum and minimum for the ${\omega }_{\mathrm{up}}$ values that are comparable to ${\omega }_{j}$ and is reduced to zero rather quickly as ${\omega }_{\mathrm{up}}$ departs from ${\omega }_{j}$, which is a common feature to SASI and convection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.1. Excitation of SASI

In order to explain the amplitude distributions of excited SASI modes in the previous section, we focus on the case with $\varphi =0$ or $\pi /2$. Paying attention to the fact that the growth rate ${\rm{\Omega }}\sim 10$ s⁻¹ of SASI is in general much smaller than its oscillation frequency $\omega \sim 100$ s⁻¹, we consider some limiting cases below.

We first take the high-frequency limit of ${\omega }_{\mathrm{up}}$: ${{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}\ll {\omega }_{j}/{\omega }_{\mathrm{up}}\ll 1$. Neglecting higher-order terms, we find for high-frequency sine-type/cosine-type perturbations (HFSPs/HFCPs), respectively, the following results:

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})\sim \left(i\displaystyle \frac{{\omega }_{j}}{{\omega }_{\mathrm{up}}}\sin \varphi +\cos \varphi \right){\omega }_{\mathrm{up}}^{-1}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}=\left\{\begin{array}{cc}{\omega }_{\mathrm{up}}^{-1} & \ (\varphi =0\,:\mathrm{HFSPs})\\ i\displaystyle \frac{{\omega }_{j}}{{\omega }_{\mathrm{up}}}{\omega }_{\mathrm{up}}^{-1} & \ (\varphi =\pi /2\,:\mathrm{HFCPs})\end{array}\right..\end{eqnarray} \tag{ 32 }$$

Note that HFCPs give much smaller amplitudes than HFSPs.

If we consider low-frequency upstream perturbations, on the other hand, ${\omega }_{j}/{\omega }_{\mathrm{up}}\gg {{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}$ and ${\omega }_{j}/{\omega }_{\mathrm{up}}\gg 1$, we obtain the following:

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})\sim -\left(i\sin \varphi +\displaystyle \frac{{\omega }_{\mathrm{up}}}{{\omega }_{j}}\cos \varphi \right){\omega }_{j}^{-1}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}=\left\{\begin{array}{cc}-\,\displaystyle \frac{{\omega }_{\mathrm{up}}}{{\omega }_{j}}{\omega }_{j}^{-1} & \ (\varphi =0\,:\mathrm{LFSPs})\\ -\,i{\omega }_{j}^{-1} & \ (\varphi =\pi /2\,:\mathrm{LFCPs})\end{array}\right.,\end{eqnarray} \tag{ 34 }$$

where LFSPs/LFCPs are the abbreviations of low-frequency sine-type/cosine-type perturbations. We note that the sine-type perturbation gives much smaller amplitudes than the cosine-type perturbation this time.

The last case we consider is the one in which the frequency of upstream perturbation is close to the oscillation frequency of SASI: ${\omega }_{j}/{\omega }_{\mathrm{up}}\sim 1\gg {{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}$. We introduce here an index, $\nu \gt 0$, to measure how close these two frequencies are to each other:

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{j}}{{\omega }_{\mathrm{up}}}=1+O\left({\left(\displaystyle \frac{{{\rm{\Omega }}}_{j}}{{\omega }_{\mathrm{up}}}\right)}^{\nu }\right).\end{eqnarray} \tag{ 35 }$$

For $\nu \gg 1$, ${\omega }_{\mathrm{up}}$ is very close to ${\omega }_{j}$, and we obtain

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})\sim \displaystyle \frac{1}{2}(\sin \varphi -i\cos \varphi ){{\rm{\Omega }}}_{j}^{-1}\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}=\left\{\begin{array}{cc}-\,\displaystyle \frac{i}{2}{{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =0\,:\mathrm{ExSFSPs})\\ \displaystyle \frac{1}{2}{{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =\pi /2\,:\mathrm{ExSFCPs})\end{array}\right.,\end{eqnarray} \tag{ 37 }$$

where the prefix ExSF means "extremely similar frequency." For $\nu \ll 1$, on the other hand, ${\omega }_{\mathrm{up}}$ is a bit more different from ${\omega }_{j}$, and we obtain

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})\sim {\left(\displaystyle \frac{{\omega }_{j}}{{{\rm{\Omega }}}_{j}}\right)}^{\nu }(i\sin \varphi +\cos \varphi ){\omega }_{j}^{-1}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}=\left\{\begin{array}{cc}{\left(\displaystyle \frac{{\omega }_{j}}{{{\rm{\Omega }}}_{j}}\right)}^{\nu }{\omega }_{j}^{-1} & \ (\varphi =0\,:\mathrm{SFSPs})\\ i{\left(\displaystyle \frac{{\omega }_{j}}{{{\rm{\Omega }}}_{j}}\right)}^{\nu }{\omega }_{j}^{-1} & \ (\varphi =\pi /2\,:\mathrm{SFCPs})\end{array}\right.,\end{eqnarray} \tag{ 39 }$$

where the prefix SF means "similar frequency."

From the comparison of the results given above, we can deduce which type of perturbations can drive what modes efficiently. In fact, we can easily get the following inequalities: ExSFCPs ∼ ExSFSPs  $\propto \,({\omega }_{j}/{{\rm{\Omega }}}_{j}){\omega }_{j}^{-1}\,\gg$  SFCPs ∼ SFSPs $\propto {({\omega }_{j}/{{\rm{\Omega }}}_{j})}^{\nu }{\omega }_{j}^{-1}\,\gtrsim$ LFCPs $\,\propto {\omega }_{j}^{-1}\,\gg$  LFSPs $\propto ({\omega }_{\mathrm{up}}^{L}/{\omega }_{j}){\omega }_{j}^{-1}\,\sim$ HFSPs  $\propto \,({\omega }_{j}/{\omega }_{\mathrm{up}}^{H}){\omega }_{j}^{-1}\,\gg$  HFCPs  $\propto \,{({\omega }_{j}/{\omega }_{\mathrm{up}}^{H})}^{2}{\omega }_{j}^{-1}$. We added here the superscripts L and H to emphasize that ${\omega }_{\mathrm{up}}^{L}\ll {\omega }_{j}\ll {\omega }_{\mathrm{up}}^{H}$. These inequalities indicate that the upstream perturbations with extremely similar frequencies are the most efficient in driving SASI. Even if the frequency is not extremely similar, $\nu \ll 1$, it is still effective. If ${\omega }_{\mathrm{up}}$ is substantially different from ${\omega }_{j}$, on the other hand, LFCPs are the most efficient in exciting SASI and HFCPs are the most inefficient. The sine-type perturbations are always in the middle, being independent of ${\omega }_{\mathrm{up}}$.

It makes sense that the upstream perturbations with the frequencies that are similar to that of an intrinsic mode are good at exciting the mode. As mentioned earlier, the resonance, which has an extra power dependence of time as $t\exp ({{\rm{\Omega }}}_{j}t)\sin ({\omega }_{j}t)$, occurs only if both the growth rate and oscillation frequency of upstream perturbations coincide with those of an intrinsic mode. This is practically impossible, and what occurs in reality is the large excitation amplitude of the intrinsic mode that has an oscillation frequency very close to that of the upstream perturbation.

If such a quasi-resonance condition is not satisfied, LFCPs are more effective because they are almost step function perturbations, which have broad spectra in frequency. In fact, the amplitude excited by LFCPs is independent of ${\omega }_{\mathrm{up}}$, since there is no characteristic timescale in the step function. On the other hand, LFSPs go to zero in the limit of ${\omega }_{\mathrm{up}}\to 0$, which is obviously ineffective to drive intrinsic modes. As a matter of fact, ${z}_{i}^{* }$ is proportional to ${\omega }_{\mathrm{up}}$, which, too, goes to zero in the same limit. It is also understandable that intrinsic modes with low oscillation frequencies are more strongly amplified by LFCPs and LFSPs.

High-frequency perturbations are inefficient in exciting intrinsic modes as is evident from the fact that the excited amplitudes go to zero in the limit of ${\omega }_{\mathrm{up}}\to \infty$. This may be interpreted as a consequence of the phase cancellation in rapid oscillation of the upstream perturbation.

4.2.2. Excitation of Convections

Next, we consider the excited amplitudes of convective modes. It is noted here again that they are characterized by the vanishing oscillation frequency, ${\omega }_{j}=0$. Neglecting higher-order terms in Equation (30) again, we obtain for the high-frequency limit, i.e., ${{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}\ll 1$,

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})\sim \left(\displaystyle \frac{{{\rm{\Omega }}}_{j}}{{\omega }_{\mathrm{up}}}\sin \varphi +\cos \varphi \right)\,{\omega }_{\mathrm{up}}^{-1}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}=\left\{\begin{array}{lc}{\omega }_{\mathrm{up}}^{-1} & \ (\varphi =0\,:\mathrm{HFSPs})\\ \displaystyle \frac{{{\rm{\Omega }}}_{j}}{{\omega }_{\mathrm{up}}}{\omega }_{\mathrm{up}}^{-1} & \ (\varphi =\pi /2\,:\mathrm{HFCPs})\end{array}\right..\end{eqnarray} \tag{ 41 }$$

For the low-frequency limit (${{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}\gg 1$), on the other hand, we find

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})\sim \left(\sin \varphi +\displaystyle \frac{{\omega }_{\mathrm{up}}}{{{\rm{\Omega }}}_{j}}\cos \varphi \right)\,{{\rm{\Omega }}}_{j}^{-1}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}=\left\{\begin{array}{cc}\displaystyle \frac{{\omega }_{\mathrm{up}}}{{{\rm{\Omega }}}_{j}}{{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =0\,:\mathrm{LFSPs})\\ {{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =\pi /2\,:\mathrm{LFCPs})\end{array}\right..\end{eqnarray} \tag{ 43 }$$

For the last case of ${{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}\sim 1$, we obtain

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})\sim \displaystyle \frac{1}{2}(\sin \varphi +\cos \varphi ){{\rm{\Omega }}}_{j}^{-1}\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}=\left\{\begin{array}{lc}\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =0\,:\mathrm{SFSPs})\\ \displaystyle \frac{1}{2}{{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =\pi /2\,:\mathrm{SFCPs})\end{array}\right..\end{eqnarray} \tag{ 45 }$$

Taking the derivative of ${z}_{i}^{* }$ with respect to ${\omega }_{\mathrm{up}}$, we see that the maximum value of ${z}_{i}^{* }$ is obtained at ${\omega }_{\mathrm{up}}/{{\rm{\Omega }}}_{j}=\cos \varphi /(1+\sin \varphi )$ as

$$\begin{eqnarray}&&{z}_{i}^{* }({s}_{j})=\displaystyle \frac{(1+\sin \varphi )}{2{{\rm{\Omega }}}_{j}}\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}=\left\{\begin{array}{lc}\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =0\,\Longleftrightarrow \,{\omega }_{\mathrm{up}}={{\rm{\Omega }}}_{j}\,:\mathrm{SFSPs})\\ {{\rm{\Omega }}}_{j}^{-1} & \ (\varphi =\pi /2\,\Longleftrightarrow \,{\omega }_{\mathrm{up}}=0\,:\mathrm{LFCPs})\end{array}\right..\end{eqnarray} \tag{ 47 }$$

It is interesting that ${\omega }_{\mathrm{up}}={{\rm{\Omega }}}_{j}$ for the sine-type perturbations, since the former is an oscillation frequency whereas the latter is a growth rate. For the cosine-type perturbations, on the other hand, the maximum is obtained in the low-frequency limit: ${\omega }_{\mathrm{up}}=0$.

From the comparison of the results given above, we obtain the following inequalities: SFCP ∼ SFSPs ∼ LFCPs $\propto \,{{\rm{\Omega }}}_{j}^{-1}\,\gg$  LFSPs  $\propto \,({\omega }_{\mathrm{up}}^{L}/{{\rm{\Omega }}}_{j}){{\rm{\Omega }}}_{j}^{-1}\,\sim$ HFSPs $\propto \,({{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}^{H}){{\rm{\Omega }}}_{j}^{-1}\,\gg$ HFCPs $\propto \,{({{\rm{\Omega }}}_{j}/{\omega }_{\mathrm{up}}^{H})}^{2}{{\rm{\Omega }}}_{j}^{-1}$. It is hence concluded that LFCPs, SFCPs, and SFSPs are more efficient in the excitation of convections as for SASI. The excited amplitudes are inversely proportional to ${{\rm{\Omega }}}_{j}$ and independent of ${\omega }_{\mathrm{up}}$, which is similar to the SASI case again.

4.3. General Perturbations

4.4. Implications for the Shock Dynamics

We investigated the links between the fluid instabilities in the iron core of CCSNe and the nonspherical fluctuations in the accreting envelopes, which are expected to be generated by violent convections in the Si/O layers. We formulated the problem as an initial boundary value problem and employed Laplace transform to find the intrinsic modes such as SASI and convection and calculate their excitations by the upstream perturbations.

We first sought the intrinsic modes, especially unstable ones, which are eigenmodes in the absence of the upstream perturbations, for background flows with different model parameters. The results are consistent with those of the previous normal mode analyses (Foglizzo et al. 2006; Yamasaki & Yamada 2007; Guilet & Foglizzo 2012): SASI, which seems to be driven by the advective-acoustic cycle, prevails in low neutrino luminosities, whereas convection becomes dominant as the luminosity increases. Then we investigated the amplitudes of the modes that are excited by the perturbations in the matter accreting on the shock. Based on our previous work on the amplification of fluctuations during the accretion (Takahashi & Yamada 2014), we approximated the perturbations just ahead of the shock by the functional form given in Equations (25)–(27). Changing the parameters, i.e., phase and frequencies of the perturbations, rather arbitrarily, we systematically studied the couplings between the upstream perturbations and the intrinsic modes for different background flows. We also gave analytic expressions to the excited amplitudes in some limiting cases.

We showed analytically that the resonance, in which the growth rate of an unstable mode acquires an extra power of t, occurs only when both of their oscillation frequencies and growth rates coincide exactly between the intrinsic mode and the upstream perturbation. Hence, it does not happen practically. What occurs actually is that, as we demonstrated both analytically and numerically, the excited amplitude becomes larger when the upstream frequency is close either to one of the oscillation frequencies for SASI or to one of the growth rates for convection. If these frequencies are not very close to each other, the excitation efficiency declines rather quickly. It is also demonstrated that the discontinuous perturbation given as a step function of time can excite various modes thanks to its broad spectrum in frequency. Fortunately, the upstream perturbations acquire temporal variations during accretion, whose timescales are not much different from those of SASI or convection, as reported in Takahashi & Yamada (2014).

The magnitude of the excitations may not be small. In fact, the violent convections in O and Si burnings are supposed to generate fluctuations at the several percent level (Bazan & Arnett 1998; Asida & Arnett 2000; Meakin & Arnett 2006, 2007; Arnett & Meakin 2011; Chatzopoulos et al. 2014; Müller et al. 2016), which will then be amplified during accretion by a factor of a few to 10 (Takahashi & Yamada 2014). They in turn may excite some intrinsic modes with amplitudes that are larger by another factor of 1–10. We may hence expect that the shock radius can be perturbed at several tens of percent initially when the upstream fluctuations hit the shock. Of course, our linear analysis will not be applicable to such large-amplitude fluctuations. It should be also noted that even if the excited amplitude is not so large, the upstream perturbation shortens the time it takes the fluid instabilities to grow to nonlinear phases. All these findings seem to be consistent with what was observed in the recent numerical simulations (Couch & Ott 2013, 2015; Fernández et al. 2014; Müller & Janka 2014; Couch et al. 2015). The latest study by Abdikamalov et al. (2016) also found that the turbulent kinetic energy can be amplified by the passage through the shock wave, although their formulation cannot treat the shock instabilities such as SASI and convection.

Although the turbulent stellar structures are booming now, other multidimensional effects such as magnetic fields and stellar rotation may also play an important role in the links between the perturbations in the envelope and the development of instabilities in the core, and should be studied next.

K.T. thanks H. Nagakura for discussions on the dynamics of CCSNe. We acknowledge the center for Computational Astrophysics, CfCA, at the National Astronomical Observatory of Japan for the use of the XC30 and the general common-use computer system. This work is partially supported by a Research Fellowship for Young Scientists from the Japan Society for the Promotion of Science, as well as by the Grants-in-Aid for Scientic Research (A) (Nos. 24244036, 24740165) and the Grants-in-Aid for Scientic Research on Innovative Areas, New Development in Astrophysics through Multi-messenger Observations of Gravitational Wave Sources (No. 24103006),

In this appendix, we give the matrices and vectors that appeared in Section 2.

The vector that describes perturbed states is defined by Equation (9). In this order of the components, the matrices in the basic equations in the form of

$$\begin{eqnarray}&&M\displaystyle \frac{\partial {\boldsymbol{y}}}{\partial t}+A^{\prime} \displaystyle \frac{\partial {\boldsymbol{y}}}{\partial r}+B^{\prime} {\boldsymbol{y}}=0,\end{eqnarray} \tag{ 48 }$$

are given as follows:

$$\begin{eqnarray}M(r)=\left(\begin{array}{cccccc}1\ & 0\ & 0\ & 0\ & 0\ & 0\ \\ 0\ & 1\ & 0\ & 0\ & 0\ & 0\ \\ 0\ & 0\ & 1\ & 0\ & 0\ & 0\ \\ -\,\displaystyle \frac{p}{\rho {v}_{r}^{2}}\ & 0\ & 0\ & \displaystyle \frac{\varepsilon }{{v}_{r}^{2}}\ & 0\ & 0\ \\ 0\ & 0\ & 0\ & 0\ & 1\ & 0\ \\ 0\ & 0\ & 0\ & 0\ & 0\ & 1\ \end{array}\right),\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}A^{\prime} (r)=\left(\begin{array}{cccccc}{v}_{r}\ & {v}_{r}\ & 0\ & 0\ & 0\ & 0\ \\ \displaystyle \frac{1}{{v}_{r}}\displaystyle \frac{\partial p}{\partial \rho }\ & {v}_{r}\ & 0\ & \displaystyle \frac{\varepsilon }{\rho {v}_{r}}\displaystyle \frac{\partial p}{\partial \varepsilon }\ & \displaystyle \frac{{Y}_{e}}{\rho {v}_{r}}\displaystyle \frac{\partial p}{\partial {Y}_{e}}\ & 0\ \\ 0\ & 0\ & {v}_{r}\ & 0\ & 0\ & 0\ \\ -\,\displaystyle \frac{p}{\rho {v}_{r}}\ & 0\ & 0\ & \displaystyle \frac{\varepsilon }{{v}_{r}}\ & 0\ & 0\ \\ 0\ & 0\ & 0\ & 0\ & {v}_{r}\ & 0\ \\ 0\ & 0\ & 0\ & 0\ & 0\ & {v}_{r}\ \end{array}\right),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}B^{\prime} (r)=\left(\begin{array}{cccccc}0\ & 0\ & -{v}_{r}\displaystyle \frac{l(l+1)}{r} & 0\ & 0\ & 0\ \\ B{^{\prime} }_{11} & 2\displaystyle \frac{{{dv}}_{r}}{{dr}}\ & 0 & \displaystyle \frac{1}{\rho {v}_{r}}\displaystyle \frac{d}{{dr}}\left(\varepsilon \displaystyle \frac{\partial p}{\partial \varepsilon }\right)\ & \displaystyle \frac{1}{\rho {v}_{r}}\displaystyle \frac{d}{{dr}}\left({Y}_{e}\displaystyle \frac{\partial p}{\partial {Y}_{e}}\right) & \,0\\ \displaystyle \frac{1}{{{rv}}_{r}}\displaystyle \frac{\partial p}{\partial \rho }\ & 0\ & \displaystyle \frac{{v}_{r}}{r}+\displaystyle \frac{{{dv}}_{r}}{{dr}} & \displaystyle \frac{\varepsilon }{r\rho {v}_{r}}\displaystyle \frac{\partial p}{\partial \varepsilon }\ & \displaystyle \frac{{Y}_{e}}{r\rho {v}_{r}}\displaystyle \frac{\partial p}{\partial {Y}_{e}}\ & 0\ \\ B{^{\prime} }_{41}\ & \displaystyle \frac{1}{{v}_{r}}\left(\displaystyle \frac{d\varepsilon }{{dr}}-\displaystyle \frac{p}{{\rho }^{2}}\displaystyle \frac{d\rho }{{dr}}\right)\ & 0 & B{^{\prime} }_{44} & B{^{\prime} }_{45}\ & 0\ \\ \displaystyle \frac{{v}_{r}}{{Y}_{e}}\displaystyle \frac{{{dY}}_{e}}{{dr}}-\displaystyle \frac{{m}_{b}}{{Y}_{e}}\displaystyle \frac{\partial \lambda }{\partial \rho }\ & \displaystyle \frac{{v}_{r}}{{Y}_{e}}\displaystyle \frac{{{dY}}_{e}}{{dr}}\ & 0 & -\displaystyle \frac{{m}_{b}\varepsilon }{\rho {Y}_{e}}\displaystyle \frac{\partial \lambda }{\partial \varepsilon }\ & \displaystyle \frac{{v}_{r}}{{Y}_{e}}\displaystyle \frac{{{dY}}_{e}}{{dr}}-\displaystyle \frac{{m}_{b}}{\rho }\displaystyle \frac{\partial \lambda }{\partial {Y}_{e}}\ & 0\ \\ 0\ & 0 & \,0\ & 0\ & 0\ & \displaystyle \frac{{v}_{r}}{r}+\displaystyle \frac{{{dv}}_{r}}{{dr}}\ \end{array}\right),\end{eqnarray} \tag{ 51 }$$

with

$$\begin{eqnarray}&&B{^{\prime} }_{11}=\displaystyle \frac{1}{\rho {v}_{r}}\displaystyle \frac{d}{{dr}}\left(\rho \displaystyle \frac{\partial p}{\partial \rho }\right)-\displaystyle \frac{1}{\rho {v}_{r}}\displaystyle \frac{{dp}}{{dr}},\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&B{^{\prime} }_{41}=\displaystyle \frac{1}{{v}_{r}}\left(\displaystyle \frac{p}{{\rho }^{2}}\displaystyle \frac{d\rho }{{dr}}-\displaystyle \frac{1}{\rho }\displaystyle \frac{d\rho }{{dr}}\displaystyle \frac{\partial p}{\partial \rho }-\displaystyle \frac{\rho }{{v}_{r}}\displaystyle \frac{\partial q}{\partial \rho }\right),\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&B{^{\prime} }_{44}=\displaystyle \frac{1}{{v}_{r}}\left(\displaystyle \frac{d\varepsilon }{{dr}}-\displaystyle \frac{\varepsilon }{{\rho }^{2}}\displaystyle \frac{d\rho }{{dr}}\displaystyle \frac{\partial p}{\partial \varepsilon }-\displaystyle \frac{\varepsilon }{{v}_{r}}\displaystyle \frac{\partial q}{\partial \varepsilon }\right),\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&B{^{\prime} }_{45}=-\displaystyle \frac{1}{{v}_{r}}\left(\displaystyle \frac{{Y}_{e}}{{\rho }^{2}}\displaystyle \frac{d\rho }{{dr}}\displaystyle \frac{\partial p}{\partial {Y}_{e}}+\displaystyle \frac{{Y}_{e}}{{v}_{r}}\displaystyle \frac{\partial q}{\partial {Y}_{e}}\right).\end{eqnarray} \tag{ 55 }$$

Note that we took $\rho ,\varepsilon$, and Y_e as independent thermodynamic quantities in the above expression, and we did not express the fixed variables for notational simplicity when we take their partial derivatives. We obtain the basic equations in the form of Equation (8) by defining $A:= -A{^{\prime} }^{-1}M$ and $B:= -A{^{\prime} }^{-1}B$ and multiplying Equation (48) with ${A}^{-1}$ from the left.

The linearized Rankine–Hugoniot relations are schematically given by

$$\begin{eqnarray}{P}^{(d)}{\boldsymbol{y}}({r}_{\mathrm{sh}},t) & = & {P}^{(u)}{\boldsymbol{z}}(t)+\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}{\boldsymbol{c}}^{\prime} \\ & & +\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}{\boldsymbol{d}}^{\prime} ,\end{eqnarray} \tag{ 56 }$$

where ${\boldsymbol{z}}(t)$ is the vector that describes the perturbed quantities just in front of the shock surface. The matrices P and vectors ${\boldsymbol{c}}^{\prime}$ and ${\boldsymbol{d}}^{\prime}$ in the equation have the following forms:

$$\begin{eqnarray}P=\left(\begin{array}{cccccc}1\ & 1\ & 0\ & 0\ & 0\ & 0\ \\ {v}_{r}+\displaystyle \frac{1}{{v}_{r}}\displaystyle \frac{\partial p}{\partial \rho }\ & 2{v}_{r}\ & 0\ & \displaystyle \frac{\varepsilon }{j}\displaystyle \frac{\partial p}{\partial \varepsilon }\ & \displaystyle \frac{{Y}_{e}}{j}\displaystyle \frac{\partial p}{\partial {Y}_{e}}\ & 0\ \\ 0\ & 0\ & {v}_{r}\ & 0\ & 0\ & 0\ \\ \displaystyle \frac{E}{\rho }+\displaystyle \frac{\partial p}{\partial \rho }\ & \displaystyle \frac{E+p+\rho {v}_{r}^{2}}{\rho }\ & 0\ & \varepsilon +\displaystyle \frac{\varepsilon }{\rho }\displaystyle \frac{\partial p}{\partial \varepsilon }\ & \displaystyle \frac{{Y}_{e}}{\rho }\displaystyle \frac{\partial p}{\partial {Y}_{e}}\ & 0\ \\ 1\ & 1\ & 0\ & 0\ & 1\ & 0\ \\ 0\ & 0\ & 0\ & 0\ & 0\ & {v}_{r}\ \end{array}\right),\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&{{\boldsymbol{c}}}^{\prime }={\left(\displaystyle \frac{{r}_{\mathrm{sh}}[[\rho ]]}{j},0,0,\displaystyle \frac{{r}_{\mathrm{sh}}[[E]]}{j},\displaystyle \frac{{r}_{\mathrm{sh}}[[\rho ]]}{j},0\right)}^{T},\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}{{\boldsymbol{d}}}^{\prime } & = & \left(0,2[[{v}_{r}]]+\displaystyle \frac{{GM}[[\rho ]]}{{r}_{\mathrm{sh}}j},\displaystyle \frac{[[p]]}{j},-\displaystyle \frac{{r}_{\mathrm{sh}}[[\rho q]]}{j}\right.\\ & & {\left.-\displaystyle \frac{{r}_{\mathrm{sh}}{m}_{b}[[\lambda ]]}{{{jY}}_{e}},0\right)}^{T},\end{eqnarray} \tag{ 59 }$$

and the superscripts $({\rm{u}})$ and $({\rm{d}})$ of P mean that it is evaluated with the background quantities above and below the shock wave, respectively. The mass flux is denoted by $j:= {\rho }_{0}^{({\rm{u}})}{v}_{r0}^{({\rm{u}})}\ ={\rho }_{0}^{({\rm{d}})}{v}_{r0}^{({\rm{d}})}$, and the bracketed symbol, $[[X]]:= {X}^{({\rm{d}})}-{X}^{({\rm{u}})}$, is a jump of a quantity X across the shock.

Defining $R:= {P}^{-1({\rm{d}})}{P}^{({\rm{u}})}$, ${\boldsymbol{c}}:= {P}^{-1({\rm{d}})}{{\boldsymbol{c}}}^{\prime }$, and ${\boldsymbol{d}}:= {P}^{-1({\rm{d}})}{{\boldsymbol{d}}}^{\prime }$ and multiplying Equation (56) by ${P}^{-1({\rm{d}})}$ from the left, we obtain the linearized Rankine–Hugoniot relations in the form of Equation (11). Then, Laplace-transforming this equation with respect to t, we obtain Equation (14).

We show here that the boundary value problem given in Equations (13)–(15) is easily solved, provided that the function f* that gives the inner boundary condition in Equation (15) is linearly dependent on ${{\boldsymbol{y}}}^{* }$. In fact, we can then obtain the solution, $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$, for a given s by integrating the ordinary differential Equation (13) twice, as discussed below.

The idea is quite simple: this is a special case of the Newton–Raphson method, in which an approximate solution ${x}^{\star }$ to $F({x}^{\star })=0$ is improved by solving the linearized equation

$$\begin{eqnarray}&&F({x}^{\star }+\delta x)\sim F({x}^{\star })+J({x}^{\star })\delta x=0,\end{eqnarray} \tag{ 60 }$$

where J is the Jacobian, ${dF}/{dx}$. In fact, F and x can be replaced with f^* and $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$, respectively, in our problem. If $F\propto x$, the linearized equation is identical to the original equation, and hence the improved approximation given by ${x}_{\mathrm{sol}}={x}^{\star }+\delta x\ ={x}^{\star }-{J}^{-1}F$ is actually the exact solution. In the following, we show that our problem is indeed linear and obtain the Jacobian matrix.

The linearity is almost obvious since we are considering a linearized system: the inner boundary condition takes in general the form given in Equation (81), and ${{\boldsymbol{y}}}^{* }({r}_{{\nu }_{e}},s)$ is proportional to $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$, from which the linear dependence of the function f^* on $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$ follows:

$$\begin{eqnarray}{f}^{* } & = & {{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)\left[(s{\boldsymbol{c}}+{\boldsymbol{d}})\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }(s)}{{r}_{\mathrm{sh}}}+R{{\boldsymbol{z}}}^{* }(s)\right]\\ & & -{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\boldsymbol{h}}}}^{* }[{{\boldsymbol{y}}}_{0}](s)+{b}^{* }(s)=0.\end{eqnarray} \tag{ 61 }$$

The Jacobian matrix is then given as

$$\begin{eqnarray}&&J=\displaystyle \frac{{{df}}^{* }}{d(\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}})}={{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)(s{\boldsymbol{c}}+{\boldsymbol{d}}),\end{eqnarray} \tag{ 62 }$$

where ${{\boldsymbol{a}}}^{* }(s)$ is known once the function f* is given; the other factor, $s{\tilde{{\rm{\Lambda }}}}^{* }(s)(s{\boldsymbol{c}}+{\boldsymbol{d}})$, can be calculated by integrating the ordinary differential equation for the modified outer boundary condition:

$$\begin{eqnarray}&&{{\boldsymbol{y}}}^{* }({r}_{\mathrm{sh}},s)=s{\boldsymbol{c}}+{\boldsymbol{d}}.\end{eqnarray} \tag{ 63 }$$

We thus come to the conclusion that we can obtain the solution of the boundary value problem for a given s by integrating the ordinary differential Equation (13) twice: once for the outer boundary condition Equation (14) with an arbitrary guess for $\delta {r}_{\mathrm{sh}}/{r}_{\mathrm{sh}}$, and the other time for the modified boundary condition, Equation (63).

Inverse Laplace transform is given by the Fourier–Mellin formula (Schiff 1999):

$$\begin{eqnarray}&&f(t)=\displaystyle \frac{1}{2\pi i}{\int }_{L}{e}^{{ts}}{ \mathcal L }[f](s){ds},\end{eqnarray} \tag{ 64 }$$

where the integral path, L, is a straight line that is parallel to the imaginary axis and lies to the right of all the poles of ${f}^{* }(s)={ \mathcal L }[f](s)$. We do not know a priori where those poles are, however, and hence it may happen that the path sits to the left of some poles. We show below what happens if we employ such a wrong path.

Let us take a function $f(t)=\exp ({\rm{\Omega }}t)\sin (\omega t+\phi )$ as an example. Its Laplace transform is

$$\begin{eqnarray}&&{f}^{* }(s)=\displaystyle \frac{1}{2i}\left(\displaystyle \frac{{e}^{i\phi }}{s-{\rm{\Omega }}-i\omega }-\displaystyle \frac{{e}^{-i\phi }}{s-{\rm{\Omega }}+i\omega }\right),\end{eqnarray} \tag{ 65 }$$

for which the poles are located at $s={\rm{\Omega }}+\pm i\omega$. For appropriate paths with $\mathrm{Re}\ s\gt {\rm{\Omega }}$, the integral, Equation (64), correctly reproduces the original function as

$$\begin{eqnarray}\left\{\begin{array}{cc}{e}^{{\rm{\Omega }}t}\sin (\omega t+\phi ) & (t\gt 0),\\ \displaystyle \frac{1}{2}\sin \phi & (t=0),\\ 0 & (t\lt 0).\end{array}\right.\end{eqnarray} \tag{ 66 }$$

Here the point is that the inversed function vanishes at $t\lt 0$, which is guaranteed by the fact that the integral path is running to the right of all poles.

If, on the other hand, the path is inappropriate with either $\mathrm{Re}\ s\lt {\rm{\Omega }}$ or $\mathrm{Re}\ s={\rm{\Omega }}$, the integral gives incorrect results. In the former case, for example, the resultant function is

$$\begin{eqnarray}\left\{\begin{array}{cc}0 & (t\gt 0),\\ -\,\displaystyle \frac{1}{2}\sin \phi & (t=0),\\ -\,{e}^{{\rm{\Omega }}t}\sin (\omega t) & (t\lt 0).\end{array}\right.\end{eqnarray} \tag{ 67 }$$

The latter case, on the other hand, gives

$$\begin{eqnarray}\left\{\begin{array}{cc}\displaystyle \frac{1}{2}{e}^{{\rm{\Omega }}t}\sin (\omega t) & (t\gt 0),\\ 0 & (t=0),\\ -\,\displaystyle \frac{1}{2}{e}^{{\rm{\Omega }}t}\sin (\omega t) & (t\lt 0).\end{array}\right.\end{eqnarray} \tag{ 68 }$$

The important thing is that the integral gives nonzero values for $t\lt 0$ in both cases.

In Figure 19, we illustrate the situation for the following function, linear combination of two growing sinusoidal functions, which is somewhat closer to the actual problem we considered in this paper:

$$\begin{eqnarray}&&f(t)={e}^{10t}\sin (100t)+{e}^{50t}\sin (500t).\end{eqnarray} \tag{ 69 }$$

Its Laplace transform is given as

$$\begin{eqnarray}&&{f}^{* }(s)=\displaystyle \frac{100}{{(s-10)}^{2}+{100}^{2}}+\displaystyle \frac{500}{{(s-50)}^{2}+{500}^{2}},\end{eqnarray} \tag{ 70 }$$

which has four single poles, as shown in the figure (marked as red crosses). We performed the integral, Equation (64), for ${f}^{* }(s)$ along three different paths with $\mathrm{Re}\ s=-10$, $\mathrm{Re}\ s=30$, and $\mathrm{Re}\ s=70$. The results are presented in the same figure. Note that only the last one is appropriate for the inverse Laplace transform, which is indeed vindicated in the figure, with the functional value vanishing at $t\lt 0$. The integrals along the other two paths return the following results:

$$\begin{eqnarray}&&-[1-\theta (t)][{e}^{10t}\sin (100t)+{e}^{50t}\sin (500t)](\mathrm{Re}\ s=-10),\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}&&{e}^{10t}\sin (100t)\theta (t)-[1-\theta (t)]{e}^{50t}\sin (500t)(\mathrm{Re}\ s=30),\end{eqnarray} \tag{ 72 }$$

where $\theta (t)$ is the Heaviside function. They are evidently wrong, having nonvanishing values at $t\lt 0$ in particular.

Zoom In Zoom Out Reset image size

The fact that the inversed function should vanish at $t\lt 0$ can be employed to judge if we have missed the pole with the greatest real part, i.e., the one corresponding to the most unstable mode.

In this section, we assume that the perturbation in the shocked region is initially nonvanishing in general. Then an additional term shows up in the Laplace-transformed linearized equations as follows:

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{y}}}^{* }}{{dr}}=({sA}+B){{\boldsymbol{y}}}^{* }-A{{\boldsymbol{y}}}_{0}(r),\end{eqnarray} \tag{ 73 }$$

where ${{\boldsymbol{y}}}_{0}(r):= {\boldsymbol{y}}(r,t=0)$ denotes the initial perturbation. The outer and inner boundary conditions are unchanged and given as Equations (14) and (15), respectively.

Using the path-ordering operator, ${ \mathcal P }$, (Peskin & Schroeder 1995), the solution of Equations (73) with the outer boundary condition, Equation (14), is symbolically given as

$$\begin{eqnarray}&&{{\boldsymbol{y}}}^{* }(r,s)\\ &&\,\,={ \mathcal P }\left[\exp \left({\displaystyle \int }_{{r}_{\mathrm{sh}}}^{r}{dr}^{\prime} ({sA}+B)\right)\right]\{{{\boldsymbol{y}}}^{* }({r}_{\mathrm{sh}},s)\\ &&\quad \left.-{\displaystyle \int }_{{r}_{\mathrm{sh}}}^{r}{dr}^{\prime} { \mathcal P }\left[\exp \left(-{\displaystyle \int }_{{r}_{\mathrm{sh}}}^{r^{\prime} }{dr}^{\prime\prime} ({sA}+B)\right)\right]A(r^{\prime} ){{\boldsymbol{y}}}_{0}(r^{\prime} )\right\}\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}&&={{\rm{\Lambda }}}^{* }(r,s){{\boldsymbol{y}}}^{* }({r}_{\mathrm{sh}},s)-{{\rm{\Lambda }}}^{* }(r,s){\int }_{{r}_{\mathrm{sh}}}^{r}{dr}^{\prime} {{\rm{\Lambda }}}^{* -1}(r^{\prime} ,s)A(r^{\prime} ){{\boldsymbol{y}}}_{0}(r^{\prime} )\end{eqnarray} \tag{ 75 }$$

$$\begin{eqnarray}&&={{\rm{\Lambda }}}^{* }(r,s){{\boldsymbol{y}}}^{* }({r}_{\mathrm{sh}},s)-{{\boldsymbol{h}}}^{* }[{{\boldsymbol{y}}}_{0}](r,s),\end{eqnarray} \tag{ 76 }$$

where we defined an n × n matrix, ${{\rm{\Lambda }}}^{* }$, with n being the number of components in ${\boldsymbol{y}}$, and a vectorial functional, ${{\boldsymbol{h}}}^{* }$, of ${{\boldsymbol{y}}}_{0}$ as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}^{* }(r,s):= { \mathcal P }\left[\exp \left({\int }_{{r}_{\mathrm{sh}}}^{r}{dr}^{\prime} ({sA}+B)\right)\right],\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}&&{{\boldsymbol{h}}}^{* }[{{\boldsymbol{y}}}_{0}](r,s):= {{\rm{\Lambda }}}^{* }(r,s){\int }_{{r}_{\mathrm{sh}}}^{r}{dr}^{\prime} {{\rm{\Lambda }}}^{* -1}(r^{\prime} ,s)A(r^{\prime} ){{\boldsymbol{y}}}_{0}(r^{\prime} ).\end{eqnarray} \tag{ 78 }$$

Evaluating Equation (76) at the neutrino sphere and using the outer boundary condition, Equation (14), we obtain the following formal solution:

$$\begin{eqnarray}{{\boldsymbol{y}}}^{* }({r}_{{\nu }_{e}},s) & = & {{\rm{\Lambda }}}^{* }({r}_{{\nu }_{e}},s)\left[(s{\boldsymbol{c}}+{\boldsymbol{d}})\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }(s)}{{r}_{\mathrm{sh}}}+R{{\boldsymbol{z}}}^{* }(s)\right]\\ & & -{{\boldsymbol{h}}}^{* }[{{\boldsymbol{y}}}_{0}]({r}_{{\nu }_{e}},s)\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&=: {\tilde{{\rm{\Lambda }}}}^{* }(s)\left[(s{\boldsymbol{c}}+{\boldsymbol{d}})\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }(s)}{{r}_{\mathrm{sh}}}+R{{\boldsymbol{z}}}^{* }(s)\right]-{\tilde{{\boldsymbol{h}}}}^{* }[{{\boldsymbol{y}}}_{0}](s).\end{eqnarray} \tag{ 80 }$$

Here and hereafter we denote the functions ${{\rm{\Lambda }}}^{* }$ and ${{\boldsymbol{h}}}^{* }$ that are evaluated at the inner boundary as ${\tilde{{\rm{\Lambda }}}}^{* }(s)$ and ${\tilde{{\boldsymbol{h}}}}^{* }[{{\boldsymbol{y}}}_{0}](s)$ for notational simplicity.

The inner boundary condition, Equation (15), should be linear with respect to ${y}_{i}^{* }$, since we are considering linear perturbations. It hence takes the following form in general:

$$\begin{eqnarray}&&{{\boldsymbol{a}}}^{* }(s)\cdot {{\boldsymbol{y}}}^{* }({r}_{{\nu }_{e}},s)+{b}^{* }(s)=0,\end{eqnarray} \tag{ 81 }$$

where ${{\boldsymbol{a}}}^{* }$ and b* are some functions of s. Then the Laplace-transformed shock radius can be formally solved from Equations (80) and (81) as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}=-\displaystyle \frac{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)R{{\boldsymbol{z}}}^{* }(s)+{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\boldsymbol{h}}}}^{* }[{{\boldsymbol{y}}}_{0}](s)-{b}^{* }(s)}{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)(s{\boldsymbol{c}}+{\boldsymbol{d}})}\end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray}&&=: {{ \mathcal I }}^{* }[{{\boldsymbol{y}}}_{0}]({{\boldsymbol{z}}}^{* },s).\end{eqnarray} \tag{ 83 }$$

Let us consider the case with ${{\boldsymbol{z}}}^{* }=0$, i.e., there is no upstream perturbation:

$$\begin{eqnarray}&&{{ \mathcal I }}^{* }[{{\boldsymbol{y}}}_{0}](0,s)=-\displaystyle \frac{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\boldsymbol{h}}}}^{* }[{{\boldsymbol{y}}}_{0}](s)-{b}^{* }(s)}{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)(s{\boldsymbol{c}}+{\boldsymbol{d}})}.\end{eqnarray} \tag{ 84 }$$

The inverse transform of this function describes the time evolution of the perturbed shock radius that is induced by the initial perturbation given in the shocked region, ${{\boldsymbol{y}}}_{0}(r)$, as well as by the fluctuation at the inner boundary represented with ${b}^{* }(s)$. Subtracting Equation (84) from Equation (82), we obtain the shock motion induced by the upstream perturbation alone, which is denoted by ${(\delta {r}_{\mathrm{sh}}/{r}_{\mathrm{sh}})}_{\mathrm{ex}}(t)={{ \mathcal L }}^{-1}[{(\delta {r}_{\mathrm{sh}}/{r}_{\mathrm{sh}})}_{\mathrm{ex}}^{* }]$:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}\right)}_{\mathrm{ex}}^{* }=\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}-{{ \mathcal I }}^{* }[{{\boldsymbol{y}}}_{0}](0,s)\end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray}&&=-\displaystyle \frac{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)R{{\boldsymbol{z}}}^{* }(s)}{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)(s{\boldsymbol{c}}+{\boldsymbol{d}})},\end{eqnarray} \tag{ 86 }$$

$$\begin{eqnarray}&&=: {{ \mathcal J }}^{* }({{\boldsymbol{z}}}^{* },s),\end{eqnarray} \tag{ 87 }$$

$$\begin{eqnarray}&&={{ \mathcal J }}^{* }({z}_{1}^{* },\ldots ,{z}_{n}^{* },s),\end{eqnarray} \tag{ 88 }$$

$$\begin{eqnarray}&&=\sum _{k=1}^{n}{{ \mathcal J }}_{k}^{* }{z}_{k}^{* }.\end{eqnarray} \tag{ 89 }$$

In the last equation we defined a set of functions,

$$\begin{eqnarray}&&{{ \mathcal J }}_{k}^{* }(s):= {{ \mathcal J }}^{* }(0,\ldots ,0,1,0,\ldots ,0,s)\ (k=1,\ldots ,n),\end{eqnarray} \tag{ 90 }$$

where the arguments of ${{ \mathcal J }}^{* }(\ldots )$ on the right-hand side are set to be zero except for the kth one, which is put to unity. Since ${ \mathcal L }[\delta (t)]=1$, ${{ \mathcal J }}_{k}(t)={{ \mathcal L }}^{-1}[{{ \mathcal J }}_{k}^{* }]$ describes $\delta {r}_{\mathrm{sh}}/{r}_{\mathrm{sh}}$ for the impulsive perturbation of the unit strength, i.e., $\delta (t)$, added only to the kth component of ${\boldsymbol{z}}(t)$. This implies that ${{ \mathcal J }}_{k}(t)$ can be regarded as a Green's function. In fact, recalling that the Laplace transform of the convolution $(f\ast g)(t)$ is the product of the Laplace transforms of the individual functions, ${ \mathcal L }[f]{ \mathcal L }[g]={f}^{* }(s){g}^{* }(s)$ (e.g., Schiff 1999), we obtain

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}\right)}_{\mathrm{ex}}=\sum _{k=1}^{n}{{ \mathcal L }}^{-1}[{{ \mathcal J }}_{k}^{* }{z}_{k}^{* }]\end{eqnarray} \tag{ 91 }$$

$$\begin{eqnarray}&&=\sum _{k=1}^{n}({{ \mathcal J }}_{k}\ast {z}_{k})(t)\end{eqnarray} \tag{ 92 }$$

$$\begin{eqnarray}&&=\sum _{k=1}^{n}{\int }_{0}^{t}{{ \mathcal J }}_{k}(t-\tau ){z}_{k}(\tau )d\tau .\end{eqnarray} \tag{ 93 }$$

We can hence obtain the evolution of the shock radius for any perturbation once we know ${{ \mathcal J }}_{k}\ (k=1,\ldots ,n)$.

The function ${{ \mathcal I }}^{* }[0,{{\boldsymbol{y}}}_{0}](s)$ can be further divided into two parts:

$$\begin{eqnarray}&&{{ \mathcal I }}^{* }[{{\boldsymbol{y}}}_{0}](0,s)={{ \mathcal K }}^{* }[{{\boldsymbol{y}}}_{0}](s)+{{ \mathcal S }}^{* }(s),\end{eqnarray} \tag{ 94 }$$

where ${{ \mathcal K }}^{* }$ and ${{ \mathcal S }}^{* }$ describe the contributions from the initial perturbation in the shocked region and from the fluctuation at the inner boundary, respectively, and given as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}\right)}_{\mathrm{ini}}^{* }:= {{ \mathcal K }}^{* }[{{\boldsymbol{y}}}_{0}](s):= -\displaystyle \frac{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\boldsymbol{h}}}}^{* }[{{\boldsymbol{y}}}_{0}](s)}{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)(s{\boldsymbol{c}}+{\boldsymbol{d}})},\end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}\right)}_{\mathrm{IB}}^{* }:= {{ \mathcal S }}^{* }(s):= \displaystyle \frac{{b}^{* }(s)}{{{\boldsymbol{a}}}^{* }(s)\cdot {\tilde{{\rm{\Lambda }}}}^{* }(s)(s{\boldsymbol{c}}+{\boldsymbol{d}})},\end{eqnarray} \tag{ 96 }$$

Then the whole evolution of perturbed shock radius is given as the sum of these contributions:

$$\begin{eqnarray}&&\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}={\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}\right)}_{\mathrm{ini}}+{\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}\right)}_{\mathrm{ex}}+{\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}\right)}_{\mathrm{IB}}.\end{eqnarray} \tag{ 97 }$$

Let us write the evolution of the perturbed shock radius as

$$\begin{eqnarray}&&\displaystyle \frac{\delta {r}_{\mathrm{sh}}}{{r}_{\mathrm{sh}}}=\sum _{j}{a}_{j}{e}^{{{\rm{\Omega }}}_{j}t}\sin ({\omega }_{j}t+{\phi }_{j}),\end{eqnarray} \tag{ 98 }$$

and determine the amplitudes a_j for a given upstream perturbation. In Equation (98) ${{\rm{\Omega }}}_{j}$ and ${\omega }_{j}$ are the growth rate and oscillation frequency of the jth mode ($j=1,2,3,\ldots$), respectively. To ensure the uniqueness of this expansion, we assume that ${\omega }_{j}\geqslant 0$ and $-\pi /2\leqslant {\phi }_{j}\lt \pi /2$ for each j. We also note that ${{\rm{\Omega }}}_{j}$ and a_j are real numbers. The amplitude a_j is related to the residue of $\delta {r}_{\mathrm{sh}}^{* }/{r}_{\mathrm{sh}}$ at the corresponding pole in the complex plane, which is directly calculated as

$$\begin{eqnarray}&&\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}=\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}\sum _{j}{ \mathcal L }[{a}_{j}{e}^{{{\rm{\Omega }}}_{j}t}\sin ({\omega }_{j}t+{\phi }_{j})]\end{eqnarray} \tag{ 99 }$$

$$\begin{eqnarray}&&=\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}\sum _{j}\displaystyle \frac{{a}_{j}}{2i}\left(\displaystyle \frac{{e}^{i{\phi }_{j}}}{s-{{\rm{\Omega }}}_{j}-i{\omega }_{j}}-\displaystyle \frac{{e}^{-i{\phi }_{j}}}{s-{{\rm{\Omega }}}_{j}+i{\omega }_{j}}\right),\end{eqnarray} \tag{ 100 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{{a}_{k}}{2i}{e}^{i{\phi }_{k}}.\end{eqnarray} \tag{ 101 }$$

On the other hand, the residue can be numerically obtained by means of the Cauchy theorem:

$$\begin{eqnarray}&&\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}=\displaystyle \frac{1}{2\pi i}{\oint }_{C}\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}{dz},\end{eqnarray} \tag{ 102 }$$

where C is any closed curve in the complex plane that includes only the kth pole inside. The integral is conducted numerically by solving the initial boundary value problem given by Equations (13)–(15) for a set of s on the contour C.

Once the residue for the kth pole is obtained from Equation (102), the amplitude a_k is derived from Equation (101) as follows:

$$\begin{eqnarray}&&{a}_{k}=-2({v}_{k}\cos {\phi }_{k}-{u}_{k}\sin {\phi }_{k})+2i({u}_{k}\cos {\phi }_{k}+{v}_{k}\sin {\phi }_{k}),\end{eqnarray} \tag{ 103 }$$

where the u_k and v_k are the real and imaginary parts of the residue, respectively. Since a_k is real, the imaginary part must vanish. Therefore, we find the result

$$\begin{eqnarray}{a}_{k} & = & -\,2\mathrm{sgn}({v}_{k})\sqrt{{u}_{k}^{2}+{v}_{k}^{2}}\\ & = & -\,2\mathrm{sgn}\left[\mathrm{Im}\left(\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}\right)\right]\left|\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}\displaystyle \frac{\delta {r}_{\mathrm{sh}}^{* }}{{r}_{\mathrm{sh}}}\right|,\end{eqnarray} \tag{ 104 }$$

and

$$\begin{eqnarray}&&\tan {\phi }_{k}=-\displaystyle \frac{{u}_{k}}{{v}_{k}}.\end{eqnarray} \tag{ 105 }$$

Corresponding to the decomposition of $\delta {r}_{\mathrm{sh}}/{r}_{\mathrm{sh}}$, Equation (97), its residue can also be divided into three parts:

$$\begin{eqnarray}| {a}_{k}| & = & 2\left|\displaystyle \sum _{j=1}^{n}\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}[{z}_{j}^{* }(s){{ \mathcal J }}_{j}^{* }(s)]\\ & & +\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}[{{ \mathcal K }}^{* }[{{\boldsymbol{y}}}_{0}](s)+{{ \mathcal S }}^{* }(s)]| \phantom{\rule{}{3.75ex}}.\end{eqnarray} \tag{ 106 }$$

If we assume that ${{\boldsymbol{z}}}^{* }(s)$ does not have a pole at $s={{\rm{\Omega }}}_{k}+i{\omega }_{k}$, which is normally the case, the first term on the right-hand side of the above equation is further calculated as

$$\begin{eqnarray}\displaystyle \sum _{j=1}^{n}\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}[{z}_{j}^{* }(s){{ \mathcal J }}_{j}^{* }(s)] & = & \displaystyle \sum _{j=1}^{n}{z}_{j}^{* }({{\rm{\Omega }}}_{k}+i{\omega }_{k})\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}{{ \mathcal J }}_{j}^{* }(s)\\ & & =: \,\displaystyle \sum _{j=1}^{n}{W}_{j}^{* }({s}_{k}).\end{eqnarray} \tag{ 107 }$$

Hence, the contribution from the upstream perturbations can be easily obtained for any z(t) once we derive the residue of ${{ \mathcal J }}_{j}^{* }(s)$ $(j=1,\ldots ,n)$. If there is no perturbation initially in the shocked region and the inner boundary does not produce fluctuations, the amplitude is solely determined by the upstream perturbation as

$$\begin{eqnarray}&&| {a}_{k}| =2\left|\sum _{j=1}^{n}{W}_{j}^{* }({s}_{k})\right|\end{eqnarray} \tag{ 108 }$$

$$\begin{eqnarray}&&=2\sqrt{\sum _{j=1}^{n}{| {W}_{j}^{* }({s}_{k})| }^{2}+\displaystyle \sum _{i\ne j}\mathrm{Re}[{W}_{i}^{* }({s}_{k})\overline{{W}_{j}^{* }({s}_{k})}]},\end{eqnarray} \tag{ 109 }$$

where the horizontal line over functions means their complex conjugates.

We here derive the excited amplitude for the generic upstream perturbation expressed in the following form:

$$\begin{eqnarray}&&{z}_{j}(t)=\sum _{\sigma }{b}_{j\sigma }\sin ({\omega }_{j\sigma }t+{\phi }_{j\sigma }).\end{eqnarray} \tag{ 110 }$$

The corresponding Laplace-transformed perturbation is given by

$$\begin{eqnarray}&&{z}_{j}^{* }(t)=\sum _{\sigma }{b}_{j\sigma }\displaystyle \frac{s\sin {\phi }_{j\sigma }+{\omega }_{j\sigma }\cos {\phi }_{j\sigma }}{{s}^{2}+{\omega }_{j\sigma }^{2}}.\end{eqnarray} \tag{ 111 }$$

Inserting the above expression into Equations (108) and (109), we obtain

$$\begin{eqnarray}&&| {a}_{k}| =2\left|\sum _{\sigma }\sum _{j=1}^{n}{W}_{j\sigma }^{* }({s}_{k})\right|\end{eqnarray} \tag{ 112 }$$

$$\begin{eqnarray}&&=2\sqrt{\sum _{\sigma }{\left|\sum _{j=1}^{n}{W}_{j\sigma }^{* }({s}_{k})\right|}^{2}+\displaystyle \sum _{\sigma \ne \rho }\mathrm{Re}\left[\sum _{i,j=1}^{n}{W}_{i\sigma }^{* }({s}_{k}){W}_{j\rho }^{* }({s}_{k})\right]},\end{eqnarray} \tag{ 113 }$$

where ${W}_{j\sigma }^{* }$ is given by

$$\begin{eqnarray}&&{W}_{j\sigma }^{* }({s}_{k})\\ &&\quad := {b}_{j\sigma }\displaystyle \frac{({{\rm{\Omega }}}_{k}+i{\omega }_{k})\sin {\phi }_{j\sigma }+{\omega }_{j\sigma }\cos {\phi }_{j\sigma }}{{({{\rm{\Omega }}}_{k}+i{\omega }_{k})}^{2}+{\omega }_{j\sigma }^{2}}\mathop{\mathrm{Res}}\limits_{s={{\rm{\Omega }}}_{k}+i{\omega }_{k}}{{ \mathcal J }}_{j}^{* }(s).\end{eqnarray} \tag{ 114 }$$

We define ${a}_{k\sigma }$ as the amplitude of the kth mode that would be excited by a single harmonic perturbation with a frequency ${\omega }_{j\sigma }$. Recalling a relation $| {a}_{k\sigma }| =2| {\sum }_{j}{W}_{j\sigma }({s}_{k})|$, we rewrite Equation (113) as

$$\begin{eqnarray}&&| {a}_{k}| =\sqrt{\sum _{\sigma }| {a}_{k\sigma }{| }^{2}+4\displaystyle \sum _{\sigma \ne \rho }\mathrm{Re}\left[\sum _{i,j=1}^{n}{W}_{i\sigma }^{* }({s}_{k}){W}_{j\rho }^{* }({s}_{k})\right]},\end{eqnarray} \tag{ 115 }$$

$$\begin{eqnarray}&&\sim \sqrt{\sum _{\sigma }| {a}_{k\sigma }{| }^{2}},\end{eqnarray} \tag{ 116 }$$

where the last approximation holds unless ${\sum }_{j}{W}_{j\sigma }({s}_{k})$ $(\sigma =1,2,\ldots )$ are correlated with each other. The excited amplitude is then expressed as a square root of the sum of individual amplitudes squared and will be approximately given by the largest ${a}_{k\sigma }$.