© The Authors 2024

1. Introduction

The explosive death of massive stars is sensitive to the development of hydrodynamical instabilities, which break the spherical symmetry of the stellar core, affect the efficiency of neutrino energy absorption, and generate turbulence (Müller 2020; Janka et al. 2016; Burrows & Vartanyan 2021). A spherically symmetric scenario based on radial motions seems possible only for the lightest progenitors (Kitaura et al. 2006; Stockinger et al. 2020). Asymmetric motions contribute to the kick and spin of the neutron star (Müller et al. 2019; Janka 2017), the emission of gravitational waves (Kotake & Kuroda 2017), and a modulation of neutrino emission (Müller 2019b; Tamborra & Murase 2019). Drago et al. (2023) considered the correlation of instability signatures in the gravitational waves and neutrino signals in order to improve the efficiency with which we can detect gravitational waves from nearby supernovae. Ultimately, the information encoded in the gravitational waves and neutrino signals can be used to recover the properties of the dying star and its explosion mechanism (Powell & Müller 2022). The computational cost of 3D numerical simulations precludes systematic coverage of the large parameter space describing the initial conditions of stellar core collapse. Understanding the underlying mechanism of hydrodynamical instabilities is essential to extrapolating the results of sparse numerical simulations, evaluating the impact of additional physical ingredients, and designing effective prescriptions for parametric studies (Müller 2019a).

Among the hydrodynamical instabilities at work during the phase of stalled accretion shock, the standing accretion shock instability (SASI; Blondin et al. 2003) is able to introduce coherent transverse motions with a large angular scale, growing over a timescale related to the advection time from the shock to the neutron star surface. The mechanism of SASI has been described as an advective-acoustic cycle between the shock and the vicinity of the proto-neutron star (Foglizzo et al. 2007; Fernández & Thompson 2009b; Scheck et al. 2008). Our analytical understanding of SASI eigenfrequencies in the WKB approximation is however restricted to the limit of a large shock radius for high-frequency overtones (Foglizzo et al. 2007), while the lowest-frequency fundamental mode is often the most unstable. A fully analytic solution including the fundamental mode was obtained only in a very idealised model, where the stationary flow is plane parallel and uniform except in a compact region of deceleration that mimics the vicinity of the neutron star (Foglizzo 2009) (hereafter F09). This toy model neglected the flow gradients that are extended all the way from the shock to the neutron star and the non-adiabatic character of the neutrino processes taking place in the vicinity of the neutron star. Despite these limitations, this model illustrated the interplay of the advective-acoustic and the purely acoustic cycles, and the expected phase mixing taking place at high frequency. The model was used by Guilet & Foglizzo (2012) to interpret the frequency spacing of the eigenspectrum in the framework of the advective-acoustic mechanism rather than a purely acoustic process. The physical understanding of SASI led Müller & Janka (2014) to define an empirical formula for the SASI oscillation period in order to interpret the modulation of the neutrino signal. Directly proportional to an approximation of the advection time from the shock to the neutron star surface, it was tested on a 2D numerical simulation of the collapse of a 25 M_⊙ progenitor. However, whether or not it can be generalised for other progenitors, as proposed by Müller (2019b), remains unclear.

The aim of the present study is to improve our understanding of the fundamental mode of SASI in spherical geometry in order to better identify the physical parameters governing its oscillation frequency. Simple cooling functions mimicking neutrino emission are used to assess the role of the cooling region. An adiabatic approximation is also used to assess the role of the advective-acoustic coupling in the radially extended region between the shock and the proto-neutron star. The adiabatic approximation is motivated by the adiabatic simulations of Blondin & Mezzacappa (2007) and by the shallow-water experiment (Foglizzo et al. 2012, 2015), which both suggest that several properties of SASI may be understood using adiabatic equations. Dunham et al. (2024) have also used the adiabatic approximation to analyse the impact of general relativistic corrections on SASI eigenfrequencies.

The set of differential equations defining the eigenfrequencies of spherical accretion are recalled in Sect. 2, with particular attention being paid to the radial extension of the non-adiabatic cooling layer and its impact on the oscillation period of SASI. The eigenfrequencies calculated in the adiabatic approximation are compared to those including neutrino losses in Sect. 3. The adiabatic model is formulated as a self-forced oscillator in Sect. 4, with eigenfrequencies defined by an integral equation. An analytic approximation of this equation is outlined and analysed in Sect. 5. Conclusions and perspectives are formulated in Sect. 6.

2. Stationary accretion with non-adiabatic cooling

2.1. Stationary flow

We use the same general framework as Blondin et al. (2003), Foglizzo et al. (2007), Blondin et al. (2017) to describe the phase where the accretion shock stalls at the radius r_sh. Neutrino absorption is neglected and neutrino emission near the proto-neutron star radius r_ns is idealised with a cooling function of the density ρ and pressure p:

$$\begin{array}{c}\hfill \mathcal{L}=-A_{\mathrm{cool}}{\rho }^{\beta -\alpha }{p}^{\alpha }.\end{array}$$(1)

The collapsing stellar core immediately after bounce is modelled as a perfect gas with an adiabatic index of γ = 4/3, dominated by the degeneracy pressure of relativistic electrons. The gravitational potential Φ ≡ −GM_ns/r is assumed to be dominated by the mass M_ns of the proto-neutron star in the Newtonian approximation for simplicity. We neglect the self-gravity of the accreting gas and the increase in M_ns with time. The dimensionless measure of the entropy S is defined as

$$\begin{array}{c}\hfill S\equiv \frac{1}{\gamma -1}log\frac{p}{{\rho }^{\gamma }}.\end{array}$$(2)

The stationary flow is described by the mass conservation, the entropy equation, and the Euler equation:

$$\begin{array}{cc}\hfill \rho v& ={\left(\frac{r_{\mathrm{sh}}}{r}\right)}^g{\rho }_{\mathrm{sh}}{v}_{\mathrm{sh}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\partial S}{\partial r}& =\frac{\mathcal{L}}{\mathit{pv}},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \frac{\partial }{\partial r}(\frac{v^2}{2}+\frac{c^2}{\gamma -1}+\mathrm{\Phi })& =\frac{\mathcal{L}}{\rho v},\hfill \end{array}$$(5)

where c ≡ (γp/ρ)^1/2 is the sound speed and v < 0 the radial velocity. The geometrical parameter g = 2 accounting for the spherical geometry is noted as a parameter in this section and Sect. 3 to keep track of its impact on the equations in the spirit of Walk et al. (2023). The subscript ‘sh’ in Eq. (3) refers to quantities immediately below the shock, and the subscript ‘1’ refers to pre-shock quantities. The jump conditions at the shock follow from the conservation of mass flux and momentum flux, taking into account the energy lost across the shock through nuclear dissociation. This latter is modelled as in Fernández & Thompson (2009b), Fernández et al. (2014) by the parameter ε, which is a measure of the energy loss Δe_disso per unit of mass in units of the pre-shock kinetic energy density v₁²/2 as in Huete et al. (2018):

$$\begin{array}{ccc}\hfill \epsilon & \equiv \hfill & \hfill \frac{\mathrm{\Delta }{e}_{\mathrm{disso}}}{v_1^2/2}.\end{array}$$(6)

The effect of ε on the post-shock Mach number ℳ_sh and the Rankine-Hugoniot conditions is described in Appendix A, following Eqs. (A.4)–(A.6) of Foglizzo et al. (2006). The pre-shock deceleration effect of pressure is neglected, with v₁² ∼ 2GM_ns/r_sh. Defining the Mach number ℳ ≡|v|/c as positive, we assume a strong adiabatic shock ℳ₁ ≫ 1 in the numerical calculations of this section.

The adiabatic compression of the post-shock flow by the gravitational potential Φ produces an inward increase in the enthalpy c²/(γ − 1) and thus an increase in the gas temperature T ∝ c² according to the Bernoulli Eq. (5). In our model of stationary accretion, the temperature profile displays a maximum at a radius denoted r_peak, where the energy losses by neutrino emission balance the adiabatic heating due to the gravitational compression. The width of the cooling layer depends a priori on the parameters (α, β) of the cooling function ℒ, the mass accretion rate, and the geometry. With α = 3/2 and β = 5/2, Walk et al. (2023) noted that both the radius r_peak ∼ 1.2r_ns of maximum temperature and the typical width Δr_peak ∼ 1.5r_ns of the temperature profile measured at half maximum seemed surprisingly independent of both the mass-accretion rate and the shock radius, and also seemed independent of the geometry. Closer inspection confirms that the location of r_peak varies by less than 1% between the cylindrical and spherical geometries, and varies by less than 0.5% with r_sh for r_sh/r_ns > 3.

Figure 1 compares r_peak to the radius of maximum deceleration denoted r_∇, with r_∇ < r_peak over a large range of parameters (α, β), except for when β approaches unity. Noting that r_∇ = r_ns for α < β (Foglizzo et al. 2007), Fig. 1 indicates that r_peak/r_∇ ∼ 1.2 for (α, β) = (3/2, 5/2), and r_peak/r_∇ ∼ 1.0 for (α, β) = (6, 1). In Sect. 2.2, the parameters (α, β) are varied continuously from (3/2, 5/2) to (6, 1) and from (3/2, 5/2) to (5, 6) in order to evaluate their impact on SASI properties.

2.2. Perturbed flow

The perturbations of the stationary flow are characterised by the wavenumbers ℓ, m of spherical harmonics and a complex eigenfrequency ω. The eigenfrequency is independent of the azimuthal wavenumber |m|≤ℓ as established in Foglizzo et al. (2007). We use the same physical variables as in Foglizzo et al. (2006) guided by analytical simplicity, namely the perturbation of the Bernoulli constant δf, the perturbed mass flux δh, and the perturbed dimensionless entropy δS. The perturbation δK is a combination of perturbed entropy δS and the quantity δw_⊥ defined from the radial component of the curl of the vorticity perturbation δw ≡ ∇ × δv:

$$\begin{array}{cc}\hfill \delta f& \equiv v\delta v_r+\frac{\delta c^2}{\gamma -1},\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \delta h& \equiv \frac{\delta v_r}{v}+\frac{\delta \rho }{\rho },\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \delta S& \equiv \frac{1}{\gamma -1}\frac{\delta c^2}{c^2}-\frac{\delta \rho }{\rho },\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill \delta w_{\perp }& \equiv r{(\mathrm{\nabla }\times \delta w)}_r,\hfill \end{array}$$(10)

$$\begin{array}{cc}& =\frac{1}{sin\theta }[\frac{\partial }{\partial \theta }(\delta w_{\phi }sin\theta )-im\delta w_{\theta }],\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill \delta K& \equiv rv\delta w_{\perp }+\ell (\ell +1)\frac{c^2}{\gamma }\delta S.\hfill \end{array}$$(12)

We introduce the quantity δA, which is defined as the divergence of the ortho-radial velocity δv_⊥ ≡ (0, δv_θ, δv_φ):

$$\begin{array}{cc}\hfill \delta A& \equiv r^2\mathrm{\nabla }·\delta v_{\perp },\hfill \end{array}$$(13)

$$\begin{array}{cc}& =\frac{r}{sin\theta }[\frac{\partial }{\partial \theta }(\delta v_{\theta }sin\theta )+im\delta v_{\phi }].\hfill \end{array}$$(14)

We establish in Appendix B that δA and δw_⊥ are related to the perturbation δv_r of the radial velocity as

$$\begin{array}{c}\hfill \delta v_r=\frac{r\delta w_{\perp }}{\ell (\ell +1)}-\frac{1}{\ell (\ell +1)}\frac{\partial \delta A}{\partial r}.\end{array}$$(15)

This equation invites us to interpret δA/ℓ(ℓ + 1) as the potential defining the compressible part of the perturbed velocity, and rδw_⊥/ℓ(ℓ + 1) as the rotational contribution to the radial velocity perturbation. Using the transverse components of the Euler equation Eqs. (B.10), (B.11) in Appendix B with Eq. (12), we note that δA is related to δK and δf as

$$\begin{array}{c}\hfill \delta K=i\omega \delta A+\ell (\ell +1)\delta f.\end{array}$$(16)

The differential system satisfied by (δA, δh, δS, δK) is as follows:

$$\begin{array}{cc}\hfill (\frac{1-{\mathcal{M}}^2}{v}\frac{\partial }{\partial r}+\frac{i\omega }{c^2})\frac{\delta A}{\ell (\ell +1)}& =-\delta h-(\gamma -1+\frac{1}{{\mathcal{M}}^2})\frac{\delta S}{\gamma }\hfill \\ \hfill & +\frac{\delta K}{v^2\ell (\ell +1)},\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill (\frac{1-{\mathcal{M}}^2}{v}\frac{\partial }{\partial r}+\frac{i\omega }{c^2})\delta h& =\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}\frac{\delta A}{\ell (\ell +1)}\hfill \\ \hfill & -\frac{i\omega }{v^2}\delta S+\frac{i\omega }{v^2{c}^2}\frac{\delta K}{\ell (\ell +1)},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill (\frac{\partial }{\partial r}-\frac{i\omega }{v})\delta S& =\delta \left(\frac{\mathcal{L}}{\mathit{pv}}\right),\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill (\frac{\partial }{\partial r}-\frac{i\omega }{v})\frac{\delta K}{\ell (\ell +1)}& =\delta \left(\frac{\mathcal{L}}{\rho v}\right).\hfill \end{array}$$(20)

This system includes explicit non-adiabatic terms only in the Eqs. (19, 20) governing δS and δK. In Eq. (18), the Lamb frequency ω_Lamb associated with the spherical harmonic of order ℓ defines the turning point of non-radial acoustic waves:

$$\begin{array}{ccc}\hfill {\omega }_{\mathrm{Lamb}}^2& \equiv \hfill & \hfill \ell (\ell +1)\frac{c^2}{r^2}(1-{\mathcal{M}}^2).\end{array}$$(21)

The boundary conditions at the shock are reformulated from Eqs. (28, 29, E.7, E.8) in Foglizzo et al. (2006) using Eq. (16):

$$\begin{array}{cc}\hfill \frac{\delta A_{\mathrm{sh}}}{\ell (\ell +1)}& =-(1-\frac{v_{\mathrm{sh}}}{v_1})v_1\mathrm{\Delta }\zeta ,\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill \delta h_{\mathrm{sh}}& =(1-\frac{v_{\mathrm{sh}}}{v_1})\frac{\mathrm{\Delta }v}{v_{\mathrm{sh}}},\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill \delta S_{\mathrm{sh}}& =\gamma \frac{v_1}{c_{\mathrm{sh}}^2}(i\omega +{\omega }_{\mathrm{\Phi }})\mathrm{\Delta }\zeta {(1-\frac{v_{\mathrm{sh}}}{v_1})}^2,\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill \delta K_{\mathrm{sh}}& =-\ell (\ell +1)\mathrm{\Delta }\zeta \frac{c_{\mathrm{sh}}^2}{\gamma }{\left[\mathrm{\nabla }S\right]}_1^{\mathrm{sh}},\hfill \end{array}$$(25)

where Δζ and Δv ≡ −iωΔζ are the shock displacement and velocity, respectively. The reference frequency ω_Φ in Eq. (24) is defined as

$$\begin{array}{c}\hfill \frac{{\omega }_{\mathrm{\Phi }}{r}_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}\equiv \frac{1}{2}\frac{v_1}{v_{\mathrm{sh}}}\frac{\frac{2r_{\mathrm{sh}}}{v_1^2}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}-4\frac{v_{\mathrm{sh}}}{v_1}}{1-\frac{v_{\mathrm{sh}}}{v_1}}+\frac{\frac{v_{\mathrm{sh}}}{v_1}}{\gamma {\mathcal{M}}_{\mathrm{sh}}^2}\frac{r_{\mathrm{sh}}{\left[\mathrm{\nabla }S\right]}_1^{\mathrm{sh}}}{{(1-\frac{v_{\mathrm{sh}}}{v_1})}^2},\end{array}$$(26)

where ω_Φ defines the threshold frequency separating the regime (ω ≪ ω_Φ) – where the phase of δS_sh is opposed to that of Δζ – from the regime (ω ≫ ω_Φ), where the phase of δS_sh is set by that of Δv.

The first two differential equations (17)–(18) are transformed into the following second-order differential equation:

$$\begin{array}{c}\hfill [{(\frac{1-{\mathcal{M}}^2}{v}\frac{\partial }{\partial r}+\frac{i\omega }{c^2})}^2+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}]\delta A=\\ \hfill \frac{\partial }{\partial r}\frac{r\delta w_{\perp }}{v}-\ell (\ell +1)\frac{\gamma -1}{\gamma }\delta \left(\frac{\mathcal{L}}{\mathit{pv}}\right).\end{array}$$(27)

We recognise an acoustic oscillator modified by advection on the left-hand side of Eq. (27), with a forcing on the right-hand side driven by vorticity perturbations and by the perturbation of neutrino cooling. Interestingly, part of this forcing can be studied analytically in the adiabatic approximation defined in Sect. 3.

The variables δv_r, δρ, δc, and δp are deduced from δA, δh, δS, and δK using Eqs. (B.4)–(B.7), as detailed in Appendix B. In particular, using Eq. (B.4) to express the lower boundary condition δv_r(r_ns) = 0, leads to

$$\begin{array}{c}\hfill {(c^2\delta h+c^2\delta S-\delta f)}_{\mathrm{ns}}=0.\end{array}$$(28)

With a vanishing sound speed c_ns = 0 resulting from the hypothesis of stationary flow, the lower boundary condition δv_r(r_ns) = 0 is equivalent to δf_ns = 0.

The boundary value problem is solved numerically using a shooting method from the shock to the inner boundary, iterating over the value of the complex eigenfrequency ω. The mathematical singularity at r_ns is overcome numerically by using logℳ as an integration variable, as in Foglizzo et al. (2007). The inner boundary is then defined as the point where the Mach number has reached a sufficiently small value (ℳ ∼ 10⁻⁹).

Foglizzo et al. (2007) noted similar properties of the SASI harmonics as functions of r_sh/r_∇ with (α, β) = (3/2, 5/2) and (6, 1) in their analysis. We extend this analysis for the fundamental mode and note in Fig. 2 that the normalised growth rate and oscillation frequency are asymptotically independent of (α, β) for r_sh/r_∇ ≫ 1, that is, when the cooling layer is thin compared to the shock radius. This surprising property is checked by varying (α, β) continuously from (3/2, 5/2) to (6, 1) and from (3/2, 5/2) to (5, 6) in Fig. 3. It is remarkable that the normalised oscillation frequency appears to be independent of the cooling process down to a very low ratio r_sh/r_∇: the following analytic formula ω_r^fit provides an approximation for ω_r within 3% for 1.8 < r_sh/r_∇ < 10:

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{adv}}^{\mathrm{adiab}}& \equiv \frac{r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}log\frac{r_{\mathrm{sh}}}{r_{\mathrm{\nabla }}},\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill R_1& \equiv log(\frac{r_{\mathrm{sh}}}{r_{\mathrm{\nabla }}}-1),\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill {\omega }_r^{\mathrm{fit}}& \equiv \frac{2\pi }{{\tau }_{\mathrm{adv}}^{\mathrm{adiab}}}(0.56773+0.28628R_1-0.031763R_1^2).\hfill \end{array}$$(31)

The adiabatic advection time ${\tau }_{\mathrm{adv}}^{\mathrm{adiab}}$ corresponds to a velocity profile increasing linearly with radius as expected asymptotically in spherical geometry for γ = 4/3 (Eq. (19) in Walk et al. 2023).

The important role of r_∇ rather than r_ns was also noted by Scheck et al. (2008), where the strength of SASI in numerical simulations seemed correlated with the abruptness of the deceleration close to the proto-neutron star, and was confirmed in the toy model used by F09 and Guilet & Foglizzo (2012).

Figures 2 and 3 also show that the detailed cooling process affects the growth rate of SASI (solid lines) more than its oscillation frequency (dashed lines) when r_sh < 4r_∇.

We note from Fig. 4 that the oscillation period T_SASI ≡ 2π/ω_r of the fundamental ℓ = 1 SASI mode can be significantly longer than the approximate advection time ${\tau }_{\mathrm{adv}}^{\mathrm{adiab}}$ for a small ratio r_sh/r_∇. The relationship between these two timescales is further discussed in reference to a simplified model in Sect. 5.3. The modest impact of the specificities of the parametrised cooling process on the fundamental SASI eigenfrequency encouraged us to apply our simple model to a more realistic physical model of core collapse in Sect. 2.3. We were also inspired by this finding to look for a deeper analytic understanding of the SASI mechanism using an adiabatic approximation, as detailed in Sect. 3.

2.3. Comparison of the perturbative model with the empirical formula proposed by Müller & Janka (2014)

The adiabatic approximation of the advection time was also used in the empirical formula (33), which describes the oscillation period T_SASI in the core-collapse simulation of a 25 M_⊙ progenitor (Müller & Janka 2014), which is inspired by the physics of the advective-acoustic cycle. r_∇ was approximated with r_ns and the time variation of the mass of the proto-neutron star was neglected:

$$\begin{array}{c}\hfill T_{\mathrm{SASI}}^{\mathrm{MJ}}\equiv 19\mathrm{ms}{\left(\frac{r_{\mathrm{sh}}}{100\mathrm{km}}\right)}^{\frac{3}{2}}log\left(\frac{r_{\mathrm{sh}}}{r_{\mathrm{ns}}}\right).\end{array}$$(32)

The SASI modulation of the neutrino signal was identified in their model s25 between t = 0.12 s and t = 0.45 s post-bounce. The shock radius r_sh decreases from 125 km to 55 km and the ratio r_sh/r_ns varies between 1.8 and 2.3 according to Figs. 3 and 6 of Müller & Janka (2014). The formula (32) well captures the $r_{\mathrm{sh}}^{3/2}$-dependence, which is the main source of variability of T_SASI during the collapse. The overall accuracy of this formula is ∼26% for model s25 as shown in Fig. 5.

The expected variation of T_SASI estimated from Eq. (31) is as follows:

$$\begin{array}{c}\hfill T_{\mathrm{SASI}}=\frac{T_{\mathrm{SASI}}^{\mathrm{MJ}}\times 1.82{\left(\frac{1.7M_{\odot }}{M_{\mathrm{ns}}}\right)}^{\frac{1}{2}}}{0.56773+0.28628R_1-0.031763R_1^2}.\end{array}$$(33)

According to Fig. 5, the overall accuracy of formula (33) applied to the model s25 is improved to ∼10%, which is remarkable given the simplicity of the perturbative model, which does not involve any adjusted parameter and neglects dissociation at the shock. The mass increase M_ns/M_⊙ ∼ 1.7 − 2 estimated from Fig. 2 in Müller & Janka (2014) contributes to a 8% decrease in T_SASI. A more physical estimate should also take into account dissociation, which can significantly increase T_SASI when r_sh is large, and the distinction between r_ns and r_∇, which can significantly decrease T_SASI when r_sh/r_ns is smallest.

Following Fig. 6, the impact of dissociation on the oscillation frequency of SASI is tentatively approximated as

$$\begin{array}{c}\hfill {\omega }_r(\epsilon )={\omega }_r(0)-{\left(\frac{2G{M}_{\mathrm{ns}}}{r_{\mathrm{sh}}^3}\right)}^{\frac{1}{2}}\epsilon (0.19817+0.74804\epsilon ).\end{array}$$(34)

We note that this analytical formula is tested only in the narrow range 1.8 < r_sh/r_∇ < 2.3 and is only meant to provide as an order of magnitude estimate. The dissociation parameter is approximated according to Fig. 12 in Huete et al. (2018):

$$\begin{array}{c}\hfill \epsilon \sim 0.5\left(\frac{r_{\mathrm{sh}}}{150\mathrm{km}}\right)\left(\frac{1.3M_{\odot }}{M_{\mathrm{ns}}}\right),\end{array}$$(35)

with a saturation at ε ∼ 0.5 due to partial dissociation of α-nuclei (Fernández & Thompson 2009a). Thus, Eq. (34) is rewritten as

$$\begin{array}{c}\hfill \frac{10\mathrm{ms}}{T_{\mathrm{SASI}}}=\frac{10\mathrm{ms}}{T_{\mathrm{SASI}}^0}-\frac{1.069\epsilon (0.19817+0.74804\epsilon )}{{\left(\frac{r_{\mathrm{sh}}}{100\mathrm{km}}\right)}^{\frac{3}{2}}{\left(\frac{1.7M_{\odot }}{M_{\mathrm{ns}}}\right)}^{\frac{1}{2}}},\end{array}$$(36)

with

$$\begin{array}{c}\hfill \frac{T_{\mathrm{SASI}}^0}{10.42\mathrm{ms}}\equiv \frac{{\left(\frac{r_{\mathrm{sh}}}{100\mathrm{km}}\right)}^{\frac{3}{2}}{\left(\frac{1.7M_{\odot }}{M_{\mathrm{ns}}}\right)}^{\frac{1}{2}}log\left(\frac{r_{\mathrm{sh}}}{r_{\mathrm{\nabla }}}\right)}{0.56773+0.28628R_1-0.031763R_1^2}.\end{array}$$(37)

According to Fig. 5, taking into account dissociation improves the accuracy of the empirical formula, but only for a limited range of prescribed ratio 1.25 < r_∇/r_ns < 1.35. Improvement of the accuracy of the present perturbative model beyond the ±10% level does not appear straightforward given its many approximations regarding the microphysics, neutrino interactions, and gravity. Its main improvement compared to the empirical formula (32) is the fact that it contains physically consistent estimates of the impact of the mass, the dissociation at the shock, and the SASI phase inferred from the perturbative analysis. Each of them can potentially affect the oscillation period by several tens of percent according to Eq. (33) and Figs. 4 and 6. By neglecting these effects, the empirical formula in Eq. (32) is not expected to maintain its ∼30% accuracy for progenitors other than s25. The accuracy of Eq. (36) should be tested on other numerical simulations of core collapse, bearing in mind that our prescription for dissociation is very crude.

3. Adiabatic model

Here, the adiabatic character of the flow refers to the region between the shock and the inner boundary, while non-adiabatic processes associated to nuclear dissociation and neutrino emission are incorporated into the boundary conditions. The production of entropy perturbations by the shock displacement Δζ is taken into account. Nuclear dissociation across the shock is formally taken into account using the parameter ε (Eq. (6)), but it is set to zero in the illustrative figures of this section. The goal of the analysis presented in this section is to gain some analytical understanding of the effect of the physical parameters involved in the SASI mechanism.

3.1. Explicit expressions in the adiabatic approximation

With ℒ = 0, the differential Equations (19)–(20) describe the advection of perturbations produced by the shock. Using the boundary conditions (24), (25) with ∇S = 0:

$$\begin{array}{cc}\hfill \delta S& =\delta S_{\mathrm{sh}}{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}r}{v}},\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill \delta K& =0.\hfill \end{array}$$(39)

Here, δS_sh is defined by Eq. (24) with a simpler expression for ω_Φ:

$$\begin{array}{c}\hfill \frac{{\omega }_{\mathrm{\Phi }}{r}_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}\equiv \frac{\frac{1}{2}\frac{v_1}{v_{\mathrm{sh}}}-g}{1-\frac{v_{\mathrm{sh}}}{v_1}}.\end{array}$$(40)

Using Eq. (12) with δK = 0 and Eq. (38) gives the explicit expression for δw_⊥ of

$$\begin{array}{c}\hfill \delta w_{\perp }=-\ell (\ell +1)\frac{c^2}{\gamma rv}\delta S.\end{array}$$(41)

The vorticity associated to the entropy perturbation is also explicitly calculated in Appendix C:

$$\begin{array}{cc}\hfill \delta w_r& =0,\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill \delta w_{\theta }& =-\frac{im{c}^2}{rvsin\theta }\frac{\delta S}{\gamma },\hfill \end{array}$$(43)

$$\begin{array}{cc}\hfill \delta w_{\phi }& =\frac{c^2}{\mathit{rv}}\frac{\partial }{\partial \theta }\frac{\delta S}{\gamma }.\hfill \end{array}$$(44)

Together with Eq. (12), Eqs. (43), (44) demonstrate that δK = 0 can be interpreted as a consequence of the baroclinic production of vorticity.

Equation (16) implies that the wave action δf/ω is directly related to δA:

$$\begin{array}{c}\hfill \frac{\delta f}{i\omega }=-\frac{\delta A}{\ell (\ell +1)}.\end{array}$$(45)

The transverse velocity components (δv_θ, δv_φ) are deduced from δA using the transverse Euler equations Eqs. (B.10), (B.11) with Eqs. (43), (44). Here, δv_r is deduced from Eqs. (15) and (41):

$$\begin{array}{cc}\hfill \delta v_r& =-\frac{c^2}{\gamma v}\delta S-\frac{1}{\ell (\ell +1)}\frac{\partial \delta A}{\partial r},\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill r\delta v_{\theta }& =-\frac{1}{\ell (\ell +1)}\frac{\partial \delta A}{\partial \theta },\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill r\delta v_{\phi }& =-\frac{\mathit{im}}{\ell (\ell +1)}\frac{\delta A}{sin\theta }.\hfill \end{array}$$(48)

The adiabatic solution can be viewed as the asymptotic limit of the model with a cooling function when r_peak ∼ r_ns, corresponding to α ≪ 1 and β ≫ 1 according to Fig. 1. For the sake of simplicity, we explore the adiabatic solution with an inner boundary defined by δv_r = 0, formulated by Eq. (28) with c_ns defined as the adiabatic sound speed at the inner boundary. This prescription could be improved to better account for non-adiabatic processes in the cooling layer. The adiabatic simulations in Blondin & Mezzacappa (2007) also used δv_r = 0 at the inner boundary.

3.2. Quantitative comparison of the eigenfrequencies with and without the region of non-adiabatic cooling

The eigenfrequencies corresponding to the adiabatic model are solved numerically and compared in Fig. 7 to the eigenfrequencies of the non-adiabatic formulation for different values of the shock radius. The overall trends suggest that the main properties of SASI can be qualitatively understood by focusing on adiabatic processes. The fundamental mode is the most unstable one for a small shock radius, and becomes dominated by higher overtones for a larger shock radius. Among the most visible differences with the non-adiabatic model, the adiabatic simplification underestimates the frequency by a factor of ≤1.3 for a large shock radius of r_sh ∼ 10r_ns. The growth rate of the most unstable mode is overestimated by a factor of ≤1.4 for a small shock radius and is underestimated by a factor of ≤1.3 for a very large shock radius.

4. Formulation of the SASI mechanism as a self-forced oscillator

4.1. Derivation of the second-order differential equation

For the sake of mathematical simplicity, we use the radial coordinate X defined as in Foglizzo (2001):

$$\begin{array}{c}\hfill \mathrm{d}\mathrm{X}\equiv \frac{v}{1-{\mathcal{M}}^2}\mathrm{d}r.\end{array}$$(49)

The second-order differential Eq. (27) is simplified to

$$\begin{array}{c}\hfill [{(\frac{\partial }{\partial X}+\frac{i\omega }{c^2})}^2+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}]\delta A=\frac{\partial }{\partial X}\frac{r\delta w_{\perp }}{v}.\end{array}$$(50)

Using Eq. (41), the forcing term of Eq. (50) is thus proportional to the entropy perturbation δS_sh produced by the shock:

$$\begin{array}{cc}\hfill [{(\frac{\partial }{\partial X}+\frac{i\omega }{c^2})}^2+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}]\frac{\delta A}{\ell (\ell +1)}& =-{\mathcal{F}}_S\delta S_{\mathrm{sh}},\hfill \end{array}$$(51)

$$\begin{array}{cc}\hfill {\mathcal{F}}_S& \equiv \frac{\partial }{\partial X}\frac{{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}r}{v}}}{\gamma {\mathcal{M}}^2}.\hfill \end{array}$$(52)

We show in Appendix C that the forced oscillator equation can be rewritten as follows, using boldface for vector quantities:

$$\begin{array}{cc}\hfill [{(\frac{\partial }{\partial X}+\frac{i\omega }{c^2})}^2+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}]\delta \mathit{L}& =\frac{\partial }{\partial X}\frac{r\delta \mathit{w}}{v},\hfill \end{array}$$(53)

$$\begin{array}{cc}& =\frac{\partial }{\partial X}\frac{rv\delta \mathit{w}}{c^2{\mathcal{M}}^2},\hfill \end{array}$$(54)

where δL ≡ r × δv is the perturbed specific angular momentum vector and δw ≡ ∇ × δv is the perturbed vorticity vector. We note from Eqs. (43), (44) that rv|δw|/c² is conserved when advected. The relation between δS_sh and δA_sh is deduced from Eqs. (22) and (24):

$$\begin{array}{c}\hfill \ell (\ell +1)\frac{c_{\mathrm{sh}}^2\delta S_{\mathrm{sh}}}{\gamma \delta A_{\mathrm{sh}}}=-(1-\frac{v_{\mathrm{sh}}}{v_1})(i\omega +{\omega }_{\mathrm{\Phi }}).\end{array}$$(55)

The components of transverse vorticity at the shock are related to the components of perturbed specific angular momentum according to Eqs. (43), (44) and Eqs. (47), (48):

$$\begin{array}{c}\hfill {\left(\frac{rv\delta w_{\theta }}{\delta L_{\theta }}\right)}_{\mathrm{sh}}={\left(\frac{rv\delta w_{\phi }}{\delta L_{\phi }}\right)}_{\mathrm{sh}}=-\ell (\ell +1)\frac{c_{\mathrm{sh}}^2\delta S_{\mathrm{sh}}}{\gamma \delta A_{\mathrm{sh}}}.\end{array}$$(56)

The boundary conditions at the shock and at the inner boundary are written in Appendix C as follows:

$$\begin{array}{cc}\hfill {\frac{\partial \delta A}{\partial X}}_{\mathrm{sh}}& =\frac{i\omega }{v_{\mathrm{sh}}^2}\delta A_{\mathrm{sh}}[1-2\frac{v_{\mathrm{sh}}}{v_1}+(\gamma -1){\mathcal{M}}_{\mathrm{sh}}^2(1-\frac{v_{\mathrm{sh}}}{v_1})]\hfill \\ \hfill & -[1+(\gamma -1){\mathcal{M}}_{\mathrm{sh}}^2]\frac{\delta A_{\mathrm{sh}}}{r_{\mathrm{sh}}{v}_{\mathrm{sh}}}(\frac{v_1}{2v_{\mathrm{sh}}}-2),\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill {\frac{\partial \delta A}{\partial X}}_{\mathrm{ns}}& =\frac{i\omega }{c_{\mathrm{ns}}^2}\delta A_{\mathrm{ns}}-\ell (\ell +1)\frac{1-{\mathcal{M}}_{\mathrm{ns}}^2}{{\mathcal{M}}_{\mathrm{ns}}^2}\frac{\delta S_{\mathrm{sh}}}{\gamma }{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{\mathrm{d}X}{v^2}}.\hfill \end{array}$$(58)

They can also be written with the variables rδv_θ, rδv_φ, and δw_θ, δw_φ, using Eqs. (43), (44) and (47), (48). Equations (51) and (53) are thus equivalent for ℓ ≥ 1. The same advective-acoustic cycle can either be considered as an entropic-acoustic cycle or as a vortical-acoustic cycle.

We define a new perturbative variable δY as follows:

$$\begin{array}{ccc}\hfill \delta Y& \equiv \hfill & \hfill \frac{\delta f}{i\omega }{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}X}{c^2}}=-\frac{\delta A}{\ell (\ell +1)}{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}X}{c^2}}.\end{array}$$(59)

Here, Y₀ is defined as the solution of the homogeneous equation associated to Eq. (51) with δS = 0 and δK = 0, and satisfying the inner boundary condition (58):

$$\begin{array}{cc}\hfill \{\frac{{\partial }^2}{\partial X^2}+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}\}{Y}_0& =0,\hfill \end{array}$$(60)

$$\begin{array}{cc}\hfill {\frac{\partial Y_0}{\partial X}}_{\mathrm{ns}}& =\frac{i\omega }{c_{\mathrm{ns}}^2}{Y}_0(r_{\mathrm{ns}}).\hfill \end{array}$$(61)

We refer to Y₀ as the acoustic structure of the post-shock cavity modified by the radial velocity, in the absence of interaction with advected perturbations. The structure of the acoustic equation and its modification by the advection velocity v are more easily recognised when rewriting Eq. (60) with the variable r:

$$\begin{array}{c}\hfill \{\frac{{\partial }^2}{\partial r^2}+\left(\frac{\partial log}{\partial r}\frac{1-{\mathcal{M}}^2}{v}\right)\frac{\partial }{\partial r}+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{c^2{(1-{\mathcal{M}}^2)}^2}\}{Y}_0=0.\end{array}$$(62)

Identifying the physics of SASI with a forced oscillator enables us to evaluate the efficiency of the coupling depending on two effects: (i) the amplitude of the forcing term (∝ 1/ℳ² in Eq. (52)), and (ii) the phase match between the forcing term (advected wavelength) and the oscillator (acoustic wavelength). The forcing amplitude is strongest where ℳ is smallest, in close proximity to the proto-neutron star, but a strong phase mixing is expected there due to the decrease in the radial wavelength (∝2π|v|/ω_r) of advected perturbations. A trade-off between effects (i) and (ii) favours advected perturbations with a sufficiently low frequency, coupled to the acoustic structure in a region sufficiently far from the proto-neutron star to minimize phase mixing.

This description as a forced oscillator was previously proposed in the context of radial Bondi accretion accelerated into a black hole (Foglizzo 2001, 2002) and in a planar adiabatic toy model of SASI with a localised region of feedback (F09). The present work is the first formulation where the radially extended character of the advective-acoustic coupling is taken into account in a shocked decelerated accretion flow in spherical geometry. We present a quantitative evaluation of the coupling efficiency in the following section using a classical resolution of the forced oscillator with the Wronskian method.

4.2. Integral equation defining the eigenfrequencies

In Appendix D, we derive the integral equation defining the eigenfrequencies, which is equivalent to the full differential system and boundary conditions. It is formulated here with the variable r:

$$\begin{array}{cc}\hfill a{\prime }_1{Y}_0^{\mathrm{sh}}+a{\prime }_2{r}_{\mathrm{sh}}{\left(\frac{\partial Y_0}{\partial r}\right)}_{\mathrm{sh}}& =-{\mathcal{M}}_{\mathrm{sh}}^2{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{i\omega }{v}\frac{\mathrm{d}r}{1-{\mathcal{M}}^2}}{Y}_0^{\mathrm{ns}}\hfill \\ \hfill & -{\int }_{\mathrm{ns}}^{\mathrm{sh}}\frac{\partial }{\partial r}\left(Y_0{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega {\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v}}\right)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{v}\mathrm{d}r}\mathrm{d}r,\hfill \end{array}$$(63)

with a′₁, a′₂ defined as

$$\begin{array}{cc}\hfill a{\prime }_1& \equiv (\gamma -1){\mathcal{M}}_{\mathrm{sh}}^2+\frac{i\omega }{i\omega +{\omega }_{\mathrm{\Phi }}}\frac{\frac{v_{\mathrm{sh}}}{v_1}}{1-\frac{v_{\mathrm{sh}}}{v_1}},\hfill \end{array}$$(64)

$$\begin{array}{cc}\hfill a{\prime }_2& \equiv -\frac{1-{\mathcal{M}}_{\mathrm{sh}}^2}{(1-\frac{v_{\mathrm{sh}}}{v_1})\frac{(i\omega +{\omega }_{\mathrm{\Phi }})r_{\mathrm{sh}}}{v_{\mathrm{sh}}}}.\hfill \end{array}$$(65)

The first two terms associated with Y₀ on the left-hand side of Eq. (63) are acoustic, and are marginally modified by advection. The terms on the right-hand side involve the phase oscillations of the advected perturbations. Despite the integration by part, we recognise in the integral the forcing term of Eq. (52) multiplied by the derivative of the homogeneous solution Y₀. As expected in Sect. 4.1 for the classical problem of a forced harmonic oscillator, this integral characterises the efficiency of the forcing, which depends on both the amplitude profile of the forcing δℱ and the matching of its phase compared to the phase of the oscillator Y₀. An analytic approximation of this integral equation is obtained in Sect. 5 in the asymptotic regime, where the acoustic radial structure is non-oscillatory. We verified numerically that the eigenfrequencies satisfying the single equation Eq. (63) are strictly the same as those obtained in Fig. 7 from the solution of the fourth-order adiabatic differential system (17), (20) with ℒ = 0, with the boundary conditions defined by Eq. (28) and Eqs. (22), (25) with ∇S = 0.

5. Analytical estimate for a large shock radius

5.1. Analytical approximations

5.1.1. Stationary flow

We focus on the case of a strong shock ℳ₁ ≫ 1 and a large shock radius r_sh ≫ r_ns, with the post-shock energy density described by the Bernoulli equation, Eq. (5), dominated by the enthalpy and the gravitational contributions:

$$\begin{array}{ccc}\hfill \frac{c^2}{\gamma -1}& \sim \hfill & \hfill \frac{G{M}_{\mathrm{ns}}}{r}.\end{array}$$(66)

The kinetic energy density associated to the radial velocity is a minor contribution, and is even more negligible when the photodissociation across the shock is taken into account. The parameter ε impacts the post-shock Mach number and velocity according to Eqs. (A.2), (A.3):

$$\begin{array}{cc}\hfill {\mathcal{M}}_{\mathrm{sh}}^2& \sim \frac{\gamma -1}{2\gamma }(1-\epsilon ),\hfill \end{array}$$(67)

$$\begin{array}{cc}\hfill \frac{v_{\mathrm{sh}}}{v_1}& \sim \frac{(\gamma -1)(1-\epsilon )}{2+(\gamma -1)(1-\epsilon )}\le \frac{\gamma -1}{\gamma +1},\hfill \end{array}$$(68)

$$\begin{array}{cc}\hfill \frac{c_{\mathrm{sh}}^2{r}_{\mathrm{sh}}}{G{M}_{\mathrm{ns}}}& \sim \frac{4\gamma (\gamma -1)(1-\epsilon )}{{[2+(\gamma -1)(1-\epsilon )]}^2}.\hfill \end{array}$$(69)

Equation (66) allows a power-law approximation of the radial velocity profile with an exponent α_v deduced from Eq. (A.11). This exponent depends on both the adiabatic index and the geometry g = 2:

$$\begin{array}{cc}\hfill {\alpha }_v& \equiv \frac{1}{\gamma -1}-g,\hfill \end{array}$$(70)

$$\begin{array}{cc}\hfill v& =v_{\mathrm{sh}}{\left(\frac{r}{r_{\mathrm{sh}}}\right)}^{{\alpha }_v},\hfill \end{array}$$(71)

$$\begin{array}{cc}\hfill \frac{{\mathcal{M}}_{\mathrm{sh}}}{\mathcal{M}}& \sim {\left(\frac{r_{\mathrm{sh}}}{r}\right)}^{{\alpha }_v+\frac{1}{2}}.\hfill \end{array}$$(72)

We note that Eq. (66) is valid only for r_ns ≤ r ≪ r_sh, and becomes inaccurate in the vicinity of the shock where the sound speed is particularly sensitive to photodissociation, as seen in Eq. (69).

5.1.2. Acoustic structure of the oscillator

Taking advantage of the low frequency of the fundamental mode for ℓ ≥ 1, we approximate in Appendix E the homogeneous solution Y₀ as a linear combination of power laws independent of the frequency for ω ≪ ω_Lamb:

$$\begin{array}{cc}\hfill Y_0& \sim {\left(\frac{r}{r_{\mathrm{ns}}}\right)}^a[{\left(\frac{r_{\mathrm{ns}}}{r}\right)}^b+\frac{b-a}{b+a}{\left(\frac{r}{r_{\mathrm{ns}}}\right)}^b],\hfill \end{array}$$(73)

$$\begin{array}{cc}\hfill \frac{\partial Y_0}{\partial r}& \sim \frac{b-a}{r_{\mathrm{ns}}}{\left(\frac{r}{r_{\mathrm{ns}}}\right)}^{a-1}[{\left(\frac{r}{r_{\mathrm{ns}}}\right)}^b-{\left(\frac{r_{\mathrm{ns}}}{r}\right)}^b],\hfill \end{array}$$(74)

where the exponents a, b are defined as

$$\begin{array}{cc}\hfill a& \equiv \frac{{\alpha }_v+1}{2},\hfill \end{array}$$(75)

$$\begin{array}{cc}\hfill b& \equiv {[{\left(\frac{{\alpha }_v+1}{2}\right)}^2+\ell (\ell +1)]}^{\frac{1}{2}}.\hfill \end{array}$$(76)

The exponents ±b correspond to the two components of the acoustic structure Y₀, where the component decreasing inwards (+b) is associated to the evanescent profile of pressure perturbations at the shock, and the component decreasing outwards (−b) is associated to the evanescent profile of pressure perturbations at the inner boundary.

This approximate solution of Y₀ satisfies the lower boundary condition only in the asymptotic limit r_ns ≪ r_sh.

5.1.3. Forcing by advected perturbations

We use the following approximation of the advection time τ_adv(r) profile:

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{adv}}(r)& \equiv {\int }_{\mathrm{sh}}^r\frac{\mathrm{d}r}{v},\hfill \end{array}$$(77)

$$\begin{array}{cc}& \sim \frac{1}{{\alpha }_v-1}\frac{r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}[{\left(\frac{r_{\mathrm{sh}}}{r}\right)}^{{\alpha }_v-1}-1]\mathrm{if}{\alpha }_v\ne 1,\hfill \end{array}$$(78)

$$\begin{array}{cc}& \sim \frac{r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}log\frac{r_{\mathrm{sh}}}{r}\mathrm{if}{\alpha }_v=1.\hfill \end{array}$$(79)

In particular, Eq. (79) is applicable for γ = 4/3. In this regime, the advection time τ_adv is asymptotically dominated by the inner region. A numerical calculation indicates that this analytical formula slightly underestimates the adiabatic advection time by less than 6% in the range 1 < r_sh/r_ns < 100.

We examine the accuracy of our approximation of ℳ and Y₀ in Fig. 8.

5.2. Analytical expression of the fundamental eigenfrequency

We approximate Eq. (63) using Eq. (72) and Eq. (73) with ℳ² ≪ 1, except near the shock. We use the approximations (73) and (74) of Y₀ and ∂Y₀/∂r and note that (∂log Y₀/∂log r)_sh ∼ 1/(b + a) with a, b defined by Eqs. (75) and (76). We define N ≡ a₂′+a₁′/(b + a) and introduce the normalised eigenfrequency Z ≡ iωr_sh/|v_sh|:

$$\begin{array}{c}\hfill N\left(Z\right)\equiv \frac{\gamma -1}{b+a}{\mathcal{M}}_{\mathrm{sh}}^2+\frac{1-{\mathcal{M}}_{\mathrm{sh}}^2-\frac{Z}{b+a}\frac{v_{\mathrm{sh}}}{v_1}}{(1-\frac{v_{\mathrm{sh}}}{v_1})(Z+\frac{{\omega }_{\mathrm{\Phi }}{r}_{\mathrm{sh}}}{|v_{\mathrm{sh}}|})}.\end{array}$$(80)

Using the variable x ≡ r/r_ns, for ℓ ≥ 1 the eigenfrequency equation, Eq. (63), is approximated by

$$\begin{array}{c}\hfill -x_{\mathrm{sh}}^{3a-b-1}{\int }_1^{x_{\mathrm{sh}}}(x^b-x^{-b}){\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}r}{v}}\frac{\mathrm{d}x}{x^{3a}}\\ \hfill =N\left(\frac{i\omega r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}\right)+\frac{2b}{\ell (\ell +1)}\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{x_{\mathrm{sh}}^{a+b}}{\mathrm{e}}^{i\omega {\tau }_{\mathrm{adv}}^{\mathrm{ns}}}.\end{array}$$(81)

The approximation of the eigenfrequency with γ = 4/3 allows an explicit calculation of the integral involved in Eq. (81). The following equation defining the eigenfrequencies for ℓ ≥ 1 is expressed using the complex amplification factor 𝒬(ω) per advective-acoustic cycle:

$$\begin{array}{cc}\hfill \mathcal{Q}(Z)& \equiv \frac{2b{\left(\frac{r_{\mathrm{sh}}}{r_{\mathrm{ns}}}\right)}^{2-b}\{1+[{(Z+2)}^2-b^2]\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{\ell (\ell +1)x_{\mathrm{sh}}^3}\}}{[1-(Z+2-b)N](Z+2+b)-\frac{Z+2-b}{x_{\mathrm{sh}}^{2b}}},\hfill \end{array}$$(82)

$$\begin{array}{cc}\hfill \mathcal{Q}\left(\frac{i\omega r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}\right){\mathrm{e}}^{i\omega {\tau }_{\mathrm{adv}}^{\mathrm{ns}}}& =1,\hfill \end{array}$$(83)

where b² = 1 + ℓ(ℓ + 1) and the advection time ${\tau }_{\mathrm{adv}}^{\mathrm{ns}}$ from the shock to the inner boundary is approximated by Eq. (79). The numerical solution of this equation is compared to the exact solution of Eq. (63) in Fig. 9.

The good agreement obtained with the analytic expression (83) demonstrates that the approximation of the flow profile v, ℳ and the approximation of the homogeneous solution Y₀ are sufficient to capture the essence of the instability in the adiabatic model and identify the leading contributions in the integral expression (63). The frequency is overestimated by up to 20% and the growth rate is overestimated by up to 26%. The region of radial propagation of acoustic waves was neglected for the sake of simplicity in the analytical approximation (73) of Y₀, without significant consequences regarding the approximation of the fundamental frequency ω_r, because the radial extension of this region is modest compared to the region affecting the phase of advected perturbations.

The larger inaccuracy of ∼76% of the simple formula ${\omega }_r\sim 2\pi /{\tau }_{\mathrm{adv}}^{\mathrm{ns}}$, also shown in Fig. 9, is mainly due to the neglect of the frequency dependence of the phase of 𝒬, denoted φ_𝒬. The complex Eq. (83) is equivalent to the following set of two real equations:

$$\begin{array}{cc}\hfill {\omega }_i& =\frac{\log |\mathcal{Q}|}{{\tau }_{\mathrm{adv}}^{\mathrm{ns}}},\hfill \end{array}$$(84)

$$\begin{array}{cc}\hfill {\omega }_r& =\frac{2\pi -{\phi }_{\mathcal{Q}}}{{\tau }_{\mathrm{adv}}^{\mathrm{ns}}}.\hfill \end{array}$$(85)

5.3. Oscillation period of the advective-acoustic cycle in the adiabatic approximation

Comparing ω_r and $2\pi /{\tau }_{\mathrm{adv}}^{\mathrm{ns}}$ in Fig. 9, it is particularly striking that the oscillation period T_SASI ≡ 2π/ω_r is significantly longer than the adiabatic advection time ${\tau }_{\mathrm{adv}}^{\mathrm{ns}}$. This contrasts with the simulations in Scheck et al. (2008) where $T_{\mathrm{SASI}}\sim {\tau }_{\mathrm{adv}}^{\mathrm{ns}}$, in line with the analytic toy model of F09. Not only does T_SASI differ from ${\tau }_{\mathrm{adv}}^{\mathrm{ns}}$ by more than 76%, but T_SASI is actually 76% longer than the longest advection timescale ${\tau }_{\mathrm{adv}}^{\mathrm{ns}}$, excluding the possibility that the oscillation period could coincide with any advection time. As explained in Scheck et al. (2008), the relation between T_SASI and ${\tau }_{\mathrm{adv}}^{\mathrm{ns}}$ does depend on the phase φ_𝒬 of the complex efficiency 𝒬.

The adiabatic approximation teaches us that, even in a simple adiabatic model, the phase shift φ_𝒬 should not be neglected a priori. This was also true in the calculation of Sect. 2.3 with non-adiabatic cooling. The value of |𝒬| and φ_𝒬 can be deduced from ω and ${\tau }_{\mathrm{adv}}^{\mathrm{ns}}$ from Eqs. (84), (85):

$$\begin{array}{cc}\hfill |\mathcal{Q}|& ={\mathrm{e}}^{{\omega }_i{\tau }_{\mathrm{adv}}^{\mathrm{ns}}},\hfill \end{array}$$(86)

$$\begin{array}{cc}\hfill {\phi }_{\mathcal{Q}}& =2\pi -{\omega }_r{\tau }_{\mathrm{adv}}^{\mathrm{ns}}.\hfill \end{array}$$(87)

This characterisation of the fundamental eigenfrequency is calculated in Fig. 10 in the adiabatic approximation and compared to the analytical estimate deduced from Eqs. (82), (83); for reference, it is also displayed in the flow with cooling, using the same adiabatic advection time, which should not be confused with the actual advection time taking into account the cooling layer. In the adiabatic flow, the winding angle of the spiral pattern of advected perturbations from the shock to the neutron star is ${\mathrm{\Phi }}_{\mathrm{adv}}\equiv 2\pi -{\phi }_{\mathcal{Q}}$. We note in Fig. 10 that φ_𝒬 < π, implying that the radial structure of advected perturbations of entropy or vorticity contains at least a change of sign, as seen in Fig. 2 of Blondin & Shaw (2007), Fig. 1 of Fernández (2010), and Fig. 10 of Buellet et al. (2023).

5.4. Comparison to the analytic model in Foglizzo (2009)

The adiabatic formulation incorporates major improvements compared to the adiabatic toy model used by F09, Sato et al. (2009), Guilet & Foglizzo (2012): the plane parallel model allows for explicit analytical solutions and a deeper understanding of the instability mechanism, but some simplifications are difficult to justify: (i) the assumption of a uniform flow between the shock and the region of deceleration implies that the vertical size of the acoustic and advected cavities are identical, thus overestimating the impact of radial acoustic propagation on the phase of the solution; (ii) the maximum strength 𝒬 ∼ 2 of the advective-acoustic cycle is arbitrarily set by the specific value (0.75) chosen for the ratio c_sh²/c_out², which limits the relative effect of the purely acoustic cycle on the most unstable modes; (iii) the optimum value of the vertical wavenumber compared to the horizontal one is set by the specific value chosen for L/H = 4 in F09 and Sato et al. (2009) and L/H = 6 in Guilet & Foglizzo (2012; iv) the width of the coupling region described by H_∇/H is also a free parameter of the model.

The present adiabatic model offers a new framework overcoming these shortcomings, and is still simple enough to allow for analytical results:

– The spherical geometry is taken into account (rather than a plane parallel approximation in F09).

– The adiabatic heating and deceleration are produced in a self-consistent manner by the radial convergence in 3D and by the proto-neutron star gravity (rather than being localised at the lower boundary by some external potential with adjustable depth c_sh²/c_out² and width H_∇/H in F09).

– The gravity produced by the neutron star at the shock is not neglected, allowing both displacement and velocity effects for the production of entropy at the shock (rather than ignoring displacement effects with ω_Φ = 0 in F09).

The adiabatic model sheds light on some objections by Blondin & Mezzacappa (2006) regarding the advective-acoustic mechanism in a non-rotating flow:

(i) The fact that the advection time may be longer than the oscillation period does not contradict the advective-acoustic mechanism.

(ii) The fact that the pressure field shows no evidence for a radially localised coupling radius is explained by the radially extended character of the coupling process, as seen in the integral in Eq. (63). The feedback should refer to the radially extended pressure structure rather than the radial propagation of an acoustic wave.

6. Conclusions and perspectives

1. The oscillation period of the fundamental SASI mode can be estimated from the shock radius r_sh, the radius of maximum deceleration r_∇, and the central mass M_ns.

2. A possible improvement to the analytic formula used by Müller & Janka (2014) is proposed, based on the perturbative analysis without any adjustment of a specific numerical simulation.

3. An adiabatic model of the shock dynamics incorporating non-adiabatic processes in the boundary conditions is able to capture the general properties of SASI eigenmodes.

4. In the adiabatic approximation, the SASI mechanism can be described as a self-forced oscillator. The forcing by advected perturbations is distributed from the shock surface to the cooling layer rather than inside the cooling layer.

5. An analytical description of the properties of the fundamental mode of SASI in the adiabatic approximation is obtained in the asymptotic limit of a large ratio r_sh/r_∇.

Improving the accuracy of the estimation of the SASI oscillation period from a perturbative analysis requires (i) a more accurate prescription for partial dissociation at the shock and in the flow, as described in Appendix A of Fernández & Thompson (2009a), and (ii) a better characterisation of the parametrised cooling function, leading to a realistic value of the ratio r_∇/r_ns.

The adiabatic framework proposed here opens a new path for analytical investigation of the physical effect of stellar rotation in the equatorial plane (Paper II), and provides a means to address the intriguing results of Walk et al. (2023), which they obtained by comparing SASI properties in cylindrical and spherical geometries.

We note that the adiabatic approximation indeed accounts for the growth rate and oscillation period within ∼30%, leaving room for additional non-adiabatic effects. By decreasing the velocity and Mach number near the proto-neutron star, the amplitude of the forcing term is increased, as is the phase mixing near the proto-neutron star. In addition to affecting the trade-off between these two effects, non-adiabatic effects modify the amplitude of δS and δK during their advection, add a source of feedback in the cooling layer, and modify the acoustic equation. A quantitative assessment of some of these effects can be obtained within the formalism of the self-forced oscillator by incorporating them into the lower boundary condition.

Acknowledgments

The anonymous referee is thanked for constructive suggestions. It is a pleasure to acknowledge stimulating discussions with Jérôme Guilet, Matteo Bugli, Thomas Janka, Bernhard Müller, Kei Kotake, Tomoya Takiwaki, Ernazar Abdikamalov, Laurie Walk, Irene Tamborra, Anne-Cécile Buellet, Sonia El Hedri, Florent Robinet and the LEAK collaboration funded by the LabEx UnivEarthS. This work has been financially supported by the PNHE.

References

Appendix A: Stationary accretion

The loss of energy through nuclear dissociation is measured by the parameter ε defined by Eq. (6). It decreases the value of the post-shock Mach number ℳ_sh below the reference adiabatic value ℳ_ad, and affects the Rankine-Hugoniot conditions according to Eqs. (A4-A6) of Foglizzo et al. (2006):

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{ad}}^2\equiv \frac{2+(\gamma -1){\mathcal{M}}_1^2}{2\gamma {\mathcal{M}}_1^2-\gamma +1},\end{array}$$(A.1)

$$\begin{array}{c}\hfill \epsilon =(1+\frac{2}{\gamma -1}\frac{1}{{\mathcal{M}}_1^2})(1-\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_{\mathrm{ad}}^2})(1-\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_1^2}),\end{array}$$(A.2)

$$\begin{array}{c}\hfill \frac{v_{\mathrm{sh}}}{v_1}=\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_1^2}\frac{1+\gamma {\mathcal{M}}_1^2}{1+\gamma {\mathcal{M}}_{\mathrm{sh}}^2},\end{array}$$(A.3)

$$\begin{array}{c}\hfill \frac{c_{\mathrm{sh}}}{c_1}=\frac{{\mathcal{M}}_{\mathrm{sh}}}{{\mathcal{M}}_1}\frac{1+\gamma {\mathcal{M}}_1^2}{1+\gamma {\mathcal{M}}_{\mathrm{sh}}^2}.\end{array}$$(A.4)

The nuclear binding energy of the iron nuclei is 8.8MeV per nucleon, or 1.77 × 10⁵²erg/M_⊙. If nuclear dissociation was complete across the shock this would translates into a dissociation parameter scaling linearly with r_sh:

$$\begin{array}{cc}\hfill \epsilon & =0.7\left(\frac{r_{\mathrm{sh}}}{150\mathrm{km}}\right)\left(\frac{1.3M_{\odot }}{M_{\mathrm{ns}}}\right).\hfill \end{array}$$(A.5)

We note that our definition of ε in Eq. (6) follows Huete et al. (2018), which differs from Eq. (4) in Fernández & Thompson (2009b) by a factor 2. A fraction of the dissociation of iron nuclei is radially distributed between the shock and the neutron star surface, as calculated in Fernández & Thompson (2009a). Taking into account incomplete dissociation, the value of ϵ formulated in Eq. (35) is a rough estimate deduced from Eq. (51) and Fig. 12 in Huete et al. (2018) for r_sh < 175km. It is smaller than Eq. (A.5) by a factor ∼1.4 and saturates at ϵ ∼ 0.5 in exploding models.

The value of ℳ_sh is expressed as a function of ℳ₁ and ε(r_sh) using Eq. (A.2):

$$\begin{array}{cc}\hfill \frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_{\mathrm{ad}}^2}& =\frac{1}{2}+\frac{1-\frac{2\epsilon }{1+\frac{2}{\gamma -1}\frac{1}{{\mathcal{M}}_1^2}}-\frac{{\mathcal{M}}_{\mathrm{ad}}^2}{2{\mathcal{M}}_1^2}}{1+{[{(1-\frac{{\mathcal{M}}_{\mathrm{ad}}^2}{{\mathcal{M}}_1^2})}^2+\frac{4\epsilon \frac{{\mathcal{M}}_{\mathrm{ad}}^2}{{\mathcal{M}}_1^2}}{1+\frac{2}{\gamma -1}\frac{1}{{\mathcal{M}}_1^2}}]}^{\frac{1}{2}}},\hfill \end{array}$$(A.6)

$$\begin{array}{cc}& \sim 1-\epsilon \mathrm{if}{\mathcal{M}}_1\gg 1\hfill \end{array}$$(A.7)

The jump condition (A.3) is thus a function of ε(r_sh) using Eq. (A.6), with the simple following expression for a strong shock:

$$\begin{array}{ccc}\hfill \frac{v_1}{v_{\mathrm{sh}}}& \sim \hfill & \hfill 1+\frac{2}{(1-\epsilon )(\gamma -1)},\end{array}$$(A.8)

$$\begin{array}{ccc}& \sim \hfill & \hfill 1+\frac{6}{1-\epsilon }.\end{array}$$(A.9)

The definition of the cooling function (1), the dimensionless entropy (2) and the mass conservation (3) are used to eliminate the variables $\rho ,p,\mathcal{M},\mathit{v}$ from the stationary equations, Eqs. (4) and (5):

$$\begin{array}{cc}\hfill \frac{\rho }{{\rho }_{\mathrm{sh}}}& ={\left(\frac{c}{c_{\mathrm{sh}}}\right)}^{\frac{2}{\gamma -1}}{\mathrm{e}}^{-(S-S_{\mathrm{sh}})},\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill r^g{c}^{\frac{\gamma +1}{\gamma -1}}\mathcal{M}& =r_{\mathrm{sh}}^g{c}_{\mathrm{sh}}^{\frac{\gamma +1}{\gamma -1}}{\mathrm{e}}^{S-S_{\mathrm{sh}}}.\hfill \end{array}$$(A.11)

Appendix B: Differential system and boundary conditions of the perturbed accretion

We rewrite the differential system satisfied by δf, δh (Eqs. (E2-E3) in Foglizzo et al. (2006)), using the radial coordinate X defined as in Foglizzo (2001):

$$\begin{array}{ccc}\hfill \mathrm{d}\mathrm{X}& \equiv \hfill & \hfill \frac{v}{1-{\mathcal{M}}^2}\mathrm{d}r.\end{array}$$(B.1)

Thus

$$\begin{array}{c}\hfill (\frac{\partial }{\partial X}+\frac{i\omega }{c^2})\frac{\delta f}{i\omega }=\delta h+(\gamma -1+\frac{1}{{\mathcal{M}}^2})\frac{\delta S}{\gamma }\\ \hfill +\frac{1-{\mathcal{M}}^2}{i\omega v}\delta \left(\frac{\mathcal{L}}{\rho v}\right),\end{array}$$(B.2)

$$\begin{array}{c}\hfill (\frac{\partial }{\partial X}+\frac{i\omega }{c^2})\delta h=\frac{i\omega }{v^2}(\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{{\omega }^2{c}^2}\delta f-\delta S)+\frac{1-{\mathcal{M}}^2}{v^2}\frac{i\delta K}{\omega r^2}\end{array}$$(B.3)

The set of equations Eqs. (A1)-(A4) in Foglizzo et al. (2007) are repeated here for completeness:

$$\begin{array}{cc}\hfill \frac{\delta v_r}{v}& =\frac{1}{1-{\mathcal{M}}^2}(\delta h+\delta S-\frac{\delta f}{c^2}),\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill \frac{\delta \rho }{\rho }& =\frac{1}{1-{\mathcal{M}}^2}(-{\mathcal{M}}^2\delta h-\delta S+\frac{\delta f}{c^2}),\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill \frac{\delta c^2}{c^2}& =\frac{\gamma -1}{1-{\mathcal{M}}^2}(\frac{\delta f}{c^2}-{\mathcal{M}}^2\delta h-{\mathcal{M}}^2\delta S),\hfill \end{array}$$(B.6)

$$\begin{array}{cc}\hfill \frac{\delta p}{\gamma p}& =\frac{1}{\gamma -1}\frac{\delta c^2}{c^2}-\frac{\delta S}{\gamma },\hfill \end{array}$$(B.7)

$$\begin{array}{cc}\hfill \delta \left(\frac{\mathcal{L}}{\rho v}\right)& =\mathrm{\nabla }S\frac{c^2}{\gamma }[(\beta -1)\frac{\delta \rho }{\rho }+\alpha \frac{\delta c^2}{c^2}-\frac{\delta v_r}{v}],\hfill \end{array}$$(B.8)

$$\begin{array}{cc}\hfill \delta \left(\frac{\mathcal{L}}{\mathit{pv}}\right)& =\mathrm{\nabla }S[(\beta -1)\frac{\delta \rho }{\rho }+(\alpha -1)\frac{\delta c^2}{c^2}-\frac{\delta v_r}{v}].\hfill \end{array}$$(B.9)

The perturbed mass conservation and transverse components of the perturbed Euler equation are:

$$\begin{array}{c}\hfill -i\omega \delta v_{\theta }+v\delta w_{\phi }+\frac{1}{r}\frac{\partial }{\partial \theta }\delta f=\frac{c^2}{\gamma }\frac{1}{r}\frac{\partial \delta S}{\partial \theta },\end{array}$$(B.10)

$$\begin{array}{c}\hfill -i\omega \delta v_{\phi }-v\delta w_{\theta }+\frac{\mathit{im}}{rsin\theta }\delta f=\frac{c^2}{\gamma }\frac{im\delta S}{rsin\theta },\end{array}$$(B.11)

We eliminate δv_θ and δv_φ in the system of Eqs. (12), (14), (B.10) and (B.11) and obtain Eq. (16).

The derivative of δA is calculated using Eqs. (16), (20) and (B.2), leading to the differential system (17-20) satisfied by δA, δh, δS, δK. The second derivative of δA is calculated using Eqs. (12), (20) and (B.3), resulting in Eq. (27).

Using Eq. (16), Eqs. (B.4-B.7) are rewritten as follows:

$$\begin{array}{cc}\hfill \frac{\delta v_r}{v}& =\frac{\delta K}{\ell (\ell +1)v^2}-\frac{1}{\ell (\ell +1)v}\frac{\partial \delta A}{\partial r}-\frac{\delta S}{\gamma {\mathcal{M}}^2},\hfill \end{array}$$(B.12)

$$\begin{array}{cc}\hfill \frac{\delta \rho }{\rho }& =\frac{1}{c^2}(v\frac{\partial }{\partial r}-i\omega )\frac{\delta A}{\ell (\ell +1)}-\frac{\gamma -1}{\gamma }\delta S,\hfill \end{array}$$(B.13)

$$\begin{array}{cc}\hfill \frac{1}{\gamma -1}\frac{\delta c^2}{c^2}& =\frac{1}{c^2}(v\frac{\partial }{\partial r}-i\omega )\frac{\delta A}{\ell (\ell +1)}+\frac{\delta S}{\gamma },\hfill \end{array}$$(B.14)

$$\begin{array}{cc}\hfill \frac{\delta p}{\gamma p}& =\frac{1}{c^2}(v\frac{\partial }{\partial r}-i\omega )\frac{\delta A}{\ell (\ell +1)},\hfill \end{array}$$(B.15)

Using Eq. (B.12) with Eq. (12), we can express δv_r with δA and δw_⊥ and obtain Eq. (15). We note that this relation is equivalent to the radial component of the vector calculus relation ∇²v = ∇(∇ ⋅ v)−∇×(∇ × v).

The baroclinic production of vorticity by the advection of entropy perturbations can be calculated in the same way as Eqs. (E17-E19) in Foglizzo et al. (2005), using the conservation of the tangential component of the velocity across the shock:

$$\begin{array}{cc}\hfill \delta v_{\theta \mathrm{sh}}& =\frac{v_1-v_{\mathrm{sh}}}{r_{\mathrm{sh}}}\frac{\partial \mathrm{\Delta }\zeta }{\partial \theta },\hfill \end{array}$$(B.16)

$$\begin{array}{cc}\hfill \delta v_{\phi \mathrm{sh}}& =\frac{v_1-v_{\mathrm{sh}}}{r_{\mathrm{sh}}}\frac{im\mathrm{\Delta }\zeta }{sin\theta }.\hfill \end{array}$$(B.17)

The vorticity produced at the shock is defined by Eqs. (B.16-B.17), the transverse components (B.10-B.11) of the Euler equation at the shock, together with the angular derivative of Eqs. (22):

$$\begin{array}{cc}\hfill \delta w_{r\mathrm{sh}}& =0,\hfill \end{array}$$(B.18)

$$\begin{array}{cc}\hfill \delta w_{\theta \mathrm{sh}}& =-\frac{c_{\mathrm{sh}}^2}{\gamma }\frac{\mathit{im}}{r_{\mathrm{sh}}{v}_{\mathrm{sh}}sin\theta }(\delta S_{\mathrm{sh}}+{\left[\mathrm{\nabla }S\right]}_1^{\mathrm{sh}}\mathrm{\Delta }\zeta ),\hfill \end{array}$$(B.19)

$$\begin{array}{cc}\hfill \delta w_{\phi \mathrm{sh}}& =\frac{c_{\mathrm{sh}}^2}{\gamma }\frac{1}{r_{\mathrm{sh}}{v}_{\mathrm{sh}}}\frac{\partial }{\partial \theta }(\delta S_{\mathrm{sh}}+{\left[\mathrm{\nabla }S\right]}_1^{\mathrm{sh}}\mathrm{\Delta }\zeta ).\hfill \end{array}$$(B.20)

Appendix C: Adiabatic model

Following Eqs. (B5)-(B7) in Foglizzo (2001), the differential equation describing the specific vorticity in a spherical adiabatic flow can be integrated as

$$\begin{array}{cc}\hfill \delta w_r& ={\left(\frac{r_{\mathrm{sh}}}{r}\right)}^2\delta w_{r\mathrm{sh}}{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}r}{v}},\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill \delta w_{\theta }& =\frac{1}{\mathit{rv}}[{(rv\delta w_{\theta })}_{\mathrm{sh}}-\frac{c^2-c_{\mathrm{sh}}^2}{sin\theta }im\frac{\delta S_{\mathrm{sh}}}{\gamma }]{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}r}{v}},\hfill \end{array}$$(C.2)

$$\begin{array}{cc}\hfill \delta w_{\phi }& =\frac{1}{\mathit{rv}}[{(rv\delta w_{\phi })}_{\mathrm{sh}}+\frac{c^2-c_{\mathrm{sh}}^2}{sin\theta }\frac{\partial }{\partial \theta }\frac{\delta S_{\mathrm{sh}}}{\gamma }]{\mathrm{e}}^{i\omega {\int }_{\mathrm{sh}}\frac{\mathrm{d}r}{v}}.\hfill \end{array}$$(C.3)

Using Eqs. (B.18-B.20) with ∇S = 0, together with Eqs. (C.1-C.3) gives the expression (42-44) of the vorticity perturbation throughout the flow.

With δY defined by Eq. (59), the expressions of δA_sh, δh_sh are deduced from Eqs. (22-24):

$$\begin{array}{c}\hfill (1-\frac{v_{\mathrm{sh}}}{v_1})\mathrm{\Delta }\zeta =\frac{\delta Y_{\mathrm{sh}}}{v_1},\end{array}$$(C.4)

$$\begin{array}{c}\hfill \delta h_{\mathrm{sh}}=-\frac{i\omega }{v_1{v}_{\mathrm{sh}}}\delta Y_{\mathrm{sh}},\end{array}$$(C.5)

$$\begin{array}{c}\hfill \frac{\delta S_{\mathrm{sh}}}{\gamma }=\frac{\delta Y_{\mathrm{sh}}}{c_{\mathrm{sh}}^2}(i\omega +{\omega }_{\mathrm{\Phi }})(1-\frac{v_{\mathrm{sh}}}{v_1}).\end{array}$$(C.6)

We rewrite Eq. (17) using the definitions of X and δY with δK = 0:

$$\begin{array}{c}\hfill \delta Y_{\mathrm{sh}}=-\frac{\delta A_{\mathrm{sh}}}{\ell (\ell +1)},\end{array}$$(C.7)

$$\begin{array}{c}\hfill {\left(\frac{\partial \delta A}{\partial r}\right)}_{\mathrm{sh}}=\frac{\ell (\ell +1)v_{\mathrm{sh}}}{1-{\mathcal{M}}_{\mathrm{sh}}^2}[\delta Y_{\mathrm{sh}}\frac{i\omega }{v_{\mathrm{sh}}^2}(\frac{v_{\mathrm{sh}}}{v_1}+{\mathcal{M}}_{\mathrm{sh}}^2)\\ \hfill -\frac{\delta S_{\mathrm{sh}}}{\gamma }(\frac{1}{{\mathcal{M}}_{\mathrm{sh}}^2}+\gamma -1)],\end{array}$$(C.8)

$$\begin{array}{c}\hfill {\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{sh}}=-\frac{1-{\mathcal{M}}_{\mathrm{sh}}^2}{\ell (\ell +1)v_{\mathrm{sh}}}{\left(\frac{\partial \delta A}{\partial r}\right)}_{\mathrm{sh}}+\frac{i\omega }{c_{\mathrm{sh}}^2}\delta Y_{\mathrm{sh}},\end{array}$$(C.9)

which results in Eq. (57).

At the inner boundary in the adiabatic approximation, we use the condition δv_r = 0:

$$\begin{array}{c}\hfill \delta h_{\mathrm{ns}}+\delta S_{\mathrm{ns}}=\frac{\delta f}{c_{\mathrm{ns}}^2}.\end{array}$$(C.10)

The entropy is simply advected from the shock (38) and we express δh with δY and ∂δY/∂X using Eq. (17):

$$\begin{array}{cc}\hfill \delta S_{\mathrm{ns}}& =\delta S_{\mathrm{sh}}{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}i\omega \frac{\mathrm{d}r}{v}},\hfill \end{array}$$(C.11)

$$\begin{array}{cc}\hfill \frac{\partial \delta Y}{\partial X}& =\delta h{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{c^2}\mathrm{d}X}\hfill \\ \hfill & +\frac{\delta S}{\gamma }{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{c^2}\mathrm{d}X}(\frac{1}{{\mathcal{M}}^2}+\gamma -1).\hfill \end{array}$$(C.12)

The inner boundary condition (C.10) is thus reduced to Eq. (58).

Multiplying Eq. (51) by m/(ω sin θ) and using Eqs. (43) and (48):

$$\begin{array}{c}\hfill [{(\frac{\partial }{\partial X}+\frac{i\omega }{c^2})}^2+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}]r\delta v_{\phi }=-\frac{\partial }{\partial X}\frac{r\delta w_{\theta }}{v}.\end{array}$$(C.13)

Equivalently, using Eq. (44) and (47) and the derivative of Eq. (51) with respect to θ leads to:

$$\begin{array}{c}\hfill [{(\frac{\partial }{\partial X}+\frac{i\omega }{c^2})}^2+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2}]r\delta v_{\theta }=\frac{\partial }{\partial X}\frac{r\delta w_{\phi }}{v},\end{array}$$(C.14)

which is equivalent to Eq. (C.13) given the spherical symmetry of the stationary flow.

Appendix D: Equation defining the eigenfrequencies in the adiabatic model

The differential equation satisfied by δY is deduced from Eqs. (51) and (59):

$$\begin{array}{c}\hfill (\frac{{\partial }^2}{\partial X^2}+\frac{{\omega }^2-{\omega }_{\mathrm{Lamb}}^2}{v^2{c}^2})\delta Y=\delta \mathcal{F},\end{array}$$(D.1)

$$\begin{array}{c}\hfill \delta \mathcal{F}\equiv {\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{c^2}\mathrm{d}X}{\mathcal{F}}_S\delta S_{\mathrm{sh}}.\end{array}$$(D.2)

We define Y₀ as the solution of the homogeneous equation satisfying the inner boundary condition of pure acoustic waves (i.e. without entropy and vorticity perturbations), and Y₋ another homogeneous solution such that their Wronskien is W:

$$\begin{array}{c}\hfill {\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{ns}}=\frac{i\omega }{c_{\mathrm{ns}}^2}{Y}_0^{\mathrm{ns}},\end{array}$$(D.3)

$$\begin{array}{c}\hfill W\equiv Y_0\frac{\partial Y_{-}}{\partial X}-Y_{-}\frac{\partial Y_0}{\partial X},\end{array}$$(D.4)

$$\begin{array}{c}\hfill \delta Y=Y_{-}(d_{-}+{\int }_{\mathrm{ns}}\frac{Y_0}{W}\delta \mathcal{F}\mathrm{d}X)-Y_0(d_0+{\int }_{\mathrm{sh}}\frac{Y_{-}}{W}\delta \mathcal{F}\mathrm{d}X),\end{array}$$(D.5)

where $\delta \mathcal{F}\equiv {\mathcal{F}}_S\delta S_{\mathrm{sh}}$ is the forcing term on the right hand side of Eq. (51) for the variable δY defined by. Using Eq. (D.5) at the upper boundary:

$$\begin{array}{c}\hfill \delta Y_{\mathrm{sh}}=Y_{-}^{\mathrm{sh}}(d_{-}+{\int }_{\mathrm{ns}}^{\mathrm{sh}}\frac{Y_0}{W}\delta \mathcal{F}\mathrm{d}X)-d_0{Y}_0^{\mathrm{sh}},\end{array}$$(D.6)

$$\begin{array}{c}\hfill {\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{sh}}={\left(\frac{\partial Y_{-}}{\partial X}\right)}_{\mathrm{sh}}(d_{-}+{\int }_{\mathrm{ns}}^{\mathrm{sh}}\frac{Y_0}{W}\delta \mathcal{F}\mathrm{d}X)-d_0{\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{sh}},\end{array}$$(D.7)

Using Eq. (D.5) at the lower boundary:

$$\begin{array}{c}\hfill \delta Y_{\mathrm{ns}}=d_{-}{Y}_{-}^{\mathrm{ns}}-Y_0^{\mathrm{ns}}(d_0-{\int }_{\mathrm{ns}}^{\mathrm{sh}}\frac{Y_{-}}{W}\delta \mathcal{F}\mathrm{d}X),\end{array}$$(D.8)

$$\begin{array}{c}\hfill {\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{ns}}=d_{-}{\left(\frac{\partial Y_{-}}{\partial X}\right)}_{\mathrm{ns}}-{\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{ns}}(d_0-{\int }_{\mathrm{ns}}^{\mathrm{sh}}\frac{Y_{-}}{W}\delta \mathcal{F}\mathrm{d}X).\end{array}$$(D.9)

Using Eq. (D.3), the lower boundary condition (58) translates into

$$\begin{array}{c}\hfill d_{-}[{\left(\frac{\partial Y_{-}}{\partial X}\right)}_{\mathrm{ns}}-\frac{i\omega }{c_{\mathrm{ns}}^2}{Y}_{-}^{\mathrm{ns}}]=\delta {\mathcal{F}}_{\mathrm{sh}}(1-{\mathcal{M}}_{\mathrm{ns}}^2)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_{\mathrm{ns}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{i\omega }{v^2}\mathrm{d}X}.\end{array}$$(D.10)

We note using Eq. (D.4) and (D.3) that

$$\begin{array}{cc}\hfill {\left(\frac{\partial Y_{-}}{\partial X}\right)}_{\mathrm{ns}}-\frac{i\omega }{c_{\mathrm{ns}}^2}{Y}_{-}^{\mathrm{ns}}& =\frac{W}{Y_0^{\mathrm{ns}}}.\hfill \end{array}$$(D.11)

Thus

$$\begin{array}{c}\hfill W{d}_{-}=Y_0^{\mathrm{ns}}\delta {\mathcal{F}}_{\mathrm{sh}}(1-{\mathcal{M}}_{\mathrm{ns}}^2)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_{\mathrm{ns}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{i\omega }{v^2}\mathrm{d}X}.\end{array}$$(D.12)

Eliminating d₀ between Eqs. (D.6) and (D.7) and using Eq. (D.4):

$$\begin{array}{c}\hfill Y_0^{\mathrm{sh}}{\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{sh}}-{\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{sh}}\delta Y_{\mathrm{sh}}=\\ \hfill [Y_0^{\mathrm{sh}}{\left(\frac{\partial Y_{-}}{\partial X}\right)}_{\mathrm{sh}}-{\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{sh}}{Y}_{-}^{\mathrm{sh}}](d_{-}+{\int }_{\mathrm{ns}}^{\mathrm{sh}}\frac{Y_0}{W}\delta \mathcal{F}\mathrm{d}X),\end{array}$$(D.13)

$$\begin{array}{c}\hfill =W{d}_{-}+{\int }_{\mathrm{ns}}^{\mathrm{sh}}{Y}_0\delta \mathcal{F}\mathrm{d}X.\end{array}$$(D.14)

The eigenfrequencies are thus defined by

$$\begin{array}{c}\hfill Y_0^{\mathrm{sh}}{\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{sh}}-{\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{sh}}\delta Y_{\mathrm{sh}}=\\ \hfill Y_0^{\mathrm{ns}}\delta {\mathcal{F}}_{\mathrm{sh}}(1-{\mathcal{M}}_{\mathrm{ns}}^2)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_{\mathrm{ns}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{i\omega }{v^2}\mathrm{d}X}+{\int }_{\mathrm{ns}}^{\mathrm{sh}}{Y}_0\delta \mathcal{F}\mathrm{d}X.\end{array}$$(D.15)

Replacing the forcing term by its expression (D.2),

$$\begin{array}{c}\hfill Y_0^{\mathrm{sh}}{\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{sh}}-{\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{sh}}\delta Y_{\mathrm{sh}}=\\ \hfill \delta {\mathcal{F}}_{\mathrm{sh}}\{Y_0^{\mathrm{ns}}(1-{\mathcal{M}}_{\mathrm{ns}}^2)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_{\mathrm{ns}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{i\omega }{v^2}\mathrm{d}X}\\ \hfill +{\int }_{\mathrm{ns}}^{\mathrm{sh}}{Y}_0{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{c^2}\mathrm{d}X}\frac{\partial }{\partial r}\left(\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{v}\mathrm{d}r}\right)\mathrm{d}r\}.\end{array}$$(D.16)

Using Eqs. (24) and (57) this equation takes the form

$$\begin{array}{c}\hfill a_1{Y}_0^{\mathrm{sh}}+a_2{r}_{\mathrm{sh}}{\left(\frac{\partial Y_0}{\partial r}\right)}_{\mathrm{sh}}+a_3{Y}_0^{\mathrm{ns}}=\\ \hfill {\int }_{\mathrm{ns}}^{\mathrm{sh}}{Y}_0{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{c^2}\mathrm{d}X}\frac{\partial }{\partial r}\left(\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{v}\mathrm{d}r}\right)\mathrm{d}r,\end{array}$$(D.17)

with a₁, a₂, a₃ defined by:

$$\begin{array}{ccc}\hfill a_1& \equiv \hfill & \hfill \frac{{\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{sh}}}{\delta {\mathcal{F}}_{\mathrm{sh}}},\end{array}$$(D.18)

$$\begin{array}{cc}& =1+(\gamma -1){\mathcal{M}}_{\mathrm{sh}}^2-\frac{i\omega }{i\omega +{\omega }_{\mathrm{\Phi }}}\frac{\frac{v_{\mathrm{sh}}}{v_1}}{1-\frac{v_{\mathrm{sh}}}{v_1}},\hfill \end{array}$$(D.19)

$$\begin{array}{ccc}\hfill a_2& \equiv \hfill & \hfill -\frac{1-{\mathcal{M}}_{\mathrm{sh}}^2}{r_{\mathrm{sh}}{v}_{\mathrm{sh}}}\frac{\delta Y_{\mathrm{sh}}}{\delta {\mathcal{F}}_{\mathrm{sh}}},\end{array}$$(D.20)

$$\begin{array}{cc}& =-\frac{1-{\mathcal{M}}_{\mathrm{sh}}^2}{(1-\frac{v_{\mathrm{sh}}}{v_1})\frac{(i\omega +{\omega }_{\mathrm{\Phi }})r_{\mathrm{sh}}}{v_{\mathrm{sh}}}},\hfill \end{array}$$(D.21)

$$\begin{array}{ccc}\hfill a_3& \equiv \hfill & \hfill -(1-{\mathcal{M}}_{\mathrm{ns}}^2)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}_{\mathrm{ns}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{i\omega }{v^2}\mathrm{d}X}.\end{array}$$(D.22)

After one integration by parts:

$$\begin{array}{c}\hfill Y_0^{\mathrm{sh}}{\left(\frac{\partial \delta Y}{\partial X}\right)}_{\mathrm{sh}}-{\left(\frac{\partial Y_0}{\partial X}\right)}_{\mathrm{sh}}\delta Y_{\mathrm{sh}}=\delta {\mathcal{F}}_{\mathrm{sh}}\{Y_0^{\mathrm{sh}}-Y_0^{\mathrm{ns}}{\mathcal{M}}_{\mathrm{sh}}^2{\mathrm{e}}^{{\int }_{\mathrm{sh}}^{\mathrm{ns}}\frac{i\omega }{v^2}\mathrm{d}X}\\ \hfill -{\int }_{\mathrm{ns}}^{\mathrm{sh}}\frac{\partial }{\partial r}\left(Y_0{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{c^2}\mathrm{d}X}\right)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\mathcal{M}}^2}{\mathrm{e}}^{{\int }_{\mathrm{sh}}\frac{i\omega }{v}\mathrm{d}r}\mathrm{d}r\}.\end{array}$$(D.23)

The equation defining the eigenfrequencies becomes Eq. (63) with a′₁ ≡ a₁ − 1 and a′₂ ≡ a₂ defined by Eqs. (64-65).

Appendix E: Approximation of the adiabatic stationary flow and the homogeneous perturbative solution

The dissociation measured by the parameter ε affects the relation between |v_sh| and the local free fall velocity, as deduced from Eq. (68) for a strong shock:

$$\begin{array}{ccc}\hfill |v_{\mathrm{sh}}|{\left(\frac{r_{\mathrm{sh}}}{2G{M}_{\mathrm{ns}}}\right)}^{\frac{1}{2}}& \sim \hfill & \hfill \frac{(\gamma -1)(1-\epsilon )}{2+(\gamma -1)(1-\epsilon )},\end{array}$$(E.1)

$$\begin{array}{cc}\hfill \frac{c_{\mathrm{sh}}^2{r}_{\mathrm{sh}}}{G{M}_{\mathrm{ns}}}& =\frac{1}{{\mathcal{M}}_{\mathrm{sh}}^2}\frac{v_{\mathrm{sh}}^2{r}_{\mathrm{sh}}}{G{M}_{\mathrm{ns}}}.\hfill \end{array}$$(E.2)

Equation (E.2) is transformed into Eq. (69) using Eqs. (E.1) and (67). We use a power law approximation of the enthalpy profile deduced from the Bernoulli equation for r ≪ r_sh and ℳ² ≪ 1:

$$\begin{array}{c}\hfill c^2[1+(\gamma -1)\frac{{\mathcal{M}}^2}{2}]=(\gamma -1)\frac{G{M}_{\mathrm{ns}}}{r}\\ \hfill +c_{\mathrm{sh}}^2[1+(\gamma -1)\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{2}]-(\gamma -1)\frac{G{M}_{\mathrm{ns}}}{r_{\mathrm{sh}}},\end{array}$$(E.3)

$$\begin{array}{c}\hfill c^2\sim (\gamma -1)\frac{G{M}_{\mathrm{ns}}}{r}.\end{array}$$(E.4)

The mass conservation and the adiabatic hypothesis in Eq. (A.11) imply the power law approximations (71) and (72) for the velocity and Mach number profiles.

For γ = 4/3 we approximate

$$\begin{array}{c}\hfill {\mathcal{M}}^2\sim {\mathcal{M}}_{\mathrm{sh}}^2{\left(\frac{r}{r_{\mathrm{sh}}}\right)}^3,\end{array}$$(E.5)

$$\begin{array}{c}\hfill v\sim v_{\mathrm{sh}}\frac{r}{r_{\mathrm{sh}}},\end{array}$$(E.6)

$$\begin{array}{c}\hfill i\omega {\int }_{\mathrm{sh}}\frac{{\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v}\sim -{\mathcal{M}}_{\mathrm{sh}}^2\frac{i\omega r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}{\int }_{\mathrm{sh}}\frac{x^2\mathrm{d}x}{1-{\mathcal{M}}_{\mathrm{sh}}^2{x}^3},\end{array}$$(E.7)

$$\begin{array}{c}\hfill \sim \frac{i\omega r_{\mathrm{sh}}}{3|v_{\mathrm{sh}}|}log\left(\frac{1-{\mathcal{M}}^2}{1-{\mathcal{M}}_{\mathrm{sh}}^2}\right).\end{array}$$(E.8)

The oscillatory phase associated to δf in Eq. (63) is made of a product of (iω) with two explicit contributions and a third contribution from the definition of δY in Eq. (59). The sum of these three contributions is:

$$\begin{array}{c}\hfill {\int }_{\mathrm{sh}}\frac{{\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v}+{\int }_{\mathrm{sh}}\frac{{\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v}+{\int }_{\mathrm{sh}}\frac{\mathrm{d}r}{v}\\ \hfill ={\int }_{\mathrm{sh}}\frac{1+{\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v}.\end{array}$$(E.9)

With α_v = 1 and $\stackrel{\sim }{r}\equiv r/r_{\mathrm{sh}}$,

$$\begin{array}{c}\hfill \frac{1+{\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{v_{\mathrm{sh}}}{v}\sim \frac{2{\mathcal{M}}_{\mathrm{sh}}^2{\tilde{r}}^2}{1-{\mathcal{M}}_{\mathrm{sh}}^2{\tilde{r}}^3}+\frac{1}{\tilde{r}}.\end{array}$$(E.10)

The integral can be estimated as follows:

$$\begin{array}{ccc}\hfill {\int }_{\mathrm{sh}}\frac{1+{\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v}& \sim \hfill & \hfill \frac{r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}(log\frac{r_{\mathrm{sh}}}{r_{\mathrm{ns}}}-\frac{2}{3}log\frac{1-{\mathcal{M}}_{\mathrm{sh}}^2}{1-{\mathcal{M}}_{\mathrm{ns}}^2}),\end{array}$$(E.11)

$$\begin{array}{ccc}& \sim \hfill & \hfill \frac{r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}(log\frac{r_{\mathrm{sh}}}{r_{\mathrm{ns}}}+\frac{2{\mathcal{M}}_{\mathrm{sh}}^2}{3})\end{array}$$(E.12)

With a strong adiabatic shock ${\mathcal{M}}_{\mathrm{sh}}^2\sim 1/8$, the correction to the advection timescale ${\tau }_{\mathrm{adv}}^{\mathrm{ns}}$ is thus of order 12% for r_sh/r_ns = 2 and 5% for r_sh/r_ns = 5. This correction is smaller if dissociation diminishes the value of ${\mathcal{M}}_{\mathrm{sh}}^2$.

The contribution of the integral inside the radial derivative in Eq. (63) is limited to a contribution of order ${\mathcal{M}}_{\mathrm{sh}}^2$:

$$\begin{array}{cc}\hfill {\int }_{\mathrm{sh}}^r\omega \frac{\mathrm{d}X}{c^2}& =\omega {\int }_{\mathrm{sh}}^r\frac{{\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v},\hfill \end{array}$$(E.13)

$$\begin{array}{ccc}& \sim \hfill & \hfill \frac{\omega r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{{\alpha }_v+2}[1-{\left(\frac{r}{r_{\mathrm{sh}}}\right)}^{{\alpha }_v+2}].\end{array}$$(E.14)

With ω_r ∼ 2π|v_sh|/r_sh for a small shock radius this phase shift is not negligible since it reaches 0.74 ∼ π/4 at the inner boundary and it is linear in ω.

We integrate Eq. (49) using Eq. (71) and neglecting ℳ² ≪ 1:

$$\begin{array}{ccc}\hfill X_{\mathrm{sh}}& \equiv \hfill & \hfill \frac{r_{\mathrm{sh}}{v}_{\mathrm{sh}}}{{\alpha }_v+1},\end{array}$$(E.15)

$$\begin{array}{ccc}\hfill \frac{X}{X_{\mathrm{sh}}}& \sim \hfill & \hfill {\left(\frac{r}{r_{\mathrm{sh}}}\right)}^{{\alpha }_v+1}.\end{array}$$(E.16)

The two solutions Y₀^± of the homogeneous equation (60) are approximated as power laws with exponents α_±.

$$\begin{array}{ccc}\hfill Y_0^{\pm }(r)& \equiv \hfill & \hfill {\left(\frac{r}{r_{\mathrm{ns}}}\right)}^{{\alpha }_{\pm }}.\end{array}$$(E.17)

Injecting Y₀^±(r) into Eq. (60),

$$\begin{array}{c}\hfill {\alpha }_{\pm }({\alpha }_{\pm }-{\alpha }_v-1){\left(\frac{r}{r_{\mathrm{sh}}}\right)}^{-2{\alpha }_v-2}=\\ \hfill -[\frac{{\omega }^2{r}_{\mathrm{sh}}^2}{c^2}-\ell (\ell +1)\frac{r_{\mathrm{sh}}^2}{r^2}]{\left(\frac{v_{\mathrm{sh}}}{v}\right)}^2.\end{array}$$(E.18)

For ℓ ≥ 1, the restoring force in Eq. (60) is independent of the frequency in the region where ω ≪ ω_Lamb. Thus, using Eq. (E.4) for ℓ ≥ 1,

$$\begin{array}{cc}\hfill {\alpha }_{\pm }({\alpha }_{\pm }-{\alpha }_v-1)& =\ell (\ell +1)-{\mathcal{M}}_{\mathrm{sh}}^2{\left(\frac{\omega r_{\mathrm{sh}}}{v_{\mathrm{sh}}}\right)}^2{\left(\frac{r}{r_{\mathrm{sh}}}\right)}^3,\hfill \end{array}$$(E.19)

$$\begin{array}{ccc}& \sim \hfill & \hfill \ell (\ell +1)\mathrm{if}\omega \ll {\omega }_{\mathrm{Lamb}}.\end{array}$$(E.20)

In this region the approximate solution is α_± ≡ a ± b with a, b defined by Eqs. (75) and (76). The homogeneous solution Y₀ for ℓ ≥ 1 is a linear combination of power laws Y₀^± satisfying the lower boundary condition (61):

$$\begin{array}{c}\hfill Y_0(r)={\left(\frac{r}{r_{\mathrm{ns}}}\right)}^{a-b}+{\mathcal{R}}_{\mathrm{ns}}{\left(\frac{r}{r_{\mathrm{ns}}}\right)}^{a+b},\end{array}$$(E.21)

$$\begin{array}{c}\hfill \frac{1-{\mathcal{M}}_{\mathrm{ns}}^2}{v_{\mathrm{ns}}}{\left(\frac{\partial Y_0}{\partial r}\right)}_{\mathrm{ns}}=\frac{i\omega }{c_{\mathrm{ns}}^2}{Y}_0^{\mathrm{ns}}.\end{array}$$(E.22)

With ℳ_ns ≪ 1 and using Eq. (72),

$$\begin{array}{c}\hfill a-b+(a+b){\mathcal{R}}_{\mathrm{ns}}={\mathcal{M}}_{\mathrm{sh}}^2{\left(\frac{r_{\mathrm{ns}}}{r_{\mathrm{sh}}}\right)}^{2+{\alpha }_v}\frac{i\omega r_{\mathrm{sh}}}{v_{\mathrm{sh}}}(1+{\mathcal{R}}_{\mathrm{ns}}),\end{array}$$(E.23)

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{ns}}=\frac{a-b-{\mathcal{M}}_{\mathrm{sh}}^{\frac{4}{\gamma +1}}{\left(\frac{r_{\mathrm{ns}}}{r_{\mathrm{sh}}}\right)}^{\frac{\gamma +3+{\alpha }_v(3-\gamma )}{\gamma +1}}\frac{i\omega r_{\mathrm{sh}}}{v_{\mathrm{sh}}}}{a+b+{\mathcal{M}}_{\mathrm{sh}}^{\frac{4}{\gamma +1}}{\left(\frac{r_{\mathrm{ns}}}{r_{\mathrm{sh}}}\right)}^{\frac{\gamma +3+{\alpha }_v(3-\gamma )}{\gamma +1}}\frac{i\omega r_{\mathrm{sh}}}{v_{\mathrm{sh}}}}.\end{array}$$(E.24)

From Eq. (E.24) we conclude that the coefficient ℛ_ns is asymptotically independent of ω when r_ns ≪ r_sh.

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{ns}}\sim -\frac{a-b}{a+b}.\end{array}$$(E.25)

The acoustic function is independent of the frequency in the asymptotic limit r_ns ≪ r_sh for low frequency perturbations ℓ ≥ 1 driven by advection such that ωr_sh/|v_sh|≤1, as described by Eqs. (73) and (74).

Using the index 0 to denote the acoustic perturbation associated to the homogeneous solution, we note from Eqs. (46) with δS = 0 and Eq. (59) that the relation between ∂Y₀/∂r and δv_r⁰ is

$$\begin{array}{cc}\hfill \delta v_r^0& =-\frac{1}{\ell (\ell +1)}\frac{\partial \delta A_0}{\partial r},\hfill \end{array}$$(E.26)

$$\begin{array}{cc}& =(\frac{\partial Y_0}{\partial r}-\frac{i\omega {\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{Y_0}{v}){\mathrm{e}}^{-{\int }_{\mathrm{sh}}\frac{i\omega {\mathcal{M}}^2}{1-{\mathcal{M}}^2}\frac{\mathrm{d}r}{v}}.\hfill \end{array}$$(E.27)

The fact that Eq. (74) imposes ∂Y₀/∂r = 0 is approximately compatible with δv_r⁰ ∼ 0 to the extent that ℳ_ns ≪ 1.

Appendix F: Approximation of the eigenfrequency equation for γ = 4/3

The advection time for α_v = 1 is approximated using Eq. (79):

$$\begin{array}{c}\hfill {\mathrm{e}}^{{\int }_{\mathrm{sh}}^r\frac{i\omega }{v_r}\mathrm{d}r}={\left(\frac{x}{x_{\mathrm{sh}}}\right)}^{-\frac{i\omega r_{\mathrm{sh}}}{|v_{\mathrm{sh}}|}}.\end{array}$$(F.1)

Defining δz ≡ iωr_sh/|v_sh|+2 − b and noting that a = 1, the leading terms in Eq. (81) for ℓ ≥ 1 are :

$$\begin{array}{c}\hfill N(\delta z-2+b)+\frac{2b{\mathcal{M}}_{\mathrm{sh}}^2}{\ell (\ell +1)}{x}_{\mathrm{sh}}^{\delta z-3}={\int }_1^{x_{\mathrm{sh}}}{\left(\frac{x_{\mathrm{sh}}}{x}\right)}^{\delta z}(\frac{1}{x^{2b}}-1)\frac{\mathrm{d}x}{x},\end{array}$$(F.2)

The integral on the right hand side can be calculated explicitly. If δz ≠ 0,

$$\begin{array}{c}\hfill {\int }_1^{x_{\mathrm{sh}}}{\left(\frac{x_{\mathrm{sh}}}{x}\right)}^{\delta z}(\frac{1}{x^{2b}}-1)\frac{\mathrm{d}x}{x}=\frac{1}{\delta z}(1-\frac{\delta z{x}_{\mathrm{sh}}^{-2b}}{2b+\delta z}-\frac{2b{x}_{\mathrm{sh}}^{\delta z}}{2b+\delta z}),\end{array}$$(F.3)

Using Eq. (F.3) in Eq. (F.2) we obtain the following equation defining δz:

$$\begin{array}{c}\hfill [1-\delta zN(b-2+\delta z)](1+\frac{\delta z}{2b})-\frac{\delta z}{2b{x}_{\mathrm{sh}}^{2b}}=\\ \hfill x_{\mathrm{sh}}^{\delta z}[1+\delta z\frac{2b{\mathcal{M}}_{\mathrm{sh}}^2}{\ell (\ell +1)x_{\mathrm{sh}}^3}(1+\frac{\delta z}{2b})].\end{array}$$(F.4)

Reintroducing the normalised eigenfrequency Z = iωr_sh/|v_sh|,

$$\begin{array}{c}\hfill \frac{x_{\mathrm{sh}}^{b-2}}{2b}\{[1-(Z+2-b)N\left(Z\right)](Z+2+b)-\frac{Z+2-b}{x_{\mathrm{sh}}^{2b}}\}=\\ \hfill x_{\mathrm{sh}}^Z\{1+[{(Z+2)}^2-b^2]\frac{{\mathcal{M}}_{\mathrm{sh}}^2}{\ell (\ell +1)x_{\mathrm{sh}}^3}\}.\end{array}$$(F.5)

The approximate advection time ${\tau }_{\mathrm{adv}}^{\mathrm{sh}}$ from the shock to the inner boundary, defined by Eq. (79), is introduced to obtain Eqs. (82-83).

All Figures