# Analytic insight into the physics of SASI I. Shock instability in a non-rotating stellar core

T Foglizzo

► To cite this version:

T Foglizzo. Analytic insight into the physics of SASI I. Shock instability in a non-rotating stellar core. Astronomy & Astrophysics - A&A, 2024, 692, pp.A196. 10.1051/0004-6361/202452300. hal-04710273

HAL Id: hal-04710273

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

Submitted on 16 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Analytic insight into the physics of the standing accretion shock instability

## I. Shock instability in a non-rotating stellar core

T. Foglizzo*

Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

Received 18 September 2024 / Accepted 23 October 2024

### ABSTRACT

**Context.** During the core collapse of a massive star, and immediately before its supernova explosion, there is amplification of asymmetric motions by the standing accretion shock instability (SASI). This imprints a frequency signature on the neutrino flux and the gravitational waves that carries direct information about the explosion process.

**Aims.** The physical interpretation of this multi-messenger signature requires a detailed understanding of the instability mechanism.

**Methods.** We carried out a perturbative analysis to characterise the properties of SASI and assess the effect of the region of neutronization above the surface of the proto-neutron star. We compared the eigenfrequencies of the most unstable modes to those obtained in an adiabatic approximation where neutrino interactions are neglected above the neutrinosphere. We solved the differential system analytically using a Wronskian method and approximated it asymptotically for a large shock radius.

**Results.** The oscillation period of SASI is well fitted with a simple analytic function of the shock radius, the radius of maximum deceleration, and the mass of the proto-neutron star. The oscillation period is weakly dependent on the parameterised cooling function, but this latter does affects the SASI growth rate. We describe the general properties of SASI eigenmodes using an adiabatic model. In this approximation, the eigenvalue problem is formulated as a self-forced oscillator. The forcing agent is the radial advection of baroclinic vorticity perturbations and entropy perturbations produced by the shock oscillation. We reduced the differential system defining the eigenfrequencies to a single integral equation. Its analytical approximation sheds light on the radially extended character of the region of advective-acoustic coupling. The simplicity of this adiabatic formalism opens new perspectives for the investigation of the effect of stellar rotation and non-adiabatic processes on SASI.

**Key words.** hydrodynamics – instabilities – neutrinos – shock waves – stars: neutron – supernovae: general

### 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

<footnote>* Corresponding author; foglizzo@cea.fr</footnote>
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_\odot$ 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 $\rho$ and pressure $p$:

$$ \mathcal{L} = -A_{\text{cool}}\rho^{\beta-\alpha}p^{\alpha}. \qquad (1) $$

The collapsing stellar core immediately after bounce is modelled as a perfect gas with an adiabatic index of $\gamma = 4/3$, dominated by the degeneracy pressure of relativistic electrons. The gravitational potential $\Phi \equiv -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

$$ S \equiv \frac{1}{\gamma - 1} \log \frac{p}{\rho^{\gamma}}. \qquad (2) $$

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

$$ \rho v = \left(\frac{r_{\text{sh}}}{r}\right)^g \rho_{\text{sh}} v_{\text{sh}}, \qquad (3) $$

$$ \frac{\partial S}{\partial r} = \frac{\mathcal{L}}{\rho v}, \qquad (4) $$

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

where $c \equiv (\gamma p/\rho)^{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 $\epsilon$, which is a measure of the energy loss $\Delta e_{\text{disso}}$ per unit of mass in units of the pre-shock kinetic energy density $v_1^2/2$ as in Huete et al. (2018):

$$ \epsilon \equiv \frac{\Delta e_{\text{disso}}}{v_1^2/2}. \qquad (6) $$

The effect of $\epsilon$ on the post-shock Mach number $M_{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_1^2 \sim 2GM_{ns}/r_{sh}$. Defining the Mach number $M \equiv |v|/c$ as positive, we assume a strong adiabatic shock $M_1 \gg 1$ in the numerical calculations of this section.

The adiabatic compression of the post-shock flow by the gravitational potential $\Phi$ produces an inward increase in the enthalpy $c^2/(\gamma - 1)$ and thus an increase in the gas temperature $T \propto c^2$ according to the Bernoulli Eq. (5). In our model of stationary accretion, the temperature profile displays a maximum at a radius denoted $r_{\text{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 $(\alpha, \beta)$ of the cooling function $\mathcal{L}$, the mass accretion rate, and the geometry. With $\alpha = 3/2$ and $\beta = 5/2$, Walk et al. (2023) noted that both the radius $r_{\text{peak}} \sim 1.2r_{\text{ns}}$ of maximum temperature and the typical width $\Delta r_{\text{peak}} \sim 1.5r_{\text{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_{\text{peak}}$ varies by less than 1% between the cylindrical and spherical geometries, and varies by less than 0.5% with $r_{\text{sh}}$ for $r_{\text{sh}}/r_{\text{ns}} > 3$.

Figure 1 compares $r_{\text{peak}}$ to the radius of maximum deceleration denoted $r_\gamma$, with $r_\gamma < r_{\text{peak}}$ over a large range of parameters $(\alpha, \beta)$, except for when $\beta$ approaches unity. Noting that $r_\gamma = r_{\text{ns}}$ for $\alpha < \beta$ (Foglizzo et al. 2007), Fig. 1 indicates that $r_{\text{peak}}/r_\gamma \sim 1.2$ for $(\alpha, \beta) = (3/2, 5/2)$, and $r_{\text{peak}}/r_\gamma \sim 1.0$ for $(\alpha, \beta) = (6, 1)$. In Sect. 2.2, the parameters $(\alpha, \beta)$ 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 $l, m$ of spherical harmonics and a complex
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Ratios $r_{\text{peak}}/r_{\text{ns}}$ (thick solid lines) and $r_{\nabla}/r_{\text{ns}}$ (thin dotted lines) depending on the coefficients $(\alpha, \beta)$ that define the cooling function in Eq. (1). The contour lines are calculated for $\gamma = 4/3$ in spherical geometry, with $r_{\text{sh}}/r_{\text{ns}} = 10$, without dissociation ($\epsilon = 0$). The values of $r_{\text{peak}}/r_{\text{ns}}$ and $r_{\nabla}/r_{\text{ns}}$ are indicated with the same colour code. The cooling parameters $(\alpha, \beta) = (3/2, 5/2)$, (6, 1) and (5, 6) used in Fig. 3 are indicated with crosses for reference and are connected by grey lines. The black dotted line marks the threshold $\beta = \alpha$ above which the advection time from $r_{\text{peak}}$ to $r_{\text{ns}}$ is finite and $r_{\nabla} = r_{\text{ns}}$.</figcaption></figure>

eigenfrequency $\omega$. The eigenfrequency is independent of the azimuthal wavenumber $|m| \le l$ 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 $\delta f$, the perturbed mass flux $\delta h$, and the perturbed dimensionless entropy $\delta S$. The perturbation $\delta K$ is a combination of perturbed entropy $\delta S$ and the quantity $\delta w_{\perp}$ defined from the radial component of the curl of the vorticity perturbation $\delta w \equiv \nabla \times \delta v$:

$$ \delta f \equiv v\delta v_r + \frac{\delta c^2}{\gamma - 1}, \qquad (7) $$

$$ \delta h \equiv \frac{\delta v_r}{v} + \frac{\delta \rho}{\rho}, \qquad (8) $$

$$ \delta S \equiv \frac{1}{\gamma - 1} \frac{\delta c^2}{c^2} - \frac{\delta \rho}{\rho}, \qquad (9) $$

$$ \delta w_{\perp} \equiv r(\nabla \times \delta w)_r, \qquad (10) \\ = \frac{1}{\sin \theta} \left[ \frac{\partial}{\partial \theta} (\delta w_{\varphi} \sin \theta) - im\delta w_{\theta} \right], \qquad (11) $$

$$ \delta K \equiv rv\delta w_{\perp} + l(l+1)\frac{c^2}{\gamma}\delta S. \qquad (12) $$

We introduce the quantity $\delta A$, which is defined as the divergence of the ortho-radial velocity $\delta v_{\perp} \equiv (0, \delta v_{\theta}, \delta v_{\varphi})$:

$$ \delta A \equiv r^2 \nabla \cdot \delta v_{\perp}, \qquad (13) $$

$$ = \frac{r}{\sin \theta} \left[ \frac{\partial}{\partial \theta} (\delta v_{\theta} \sin \theta) + im\delta v_{\varphi} \right]. \qquad (14) $$

We establish in Appendix B that $\delta A$ and $\delta w_{\perp}$ are related to the perturbation $\delta v_r$ of the radial velocity as

$$ \delta v_r = \frac{r\delta w_{\perp}}{l(l+1)} - \frac{1}{l(l+1)} \frac{\partial \delta A}{\partial r}. \qquad (15) $$

This equation invites us to interpret $\delta A/\ell(l+1)$ as the potential defining the compressible part of the perturbed velocity, and

$r\delta w_{\perp}/\ell(l+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 $\delta A$ is related to $\delta K$ and $\delta f$ as

$$ \delta K = i\omega\delta A + \ell(l+1)\delta f. \qquad (16) $$

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

$$ \left( \frac{1 - M^2}{v} \frac{\partial}{\partial r} + \frac{i\omega}{c^2} \right) \frac{\delta A}{\ell(l+1)} = -\delta h - \left( \gamma - 1 + \frac{1}{M^2} \right) \frac{\delta S}{\gamma} + \frac{\delta K}{v^2 \ell(l+1)}, \qquad (17) $$

$$ \left( \frac{1 - M^2}{v} \frac{\partial}{\partial r} + \frac{i\omega}{c^2} \right) \delta h = \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \frac{\delta A}{\ell(l+1)} - \frac{i\omega}{v^2} \delta S + \frac{i\omega}{v^2 c^2} \frac{\delta K}{\ell(l+1)}, \qquad (18) $$

$$ \left( \frac{\partial}{\partial r} - \frac{i\omega}{v} \right) \delta S = \delta \left( \frac{\mathcal{L}}{\rho v} \right), \qquad (19) $$

$$ \left( \frac{\partial}{\partial r} - \frac{i\omega}{v} \right) \frac{\delta K}{\ell(l+1)} = \delta \left( \frac{\mathcal{L}}{\rho v} \right). \qquad (20) $$

This system includes explicit non-adiabatic terms only in the Eqs. (19, 20) governing $\delta S$ and $\delta K$. In Eq. (18), the Lamb frequency $\omega_{\text{Lamb}}$ associated with the spherical harmonic of order $l$ defines the turning point of non-radial acoustic waves:

$$ \omega_{\text{Lamb}}^2 = \ell(l+1) \frac{c^2}{r^2} (1 - M^2). \qquad (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):

$$ \frac{\delta A_{\text{sh}}}{\ell(l+1)} = - \left(1 - \frac{v_{\text{sh}}}{v_1}\right) v_1 \Delta \zeta, \qquad (22) $$

$$ \delta h_{\text{sh}} = \left(1 - \frac{v_{\text{sh}}}{v_1}\right) \frac{\Delta v}{v_{\text{sh}}}, \qquad (23) $$

$$ \delta S_{\text{sh}} = \gamma \frac{v_1}{c_{\text{sh}}^2} (i\omega + \omega_{\Phi}) \Delta \zeta \left(1 - \frac{v_{\text{sh}}}{v_1}\right)^2, \qquad (24) $$

$$ \delta K_{\text{sh}} = -\ell(l+1)\Delta\zeta \frac{c_{\text{sh}}^2}{\gamma} [\nabla S]_1^{\text{sh}}, \qquad (25) $$

where $\Delta\zeta$ and $\Delta v \equiv -i\omega\Delta\zeta$ are the shock displacement and velocity, respectively. The reference frequency $\omega_{\Phi}$ in Eq. (24) is defined as

$$ \frac{\omega_{\Phi} r_{\text{sh}}}{|v_{\text{sh}}|} = \frac{1}{2} \frac{v_1}{v_{\text{sh}}} \frac{\frac{2r_{\text{sh}}}{v_1^2} \frac{d\Phi}{dr} - 4 \frac{v_{\text{sh}}}{v_1}}{1 - \frac{v_{\text{sh}}}{v_1}} + \frac{\frac{v_{\text{sh}}}{v_1}}{\gamma M_{\text{sh}}^2} \frac{r_{\text{sh}} [\nabla S]_1^{\text{sh}}}{(1 - \frac{v_{\text{sh}}}{v_1})^2}, \qquad (26) $$

where $\omega_{\Phi}$ defines the threshold frequency separating the regime ($\omega \ll \omega_{\Phi}$) – where the phase of $\delta S_{\text{sh}}$ is opposed to that of $\Delta\zeta$ – from the regime ($\omega \gg \omega_{\Phi}$), where the phase of $\delta S_{\text{sh}}$ is set by that of $\Delta v$.

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

$$ \left[ \left( \frac{1 - M^2}{v} \frac{\partial}{\partial r} + \frac{i\omega}{c^2} \right)^2 + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right] \delta A = \frac{\partial}{\partial r} \frac{r\delta w_{\perp}}{v} - \ell(l+1) \frac{\gamma - 1}{\gamma} \delta \left( \frac{\mathcal{L}}{\rho v} \right). \qquad (27) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Growth rate (solid lines) and oscillation frequency (dashed lines) of the fundamental SASI mode, normalised by $|v_{sh}|/r_{sh}$, calculated with $\epsilon = 0$ for the cooling parameters (3/2, 5/2) (purple) and (6, 1) (orange), and displayed as functions of $r_{sh}/r_{\nabla}$. The analytic fitting formula (31) for $\omega_r$ is displayed with a green dashed line.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Growth rate (solid lines) and oscillation frequency (dashed lines) of the fundamental SASI mode, normalised by $|v_{sh}|/r_{sh}$, calculated with $r_{sh}/r_{\nabla} = 2$ (purple lines) and $r_{sh}/r_{\nabla} = 4$ (orange lines) for a range of cooling parameters $(\alpha, \beta)$ varying continuously from (3/2, 5/2) to (6, 1) (thin lines) and from (3/2, 5/2) to (5, 6) (thick lines), as indicated in Fig. 1.</figcaption></figure>

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 $\delta v_r$, $\delta \rho$, $\delta c$, and $\delta \rho$ are deduced from $\delta A$, $\delta h$, $\delta S$, and $\delta K$ using Eqs. (B.4)–(B.7), as detailed in Appendix B. In particular, using Eq. (B.4) to express the lower boundary condition $\delta v_r(r_{ns}) = 0$, leads to

$$ (c^2\delta h + c^2\delta S - \delta f)_{ns} = 0. \qquad (28) $$

With a vanishing sound speed $c_{ns} = 0$ resulting from the hypothesis of stationary flow, the lower boundary condition $\delta v_r(r_{ns}) = 0$ is equivalent to $\delta 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 $\omega$. The mathematical singularity at $r_{ns}$ is overcome numerically by using $\log M$ 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 ($M \sim 10^{-9}$).

Foglizzo et al. (2007) noted similar properties of the SASI harmonics as functions of $r_{sh}/r_{\nabla}$ with $(\alpha, \beta) = (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 $(\alpha, \beta)$ for $r_{sh}/r_{\nabla} \gg 1$, that is, when the cooling layer is thin compared to the shock radius. This surprising property is checked by varying $(\alpha, \beta)$ 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_{\nabla}$: the following analytic formula $\omega_r^{\text{fit}}$ provides an approximation for $\omega_r$ within 3% for $1.8 < r_{sh}/r_{\nabla} < 10$:

$$ \tau_{\text{adv}}^{\text{adiab}} = \frac{r_{\text{sh}}}{|v_{\text{sh}}|} \log \frac{r_{\text{sh}}}{r_{\nabla}}, \qquad (29) $$

$$ R_1 = \log \left( \frac{r_{\text{sh}}}{r_{\nabla}} - 1 \right), \qquad (30) $$

$$ \omega_r^{\text{fit}} = \frac{2\pi}{\tau_{\text{adv}}^{\text{adiab}}} \left( 0.56773 + 0.28628R_1 - 0.031763R_1^2 \right). \qquad (31) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Variation in the ratio $T_{SASI}/\tau_{\text{adv}}^{\text{adiab}}$ associated with the fundamental SASI mode, $\ell = 1$, calculated for the cooling parameters (3/2, 5/2) (purple) and (6, 1) (orange) and displayed as functions of $r_{sh}/r_{\nabla}$. The vertical lines highlight the range of $r_{sh}/r_{ns}$ in the analysis of model s25 in Müller & Janka (2014). The analytic fitting formula (31) is displayed with a dashed green line.</figcaption></figure>

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

The important role of $r_{\nabla}$ 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_{\nabla}$.

We note from Fig. 4 that the oscillation period $T_{SASI} \equiv 2\pi/\omega_r$ of the fundamental $\ell = 1$ SASI mode can be significantly longer than the approximate advection time $\tau_{\text{adv}}^{\text{adiab}}$ for a small ratio $r_{sh}/r_{\nabla}$. 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
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Estimate of the relative error between the oscillation period of the neutrino signal in model s25 in Müller & Janka (2014) and in analytical estimates. The empirical formula (32) is shown with a purple line. The result of the perturbative calculation is shown for $r_ν = r_{ns}$ without dissociation (orange line, Eq. (33)), and for $r_ν = 1.3r_{ns}$ with dissociation prescribed by Eq. (35) (green line, Eq. (36)). The sensitivity to the estimate of $r_ν$ is shown with dashed lines for $r_ν = 1.25r_{ns}$ (khaki) and $r_ν = 1.35r_{ns}$ (blue).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Effect of dissociation, measured by $\epsilon$, on the oscillation frequency of SASI for $r_{sh}/r_ν = 1.8$ (dashed lines), 1 (solid lines), and 2.3 (dotted lines) for $(\alpha, \beta) = (3/2, 5/2)$ (purple lines) and $(\alpha, \beta) = (6, 1)$ (orange lines). The analytical formula (34) is indicated with crosses.</figcaption></figure>

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_\odot$ 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:

$$T_{SASI}^{\text{MJ}} = 19 \text{ ms} \left( \frac{r_{\text{sh}}}{100 \text{ km}} \right)^{\frac{3}{2}} \log \left( \frac{r_{\text{sh}}}{r_{\text{ns}}} \right). \quad (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_{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:

$$T_{SASI} = \frac{T_{SASI}^{\text{MJ}} \times 1.82 \left( \frac{1.7 M_{\odot}}{M_{\text{ns}}} \right)^{\frac{1}{2}}}{0.56773 + 0.28628 R_1 - 0.031763 R_1^2}. \quad (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_\odot \sim 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

$$\omega_r(\epsilon) = \omega_r(0) - \left( \frac{2GM_{\text{ns}}}{r_{\text{sh}}^3} \right)^{\frac{1}{2}} \epsilon (0.19817 + 0.74804\epsilon). \quad (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):

$$\epsilon \sim 0.5 \left( \frac{r_{\text{sh}}}{150 \text{ km}} \right) \left( \frac{1.3 M_{\odot}}{M_{\text{ns}}} \right)^{\frac{1}{2}}, \quad (35)$$

with a saturation at $\epsilon \sim 0.5$ due to partial dissociation of $\alpha$-nuclei (Fernández & Thompson 2009a). Thus, Eq. (34) is rewritten as

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

with

$$\frac{T_{\text{SASI}}^{0}}{10.42 \text{ ms}} = \frac{\left(\frac{r_{\text{sh}}}{100 \text{ km}}\right)^{\frac{3}{2}} \left(\frac{1.7 M_{\odot}}{M_{\text{ns}}}\right)^{\frac{1}{2}} \log\left(\frac{r_{\text{sh}}}{r_{\nu}}\right)}{0.56773 + 0.28628 R_1 - 0.031763 R_1^2}. \quad (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 $\Delta\zeta$ is taken into account. Nuclear dissociation across the shock is formally taken into account using the parameter $\epsilon$ (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 $\mathcal{L} = 0$, the differential Equations (19)–(20) describe the advection of perturbations produced by the shock. Using the boundary conditions (24), (25) with $\nabla S = 0$:

$$ \delta S = \delta S_{\text{sh}} e^{i\omega \int_{\text{sh}} \frac{dr}{v}}, \quad (38) $$

$$ \delta K = 0. \quad (39) $$

Here, $\delta S_{\text{sh}}$ is defined by Eq. (24) with a simpler expression for $\omega_{\Phi}$:

$$ \frac{\omega_{\Phi} r_{\text{sh}}}{|v_{\text{sh}}|} \equiv \frac{\frac{1}{2} \frac{v_1}{v_{\text{sh}}} - g}{1 - \frac{v_{\text{sh}}}{v_1}}. \quad (40) $$

Using Eq. (12) with $\delta K = 0$ and Eq. (38) gives the explicit expression for $\delta w_{\perp}$ of

$$ \delta w_{\perp} = -\ell(\ell + 1) \frac{c^2}{\gamma r v} \delta S. \quad (41) $$

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

$$ \delta w_r = 0, \quad (42) $$

$$ \delta w_{\theta} = - \frac{imc^2}{rv \sin \theta} \frac{\delta S}{\gamma}, \quad (43) $$

$$ \delta w_{\varphi} = \frac{c^2}{rv} \frac{\partial}{\partial \theta} \frac{\delta S}{\gamma}. \quad (44) $$

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

Equation (16) implies that the wave action $\delta f/\omega$ is directly related to $\delta A$:

$$ \frac{\delta f}{i\omega} = - \frac{\delta A}{\ell(\ell + 1)}. \quad (45) $$

The transverse velocity components ($\delta v_{\theta}$, $\delta v_{\varphi}$) are deduced from $\delta A$ using the transverse Euler equations Eqs. (B.10), (B.11) with Eqs. (43), (44). Here, $\delta v_r$ is deduced from Eqs. (15) and (41):

$$ \delta v_r = - \frac{c^2}{\gamma v} \delta S - \frac{1}{\ell(\ell + 1)} \frac{\partial \delta A}{\partial r}, \quad (46) $$

$$ r \delta v_{\theta} = - \frac{1}{\ell(\ell + 1)} \frac{\partial \delta A}{\partial \theta}, \quad (47) $$

$$ r \delta v_{\varphi} = - \frac{im}{\ell(\ell + 1)} \frac{\delta A}{\sin \theta}. \quad (48) $$

The adiabatic solution can be viewed as the asymptotic limit of the model with a cooling function when $r_{\text{peak}} \sim r_{\text{ns}}$, corresponding to $\alpha \ll 1$ and $\beta \gg 1$ according to Fig. 1. For the sake of simplicity, we explore the adiabatic solution with an inner boundary defined by $\delta v_r = 0$, formulated by Eq. (28) with $c_{\text{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 $\delta 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 $\le 1.3$ for a large shock radius of $r_{\text{sh}} \sim 10r_{\text{ns}}$. The growth rate of the most unstable mode is overestimated by a factor of $\le 1.4$ for a small shock radius and is underestimated by a factor of $\le 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):

$$ dX \equiv \frac{v}{1 - M^2} dr. \quad (49) $$

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

$$ \left[ \left( \frac{\partial}{\partial X} + \frac{i\omega}{c^2} \right)^2 + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right] \delta A = \frac{\partial}{\partial X} \frac{r \delta w_{\perp}}{v}. \quad (50) $$

Using Eq. (41), the forcing term of Eq. (50) is thus proportional to the entropy perturbation $\delta S_{\text{sh}}$ produced by the shock:

$$ \left[ \left( \frac{\partial}{\partial X} + \frac{i\omega}{c^2} \right)^2 + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right] \frac{\delta A}{\ell(\ell + 1)} = - \mathcal{F}_S \delta S_{\text{sh}}, \quad (51) $$

$$ \mathcal{F}_S \equiv \frac{\partial}{\partial X} \frac{e^{i\omega \int_{\text{sh}} \frac{dr}{v}}}{\gamma M^2}. \quad (52) $$

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

$$ \left[ \left( \frac{\partial}{\partial X} + \frac{i\omega}{c^2} \right)^2 + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right] \delta \mathbf{L} = \frac{\partial}{\partial X} \frac{r \delta \mathbf{w}}{v}, \quad (53) $$

$$ = \frac{\partial}{\partial X} \frac{r v \delta \mathbf{w}}{c^2 M^2}, \quad (54) $$
where $\delta L \equiv r \times \delta v$ is the perturbed specific angular momentum vector and $\delta w \equiv \nabla \times \delta v$ is the perturbed vorticity vector. We note from Eqs. (43), (44) that $rv|\delta w|/c^2$ is conserved when advected. The relation between $\delta S_{\text{sh}}$ and $\delta A_{\text{sh}}$ is deduced from Eqs. (22) and (24):

$$ \ell(\ell + 1) \frac{c_{\text{sh}}^2 \delta S_{\text{sh}}}{\gamma \delta A_{\text{sh}}} = - \left(1 - \frac{v_{\text{sh}}}{v_1}\right) (i\omega + \omega_{\Phi}). \quad (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):

$$ \left( \frac{rv\delta w_{\theta}}{\delta L_{\theta}} \right)_{\text{sh}} = \left( \frac{rv\delta w_{\varphi}}{\delta L_{\varphi}} \right)_{\text{sh}} = -\ell(\ell + 1) \frac{c_{\text{sh}}^2 \delta S_{\text{sh}}}{\gamma \delta A_{\text{sh}}}. \quad (56) $$

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

$$ \frac{\partial \delta A}{\partial X} \bigg|_{\text{sh}} = \frac{i\omega}{v_{\text{sh}}^2} \delta A_{\text{sh}} \left[ 1 - 2 \frac{v_{\text{sh}}}{v_1} + (\gamma - 1) M_{\text{sh}}^2 \left( 1 - \frac{v_{\text{sh}}}{v_1} \right) \right] \\ - \left[ 1 + (\gamma - 1) M_{\text{sh}}^2 \right] \frac{\delta A_{\text{sh}}}{r_{\text{sh}} v_{\text{sh}}} \left( \frac{v_1}{2v_{\text{sh}}} - 2 \right), \quad (57) $$

$$ \frac{\partial \delta A}{\partial X} \bigg|_{\text{ns}} = \frac{i\omega}{c_{\text{ns}}^2} \delta A_{\text{ns}} - \ell(\ell + 1) \frac{1 - M_{\text{ns}}^2}{M_{\text{ns}}^2} \frac{\delta S_{\text{sh}}}{\gamma} e^{i\omega \int_{\text{sh}}^{\text{ns}} \frac{dX}{v}}. \quad (58) $$

They can also be written with the variables $r\delta v_{\theta}$, $r\delta v_{\varphi}$, and $\delta w_{\theta}$, $\delta w_{\varphi}$, using Eqs. (43), (44) and (47), (48). Equations (51) and (53) are thus equivalent for $\ell \ge 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 $\delta Y$ as follows:

$$ \delta Y \equiv \frac{\delta f}{i\omega} e^{i\omega \int_{\text{sh}}^{\text{sh}} \frac{dX}{v^2}} = - \frac{\delta A}{\ell(\ell + 1)} e^{i\omega \int_{\text{sh}}^{\text{sh}} \frac{dX}{v^2}}. \quad (59) $$

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

$$ \left\{ \frac{\partial^2}{\partial X^2} + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right\} Y_0 = 0, \quad (60) $$

$$ \frac{\partial Y_0}{\partial X} \bigg|_{\text{ns}} = \frac{i\omega}{c_{\text{ns}}^2} Y_0(r_{\text{ns}}). \quad (61) $$

We refer to $Y_0$ 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$:

$$ \left\{ \frac{\partial^2}{\partial r^2} + \left( \frac{\partial \log \left( 1 - M^2 \right)}{\partial r} \right) \frac{\partial}{\partial r} + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{c^2 (1 - M^2)^2} \right\} Y_0 = 0. \quad (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 ($\propto 1/M^2$ 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 $M$ 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 ($\propto 2\pi|v|/\omega_r$) of advected perturbations. A trade-off between effects (i) and (ii) favours advected perturbations with

<figure><img src="image_8.png" /><figcaption>Fig. 7. Oscillation frequency (upper plot) and growth rate (lower plot) of the modes $\ell = 1$ calculated in units of the post-shock frequency $v_{\text{sh}}/r_{\text{sh}}$, for $\gamma = 4/3$, as a function of the shock distance in the model with cooling using $(\alpha, \beta) = (3/2, 5/2)$ (dashed lines) and in the adiabatic approximation (solid lines). The fundamental mode (in red) and the first three overtones (green, orange, blue) are displayed. The grey horizontal line in the upper plot indicates the Lamb frequency at the shock. The fundamental mode becomes dominated by higher overtones as its frequency becomes too low for acoustic propagation ($\omega_r < \omega_{\text{Lamb}}^{\text{sh}}$).</figcaption></figure>

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$:

$$a'_1 Y_0^{\text{sh}} + a'_2 r_{\text{sh}} \left(\frac{\partial Y_0}{\partial r}\right)_{\text{sh}} = -M_{\text{sh}}^2 e^{\int_{Y_0}^{\text{sh}} \frac{i\omega}{\nu} \frac{dr}{1-M^2}} Y_0^{\text{ns}} \\ - \int_{\text{ns}}^{\text{sh}} \frac{\partial}{\partial r} \left(Y_0 e^{\int_{\text{sh}}^{\text{sh}} \frac{i\omega M^2}{1-M^2} \frac{dr}{\nu}}\right) \frac{M_{\text{sh}}^2}{M^2} e^{\int_{\text{sh}}^{\text{sh}} \frac{i\omega}{\nu} dr} dr, \qquad (63)$$

with $a'_1, a'_2$ defined as

$$a'_1 = (\gamma - 1)M_{\text{sh}}^2 + \frac{i\omega}{i\omega + \omega_{\Phi}} \frac{\frac{\nu_{\text{sh}}}{v_1}}{1 - \frac{\nu_{\text{sh}}}{v_1}}, \qquad (64)$$

$$a'_2 = - \frac{1 - M_{\text{sh}}^2}{\left(1 - \frac{\nu_{\text{sh}}}{v_1}\right) \frac{(i\omega + \omega_{\Phi})r_{\text{sh}}}{v_{\text{sh}}}}. \qquad (65)$$

The first two terms associated with $Y_0$ 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_0$. 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 $\delta F$ and the matching of its phase compared to the phase of the oscillator $Y_0$. 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 $\mathcal{L} = 0$, with the boundary conditions defined by Eq. (28) and Eqs. (22), (25) with $\nabla 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 $M_1 \gg 1$ and a large shock radius $r_{\text{sh}} \gg r_{\text{ns}}$, with the post-shock energy density described by the Bernoulli equation, Eq. (5), dominated by the enthalpy and the gravitational contributions:

$$\frac{c^2}{\gamma - 1} \sim \frac{GM_{\text{ns}}}{r}. \qquad (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 $\varepsilon$ impacts the post-shock Mach number and velocity according to Eqs. (A.2), (A.3):

$$M_{\text{sh}}^2 \sim \frac{\gamma - 1}{2\gamma} (1 - \varepsilon), \qquad (67)$$

$$\frac{\nu_{\text{sh}}}{v_1} \sim \frac{(\gamma - 1)(1 - \varepsilon)}{2 + (\gamma - 1)(1 - \varepsilon)} \le \frac{\gamma - 1}{\gamma + 1}, \qquad (68)$$

$$\frac{c_{\text{sh}}^2 r_{\text{sh}}}{GM_{\text{ns}}} \sim \frac{4\gamma(\gamma - 1)(1 - \varepsilon)}{[2 + (\gamma - 1)(1 - \varepsilon)]^2}. \qquad (69)$$

Equation (66) allows a power-law approximation of the radial velocity profile with an exponent $\alpha_v$ deduced from Eq. (A.11).

This exponent depends on both the adiabatic index and the geometry $g = 2$:

$$\alpha_v \equiv \frac{1}{\gamma - 1} - g, \qquad (70)$$

$$v = v_{\text{sh}} \left( \frac{r}{r_{\text{sh}}} \right)^{\alpha_v}, \qquad (71)$$

$$\frac{M_{\text{sh}}}{M} \sim \left( \frac{r_{\text{sh}}}{r} \right)^{\alpha_v + \frac{1}{2}}. \qquad (72)$$

We note that Eq. (66) is valid only for $r_{\text{ns}} \le r \ll r_{\text{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 $\ell \ge 1$, we approximate in Appendix E the homogeneous solution $Y_0$ as a linear combination of power laws independent of the frequency for $\omega \ll \omega_{\text{Lamb}}$:

$$Y_0 \sim \left(\frac{r}{r_{\text{ns}}}\right)^a \left[\left(\frac{r_{\text{ns}}}{r}\right)^b + \frac{b-a}{b+a}\left(\frac{r}{r_{\text{ns}}}\right)^b\right], \qquad (73)$$

$$\frac{\partial Y_0}{\partial r} \sim \frac{b-a}{r_{\text{ns}}} \left(\frac{r}{r_{\text{ns}}}\right)^{a-1} \left[\left(\frac{r}{r_{\text{ns}}}\right)^b - \left(\frac{r_{\text{ns}}}{r}\right)^b\right], \qquad (74)$$

where the exponents $a, b$ are defined as

$$a \equiv \frac{\alpha_v + 1}{2}, \qquad (75)$$

$$b \equiv \left[ \left( \frac{\alpha_v + 1}{2} \right)^2 + \ell(\ell + 1) \right]^{\frac{1}{2}}. \qquad (76)$$

The exponents $\pm b$ correspond to the two components of the acoustic structure $Y_0$, 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_0$ satisfies the lower boundary condition only in the asymptotic limit $r_{\text{ns}} \ll r_{\text{sh}}$.

### 5.1.3. Forcing by advected perturbations

We use the following approximation of the advection time $\tau_{\text{adv}}(r)$ profile:

$$\tau_{\text{adv}}(r) \equiv \int_{\text{sh}}^{r} \frac{dr}{v}, \qquad (77)$$

$$\sim \frac{1}{\alpha_v - 1} \frac{r_{\text{sh}}}{|v_{\text{sh}}|} \left[ \left( \frac{r_{\text{sh}}}{r} \right)^{\alpha_v - 1} - 1 \right] \quad \text{if } \alpha_v \neq 1, \qquad (78)$$

$$\sim \frac{r_{\text{sh}}}{|v_{\text{sh}}|} \log \frac{r_{\text{sh}}}{r} \quad \text{if } \alpha_v = 1. \qquad (79)$$

In particular, Eq. (79) is applicable for $\gamma = 4/3$. In this regime, the advection time $\tau_{\text{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_{\text{sh}}/r_{\text{ns}} < 100$.

We examine the accuracy of our approximation of $M$ and $Y_0$ in Fig. 8.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Solution of the homogeneous equation $Y_0$ associated with the fundamental mode $\ell = 1$ (solid orange line) for $r_{sh}/r_{ns} = 5$. This solution is very well approximated analytically by Eq. (73) (dashed orange line with crosses). The radial profile of the amplitude of the forcing term in Eq. (63) (purple solid line) is compared to its analytical approximation (dashed purple line with crosses) using Eqs. (72) and (74). The Mach number profile in the adiabatic model (green solid line) is compared to its analytical approximation (Eq. (72), green dashed line with crosses). The Mach profile of the flow including cooling is shown for reference (dark green line).</figcaption></figure>

## 5.2. Analytical expression of the fundamental eigenfrequency

We approximate Eq. (63) using Eq. (72) and Eq. (73) with $M^2 \ll 1$, except near the shock. We use the approximations (73) and (74) of $Y_0$ and $\partial Y_0/\partial r$ and note that $(\partial \log Y_0/\partial \log r)_{sh} \sim 1/(b+a)$ with $a, b$ defined by Eqs. (75) and (76). We define $N \equiv a'_2 + a'_1/(b+a)$ and introduce the normalised eigenfrequency $Z \equiv i\omega r_{sh}/|v_{sh}|$:

$$N(Z) \equiv \frac{\gamma - 1}{b + a} M_{sh}^2 + \frac{1 - M_{sh}^2 - \frac{Z}{b+a} \frac{v_{sh}}{v_1}}{\left(1 - \frac{v_{sh}}{v_1}\right) \left(Z + \frac{\omega_{\Phi} r_{sh}}{|v_{sh}|}\right)}. \quad (80)$$

Using the variable $x \equiv r/r_{ns}$, for $\ell \ge 1$ the eigenfrequency equation, Eq. (63), is approximated by

$$-x_{\text{sh}}^{3a-b-1} \int_{1}^{x_{\text{sh}}} (x^b - x^{-b}) e^{i\omega \int_{x_{\text{sh}}}^{x} \frac{dr}{r}} \frac{dx}{x^{3a}} \\ = N \left( \frac{i\omega r_{\text{sh}}}{|v_{\text{sh}}|} \right) + \frac{2b}{\ell(\ell+1)} \frac{M_{\text{sh}}^2}{x_{\text{sh}}^{a+b}} e^{i\omega r_{\text{adv}}^{\text{ns}}}. \quad (81)$$

The approximation of the eigenfrequency with $\gamma = 4/3$ allows an explicit calculation of the integral involved in Eq. (81). The following equation defining the eigenfrequencies for $\ell \ge 1$ is expressed using the complex amplification factor $Q(\omega)$ per advective-acoustic cycle:

$$Q(Z) \equiv \frac{2b \left(\frac{r_{\text{sh}}}{r_{\text{ns}}}\right)^{2-b} \left\{1 + \left[(Z+2)^2 - b^2\right] \frac{M_{\text{sh}}^2}{\ell(\ell+1)x_{\text{sh}}^3}\right\}}{[1 - (Z+2-b)N](Z+2+b) - \frac{Z+2-b}{x_{\text{sh}}^{2b}}}, \quad (82)$$

$$Q\left(\frac{i\omega r_{\text{sh}}}{|v_{\text{sh}}|}\right) e^{i\omega r_{\text{adv}}^{\text{ns}}} = 1, \quad (83)$$

where $b^2 = 1+\ell(\ell+1)$ and the advection time $\tau_{\text{adv}}^{\text{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.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Eigenfrequency of the adiabatic solution (green lines) based on the integral formulation (63) compared to its integrated approximation (83) (orange lines), which is based on the approximate description of the flow profile $v$, $M$ and acoustic structure $Y_0$. The estimate of the oscillation frequency $2\pi/\tau_{\text{adv}}^{\text{ns}}$ neglecting the phase of $Q$ is shown as a dashed blue line. For reference, the eigenfrequency of the flow with cooling using $(\alpha, \beta) = (3/2, 5/2)$ is also shown (purple lines with crosses).</figcaption></figure>

The good agreement obtained with the analytic expression (83) demonstrates that the approximation of the flow profile $v$, $M$ and the approximation of the homogeneous solution $Y_0$ 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_0$, without significant consequences regarding the approximation of the fundamental frequency $\omega_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_{\text{adv}}^{\text{ns}}$, also shown in Fig. 9, is mainly due to the neglect of the frequency dependence of the phase of $Q$, denoted $\varphi_Q$. The complex Eq. (83) is equivalent to the following set of two real equations:

$$\omega_i = \frac{\log|Q|}{\tau_{\text{adv}}^{\text{ns}}}, \quad (84)$$

$$\omega_r = \frac{2\pi - \varphi_Q}{\tau_{\text{adv}}^{\text{ns}}}. \quad (85)$$

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

Comparing $\omega_r$ and $2\pi/\tau_{\text{adv}}^{\text{ns}}$ in Fig. 9, it is particularly striking that the oscillation period $T_{\text{SASI}} \equiv 2\pi/\omega_r$ is significantly longer than the adiabatic advection time $\tau_{\text{adv}}^{\text{ns}}$. This contrasts with the simulations in Scheck et al. (2008) where $T_{\text{SASI}} \sim \tau_{\text{adv}}^{\text{ns}}$, in line with the analytic toy model of F09. Not only does $T_{\text{SASI}}$ differ from $\tau_{\text{adv}}^{\text{ns}}$ by more than 76%, but $T_{\text{SASI}}$ is actually 76% longer than the longest advection timescale $\tau_{\text{adv}}^{\text{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_{\text{SASI}}$ and $\tau_{\text{adv}}^{\text{ns}}$ does depend on the phase $\varphi_Q$ of the complex efficiency $Q$.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Amplitude $|Q|$ (solid lines) and normalised phase $\varphi_Q/2\pi$ (dashed lines) of the complex amplification factor associated through Eqs. (86), (87) to the fundamental eigenfrequency $\omega$ calculated with cooling (purple lines with crosses) and in the adiabatic approximation (green lines), using the same adiabatic estimate of the advection time $\tau_{\text{adv}}^{\text{ns}}$ as that in Eq. (79). The orange line shows $|Q|$ and $\varphi_Q/2\pi$ deduced from the analytic estimate (orange lines) using Eq. (82) with $\omega$ satisfying Eq. (83).</figcaption></figure>

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

$$|Q| = e^{\omega_r \tau_{\text{adv}}^{\text{ns}}}, \qquad (86)$$

$$\varphi_Q = 2\pi - \omega_r \tau_{\text{adv}}^{\text{ns}}. \qquad (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 $\Phi_{\text{adv}} \equiv 2\pi - \varphi_Q$. We note in Fig. 10 that $\varphi_Q < \pi$, 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 $Q \sim 2$ of the advective-acoustic cycle is arbitrarily set by the specific value (0.75) chosen for the ratio $c_{\text{sh}}^2/c_{\text{out}}^2$, 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_{\nabla}/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_{\text{sh}}^2/c_{\text{out}}^2$ and width $H_{\nabla}/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 $\omega_{\Phi} = 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_{\text{sh}}$, the radius of maximum deceleration $r_{\nabla}$, and the central mass $M_{\text{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_{\text{sh}}/r_{\nabla}$.

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_{\nabla}/r_{\text{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 $\delta S$ and $\delta 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.

*Acknowledgements.* 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

Blondin, J. M., & Mezzacappa, A. 2006, ApJ, 642, 401

Blondin, J. M., & Mezzacappa, A. 2007, Nature, 445, 58

Blondin, J. M., & Shaw, S. 2007, ApJ, 656, 366

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

Blondin, J. M., Gipson, E., Harris, S., & Mezzacappa, A. 2017, ApJ, 835, 170

Buellet, A. C., Foglizzo, T., Guilet, J., & Abdikamalov, E. 2023, A&A, 674, A205

Burrows, A., & Vartanyan, D. 2021, Nature, 589, 29

Drago, M., Andresen, H., Di Palma, I., Tamborra, I., & Torres-Forné, A. 2023, Phys. Rev. D, 108, 103036

Dunham, S. J., Endeve, E., Mezzacappa, A., et al. 2024, ApJ, 964, 38

Fernández, R. 2010, ApJ, 725, 1563

Fernández, R., & Thompson, C. 2009a, ApJ, 703, 1464

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

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

Foglizzo, T. 2001, A&A, 368, 311

Foglizzo, T. 2002, A&A, 392, 353

Foglizzo, T. 2009, ApJ, 694, 820

Foglizzo, T., Galletti, P., & Ruffert, M. 2005, A&A, 435, 397

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

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

Foglizzo, T., Masset, F., Guilet, J., & Durand, G. 2012, Phys. Rev. Lett., 108, 051103

Foglizzo, T., Kazeroni, R., Guilet, J., et al. 2015, PASA, 32, e009

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

Huete, C., Abdikamalov, E., & Radice, D. 2018, MNRAS, 475, 3305

Janka, H.-T. 2017, ApJ, 837, 84

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

Kitaura, F. S., Janka, H. T., & Hillebrandt, W. 2006, A&A, 450, 345

Kotake, K., & Kuroda, T. 2017, in Handbook of Supernovae, eds. A. W. Alsabti, & P. Murdin, 1671

Müller, B. 2019a, MNRAS, 487, 5304

Müller, B. 2019b, Ann. Rev. Nucl. Part. Sci., 69, 253

Müller, B. 2020, Liv. Rev. Comput. Astrophys., 6, 3

Müller, B., & Janka, H.-T. 2014, ApJ, 788, 82

Müller, B., Tauris, T. M., Heger, A., et al. 2019, MNRAS, 484, 3307

Powell, J., & Müller, B. 2022, Phys. Rev. D, 105, 063018

Sato, J., Foglizzo, T., & Fromang, S. 2009, ApJ, 694, 833

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

Stockinger, G., Janka, H. T., Kresse, D., et al. 2020, MNRAS, 496, 2039

Tamborra, I., & Murase, K. 2019, Supernovae. Ser.: Space Sci. Ser. ISSI, 68, 87

Walk, L., Foglizzo, T., & Tamborra, I. 2023, Phys. Rev. D, 107, 063014
## Appendix A: Stationary accretion

The loss of energy through nuclear dissociation is measured by the parameter $\epsilon$ defined by Eq. (6). It decreases the value of the post-shock Mach number $M_{sh}$ below the reference adiabatic value $M_{ad}$, and affects the Rankine-Hugoniot conditions according to Eqs. (A4-A6) of Foglizzo et al. (2006):

$$M_{\text{ad}}^2 = \frac{2 + (\gamma - 1)M_1^2}{2\gamma M_1^2 - \gamma + 1}, \quad (\text{A.1})$$

$$\epsilon = \left(1 + \frac{2}{\gamma - 1} \frac{1}{M_1^2}\right) \left(1 - \frac{M_{\text{sh}}^2}{M_{\text{ad}}^2}\right) \left(1 - \frac{M_{\text{sh}}^2}{M_1^2}\right), \quad (\text{A.2})$$

$$\frac{v_{\text{sh}}}{v_1} = \frac{M_{\text{sh}}^2}{M_1^2} \frac{1 + \gamma M_1^2}{1 + \gamma M_{\text{sh}}^2}, \quad (\text{A.3})$$

$$\frac{c_{\text{sh}}}{c_1} = \frac{M_{\text{sh}}}{M_1} \frac{1 + \gamma M_1^2}{1 + \gamma M_{\text{sh}}^2}. \quad (\text{A.4})$$

The nuclear binding energy of the iron nuclei is 8.8MeV per nucleon, or $1.77 \times 10^{52}$ erg/$M_{\odot}$. If nuclear dissociation was complete across the shock this would translates into a dissociation parameter scaling linearly with $r_{\text{sh}}$:

$$\epsilon = 0.7 \left( \frac{r_{\text{sh}}}{150\text{km}} \right) \left( \frac{1.3 M_{\odot}}{M_{\text{ns}}} \right). \quad (\text{A.5})$$

We note that our definition of $\epsilon$ 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 $\epsilon$ formulated in Eq. (35) is a rough estimate deduced from Eq. (51) and Fig. 12 in Huete et al. (2018) for $r_{\text{sh}} < 175$km. It is smaller than Eq. (A.5) by a factor ~ 1.4 and saturates at $\epsilon \sim 0.5$ in exploding models.

The value of $M_{sh}$ is expressed as a function of $M_1$ and $\epsilon(r_{\text{sh}})$ using Eq. (A.2):

$$\frac{M_{\text{sh}}^2}{M_{\text{ad}}^2} = \frac{1}{2} + \frac{1 - \frac{2\epsilon}{1 + \frac{2}{\gamma-1} \frac{1}{M_1^2}} - \frac{M_{\text{ad}}^2}{2M_1^2}}{1 + \left[ \left(1 - \frac{M_{\text{ad}}^2}{M_1^2}\right)^2 + \frac{4\epsilon \frac{M_{\text{ad}}^2}{M_1^2}}{1 + \frac{2}{\gamma-1} \frac{1}{M_1^2}} \right]^{\frac{1}{2}}}, \quad (\text{A.6}) \\ \sim 1 - \epsilon \text{ if } M_1 \gg 1 \quad (\text{A.7})$$

The jump condition (A.3) is thus a function of $\epsilon(r_{\text{sh}})$ using Eq. (A.6), with the simple following expression for a strong shock:

$$\frac{v_1}{v_{\text{sh}}} \sim 1 + \frac{2}{(1 - \epsilon)(\gamma - 1)}, \quad (\text{A.8})$$

$$\sim 1 + \frac{6}{1 - \epsilon}. \quad (\text{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, M, v$ from the stationary equations, Eqs. (4) and (5):

$$\frac{\rho}{\rho_{\text{sh}}} = \left( \frac{c}{c_{\text{sh}}} \right)^{\frac{2}{\gamma-1}} e^{-(S-S_{\text{sh}})}, \quad (\text{A.10})$$

$$r^{\gamma} c^{\frac{\gamma+1}{\gamma-1}} M = r_{\text{sh}}^{\gamma} c_{\text{sh}}^{\frac{\gamma+1}{\gamma-1}} e^{S-S_{\text{sh}}}. \quad (\text{A.11})$$

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

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

$$dX \equiv \frac{v}{1 - M^2} dr. \quad (\text{B.1})$$

Thus

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

$$\left(\frac{\partial}{\partial X} + \frac{i\omega}{c^2}\right) \delta h = \frac{i\omega}{v^2} \left( \frac{\omega^2 - \omega_{\text{Lamb}}^2}{\omega^2 c^2} \delta f - \delta S \right) + \frac{1 - M^2}{v^2} \frac{i\delta K}{\omega r^2} \quad (\text{B.3})$$

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

$$\frac{\delta v_r}{v} = \frac{1}{1 - M^2} \left( \delta h + \delta S - \frac{\delta f}{c^2} \right), \quad (\text{B.4})$$

$$\frac{\delta \rho}{\rho} = \frac{1}{1 - M^2} \left( -M^2 \delta h - \delta S + \frac{\delta f}{c^2} \right), \quad (\text{B.5})$$

$$\frac{\delta c^2}{c^2} = \frac{\gamma - 1}{1 - M^2} \left( \frac{\delta f}{c^2} - M^2 \delta h - M^2 \delta S \right), \quad (\text{B.6})$$

$$\frac{\delta p}{\gamma p} = \frac{1}{\gamma - 1} \frac{\delta c^2}{c^2} - \frac{\delta S}{\gamma}, \quad (\text{B.7})$$

$$\delta \left(\frac{\mathcal{L}}{\rho v}\right) = \nabla S \frac{c^2}{\gamma} \left[ \left(\beta - 1\right) \frac{\delta \rho}{\rho} + \alpha \frac{\delta c^2}{c^2} - \frac{\delta v_r}{v} \right], \quad (\text{B.8})$$

$$\delta \left(\frac{\mathcal{L}}{\rho v}\right) = \nabla S \left[ \left(\beta - 1\right) \frac{\delta \rho}{\rho} + \left(\alpha - 1\right) \frac{\delta c^2}{c^2} - \frac{\delta v_r}{v} \right]. \quad (\text{B.9})$$

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

$$-i\omega\delta v_{\theta} + v\delta w_{\varphi} + \frac{1}{r}\frac{\partial}{\partial\theta}\delta f = \frac{c^2}{\gamma}\frac{1}{r}\frac{\partial\delta S}{\partial\theta}, \quad (\text{B.10})$$

$$-i\omega\delta v_{\varphi} - v\delta w_{\theta} + \frac{im}{r\sin\theta}\delta f = \frac{c^2}{\gamma}\frac{im\delta S}{r\sin\theta}, \quad (\text{B.11})$$

We eliminate $\delta v_{\theta}$ and $\delta v_{\varphi}$ in the system of Eqs. (12), (14), (B.10) and (B.11) and obtain Eq. (16).

The derivative of $\delta A$ is calculated using Eqs. (16), (20) and (B.2), leading to the differential system (17-20) satisfied by $\delta A, \delta h, \delta S, \delta K$. The second derivative of $\delta 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:

$$\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 M^2}, \quad (\text{B.12})$$

$$\frac{\delta \rho}{\rho} = \frac{1}{c^2} \left( v \frac{\partial}{\partial r} - i\omega \right) \frac{\delta A}{\ell(\ell + 1)} - \frac{\gamma - 1}{\gamma} \delta S, \quad (\text{B.13})$$

$$\frac{1}{\gamma - 1} \frac{\delta c^2}{c^2} = \frac{1}{c^2} \left( v \frac{\partial}{\partial r} - i\omega \right) \frac{\delta A}{\ell(\ell + 1)} + \frac{\delta S}{\gamma}, \quad (\text{B.14})$$

$$\frac{\delta p}{\gamma p} = \frac{1}{c^2} \left( v \frac{\partial}{\partial r} - i\omega \right) \frac{\delta A}{\ell(\ell + 1)}, \quad (\text{B.15})$$

Using Eq. (B.12) with Eq. (12), we can express $\delta v_r$ with $\delta A$ and $\delta w_{\perp}$ and obtain Eq. (15). We note that this relation is equivalent to the radial component of the vector calculus relation
$$\nabla^2 v = \nabla(\nabla \cdot v) - \nabla \times (\nabla \times 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:

$$\delta v_{\theta sh} = \frac{v_1 - v_{sh}}{r_{sh}} \frac{\partial \Delta \zeta}{\partial \theta}, \qquad (B.16)$$

$$\delta v_{\varphi sh} = \frac{v_1 - v_{sh}}{r_{sh}} \frac{im\Delta\zeta}{\sin\theta}. \qquad (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):

$$\delta w_{rsh} = 0, \qquad (B.18)$$

$$\delta w_{\theta sh} = -\frac{c_{sh}^2}{\gamma} \frac{im}{r_{sh}v_{sh}\sin\theta} (\delta S_{sh} + [\nabla S]_1^{sh} \Delta\zeta), \qquad (B.19)$$

$$\delta w_{\varphi sh} = \frac{c_{sh}^2}{\gamma} \frac{1}{r_{sh}v_{sh}} \frac{\partial}{\partial \theta} (\delta S_{sh} + [\nabla S]_1^{sh} \Delta\zeta). \qquad (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

$$\delta w_r = \left(\frac{r_{sh}}{r}\right)^2 \delta w_{rsh} e^{i\omega \int_{sh} \frac{dr}{v}}, \qquad (C.1)$$

$$\delta w_{\theta} = \frac{1}{rv} \left[ \left( rv\delta w_{\theta} \right)_{sh} - \frac{c^2 - c_{sh}^2}{\sin\theta} \frac{im}{\gamma} \delta S_{sh} \right] e^{i\omega \int_{sh} \frac{dr}{v}}, \qquad (C.2)$$

$$\delta w_{\varphi} = \frac{1}{rv} \left[ \left( rv\delta w_{\varphi} \right)_{sh} + \frac{c^2 - c_{sh}^2}{\sin\theta} \frac{\partial}{\partial\theta} \frac{\delta S_{sh}}{\gamma} \right] e^{i\omega \int_{sh} \frac{dr}{v}}. \qquad (C.3)$$

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

With $\delta Y$ defined by Eq. (59), the expressions of $\delta A_{sh}$, $\delta h_{sh}$ are deduced from Eqs. (22-24):

$$\left(1 - \frac{v_{sh}}{v_1}\right) \Delta\zeta = \frac{\delta Y_{sh}}{v_1}, \qquad (C.4)$$

$$\delta h_{sh} = -\frac{i\omega}{v_1 v_{sh}} \delta Y_{sh}, \qquad (C.5)$$

$$\frac{\delta S_{sh}}{\gamma} = \frac{\delta Y_{sh}}{c_{sh}^2} (i\omega + \omega_{\Phi}) \left(1 - \frac{v_{sh}}{v_1}\right). \qquad (C.6)$$

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

$$\delta Y_{sh} = -\frac{\delta A_{sh}}{\ell(\ell + 1)}, \qquad (C.7)$$

$$\left(\frac{\partial \delta A}{\partial r}\right)_{sh} = \frac{\ell(\ell + 1)v_{sh}}{1 - M_{sh}^2} \left[ \delta Y_{sh} \frac{i\omega}{v_{sh}^2} \left( \frac{v_{sh}}{v_1} + M_{sh}^2 \right) - \frac{\delta S_{sh}}{\gamma} \left( \frac{1}{M_{sh}^2} + \gamma - 1 \right) \right], \qquad (C.8)$$

$$\left(\frac{\partial \delta Y}{\partial X}\right)_{sh} = -\frac{1 - M_{sh}^2}{\ell(\ell + 1)v_{sh}} \left(\frac{\partial \delta A}{\partial r}\right)_{sh} + \frac{i\omega}{c_{sh}^2} \delta Y_{sh}, \qquad (C.9)$$

which results in Eq. (57).

At the inner boundary in the adiabatic approximation, we use the condition $\delta v_r = 0$:

$$\delta h_{ns} + \delta S_{ns} = \frac{\delta f}{c_{ns}^2}. \qquad (C.10)$$

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

$$\delta S_{ns} = \delta S_{sh} e^{\int_{sh}^{ns} i\omega \frac{dr}{v}}, \qquad (C.11)$$

$$\frac{\partial \delta Y}{\partial X} = \delta h e^{\int_{sh}^{sh} \frac{i\omega}{c^2} dX} + \frac{\delta S}{\gamma} e^{\int_{sh}^{sh} \frac{i\omega}{c^2} dX} \left( \frac{1}{M^2} + \gamma - 1 \right). \qquad (C.12)$$

The inner boundary condition (C.10) is thus reduced to Eq. (58). Multiplying Eq. (51) by $m/(\omega \sin \theta)$ and using Eqs. (43) and (48):

$$\left[ \left( \frac{\partial}{\partial X} + \frac{i\omega}{c^2} \right)^2 + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right] r \delta v_{\varphi} = - \frac{\partial}{\partial X} \frac{r \delta w_{\theta}}{v}. \qquad (C.13)$$

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

$$\left[ \left( \frac{\partial}{\partial X} + \frac{i\omega}{c^2} \right)^2 + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right] r \delta v_{\theta} = \frac{\partial}{\partial X} \frac{r \delta w_{\varphi}}{v}, \qquad (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 $\delta Y$ is deduced from Eqs. (51) and (59):

$$\left( \frac{\partial^2}{\partial X^2} + \frac{\omega^2 - \omega_{\text{Lamb}}^2}{v^2 c^2} \right) \delta Y = \delta \mathcal{F}, \qquad (D.1)$$

$$\delta \mathcal{F} \equiv e^{\int_{sh}^{sh} \frac{i\omega}{c^2} dX} \mathcal{F}_S \delta S_{sh}. \qquad (D.2)$$

We define $Y_0$ 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$:

$$\left(\frac{\partial Y_0}{\partial X}\right)_{ns} = \frac{i\omega}{c_{ns}^2} Y_0^{ns}, \qquad (D.3)$$

$$W \equiv Y_0 \frac{\partial Y_-}{\partial X} - Y_- \frac{\partial Y_0}{\partial X}, \qquad (D.4)$$

$$\delta Y = Y_{-} \left( d_{-} + \int_{ns}^{sh} \frac{Y_{0}}{W} \delta \mathcal{F} dX \right) - Y_{0} \left( d_{0} + \int_{sh}^{sh} \frac{Y_{-}}{W} \delta \mathcal{F} dX \right), \qquad (D.5)$$

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

$$\delta Y_{sh} = Y_{-}^{sh} \left( d_{-} + \int_{ns}^{sh} \frac{Y_{0}}{W} \delta \mathcal{F} dX \right) - d_{0} Y_{0}^{sh}, \qquad (D.6)$$

$$\left(\frac{\partial \delta Y}{\partial X}\right)_{sh} = \left(\frac{\partial Y_-}{\partial X}\right)_{sh} \left( d_- + \int_{ns}^{sh} \frac{Y_0}{W} \delta \mathcal{F} dX \right) - d_0 \left(\frac{\partial Y_0}{\partial X}\right)_{sh}, \qquad (D.7)$$
Using Eq. (D.5) at the lower boundary:

$$ \delta Y_{\text{ns}} = d_{-} Y_{-}^{\text{ns}} - Y_{0}^{\text{ns}} \left( d_{0} - \int_{\text{ns}}^{\text{sh}} \frac{Y_{-}}{W} \delta \mathcal{F} \, dX \right), \quad (\text{D.8}) $$

$$ \left( \frac{\partial \delta Y}{\partial X} \right)_{\text{ns}} = d_{-} \left( \frac{\partial Y_{-}}{\partial X} \right)_{\text{ns}} - \left( \frac{\partial Y_{0}}{\partial X} \right)_{\text{ns}} \left( d_{0} - \int_{\text{ns}}^{\text{sh}} \frac{Y_{-}}{W} \delta \mathcal{F} \, dX \right). \quad (\text{D.9}) $$

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

$$ d_{-} \left[ \left( \frac{\partial Y_{-}}{\partial X} \right)_{\text{ns}} - \frac{i\omega}{c_{\text{ns}}^2} Y_{-}^{\text{ns}} \right] = \delta \mathcal{F}_{\text{sh}} (1 - M_{\text{ns}}^2) \frac{M_{\text{sh}}^2}{M_{\text{ns}}^2} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dX} . \quad (\text{D.10}) $$

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

$$ \left( \frac{\partial Y_{-}}{\partial X} \right)_{\text{ns}} - \frac{i\omega}{c_{\text{ns}}^2} Y_{-}^{\text{ns}} = \frac{W}{Y_{0}^{\text{ns}}}. \quad (\text{D.11}) $$

Thus

$$ W d_{-} = Y_{0}^{\text{ns}} \delta \mathcal{F}_{\text{sh}} (1 - M_{\text{ns}}^{2}) \frac{M_{\text{sh}}^{2}}{M_{\text{ns}}^{2}} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dX}. \quad (\text{D.12}) $$

Eliminating $d_0$ between Eqs. (D.6) and (D.7) and using Eq. (D.4):

$$ Y_{0}^{\text{sh}} \left( \frac{\partial \delta Y}{\partial X} \right)_{\text{sh}} - \left( \frac{\partial Y_{0}}{\partial X} \right)_{\text{sh}} \delta Y_{\text{sh}} = \left[ Y_{0}^{\text{sh}} \left( \frac{\partial Y}{\partial X} \right)_{\text{sh}} - \left( \frac{\partial Y_{0}}{\partial X} \right)_{\text{sh}} Y_{-}^{\text{sh}} \right] \left( d_{-} + \int_{\text{ns}}^{\text{sh}} \frac{Y_{0}}{W} \delta \mathcal{F} \, dX \right), \quad (\text{D.13}) $$

$$ = W d_{-} + \int_{\text{ns}}^{\text{sh}} Y_{0} \delta \mathcal{F} \, dX. \quad (\text{D.14}) $$

The eigenfrequencies are thus defined by

$$ Y_{0}^{\text{sh}} \left( \frac{\partial \delta Y}{\partial X} \right)_{\text{sh}} - \left( \frac{\partial Y_{0}}{\partial X} \right)_{\text{sh}} \delta Y_{\text{sh}} = Y_{0}^{\text{ns}} \delta \mathcal{F}_{\text{sh}} (1 - M_{\text{ns}}^{2}) \frac{M_{\text{sh}}^{2}}{M_{\text{ns}}^{2}} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dX} + \int_{\text{ns}}^{\text{sh}} Y_{0} \delta \mathcal{F} \, dX. \quad (\text{D.15}) $$

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

$$ Y_{0}^{\text{sh}} \left( \frac{\partial \delta Y}{\partial X} \right)_{\text{sh}} - \left( \frac{\partial Y_{0}}{\partial X} \right)_{\text{sh}} \delta Y_{\text{sh}} = \delta \mathcal{F}_{\text{sh}} \left\{ Y_{0}^{\text{ns}} (1 - M_{\text{ns}}^{2}) \frac{M_{\text{sh}}^{2}}{M_{\text{ns}}^{2}} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dX} + \int_{\text{ns}}^{\text{sh}} Y_{0} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{c^2} dX} \frac{\partial}{\partial r} \left( \frac{M_{\text{sh}}^{2}}{M_{\text{ns}}^{2}} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dr} \right) dr \right\}. \quad (\text{D.16}) $$

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

$$ a_1 Y_0^{\text{sh}} + a_2 r_{\text{sh}} \left( \frac{\partial Y_0}{\partial r} \right)_{\text{sh}} + a_3 Y_0^{\text{ns}} = \int_{\text{ns}}^{\text{sh}} Y_0 e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{c^2} dX} \frac{\partial}{\partial r} \left( \frac{M_{\text{sh}}^2}{M_{\text{ns}}^2} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dr} \right) dr, \quad (\text{D.17}) $$

with $a_1, a_2, a_3$ defined by:

$$ a_1 \equiv \frac{\left( \frac{\partial \delta Y}{\partial X} \right)_{\text{sh}}}{\delta \mathcal{F}_{\text{sh}}}, \quad (\text{D.18}) $$

$$ = 1 + (\gamma - 1) M_{\text{sh}}^2 - \frac{i\omega}{i\omega + \omega_{\Phi}} \frac{\frac{v_{\text{sh}}}{v_1}}{1 - \frac{v_{\text{sh}}}{v_1}}, \quad (\text{D.19}) $$

$$ a_2 \equiv - \frac{1 - M_{\text{sh}}^2}{r_{\text{sh}} v_{\text{sh}}} \frac{\delta Y_{\text{sh}}}{\delta \mathcal{F}_{\text{sh}}}, \quad (\text{D.20}) $$

$$ = - \frac{1 - M_{\text{sh}}^2}{\left(1 - \frac{v_{\text{sh}}}{v_1}\right) \frac{(i\omega + \omega_{\Phi})r_{\text{sh}}}{v_{\text{sh}}}}, \quad (\text{D.21}) $$

$$ a_3 \equiv -(1 - M_{\text{ns}}^2) \frac{M_{\text{sh}}^2}{M_{\text{ns}}^2} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dX}. \quad (\text{D.22}) $$

After one integration by parts:

$$ Y_{0}^{\text{sh}} \left( \frac{\partial \delta Y}{\partial X} \right)_{\text{sh}} - \left( \frac{\partial Y_{0}}{\partial X} \right)_{\text{sh}} \delta Y_{\text{sh}} = \delta \mathcal{F}_{\text{sh}} \left\{ Y_{0}^{\text{sh}} - Y_{0}^{\text{ns}} M_{\text{sh}}^{2} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dX} - \int_{\text{ns}}^{\text{sh}} \frac{\partial}{\partial r} \left( Y_{0} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dX} \right) \frac{M_{\text{sh}}^{2}}{M_{\text{ns}}^{2}} e^{\int_{\text{sh}}^{\text{ns}} \frac{i\omega}{v} dr} dr \right\} \quad (\text{D.23}) $$

The equation defining the eigenfrequencies becomes Eq. (63) with $a'_1 \equiv a_1 - 1$ and $a'_2 \equiv a_2$ defined by Eqs. (64-65).

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

The dissociation measured by the parameter $\epsilon$ affects the relation between $|v_{\text{sh}}|$ and the local free fall velocity, as deduced from Eq. (68) for a strong shock:

$$ |v_{\text{sh}}| \left( \frac{r_{\text{sh}}}{2GM_{\text{ns}}} \right)^{\frac{1}{2}} \sim \frac{(\gamma - 1)(1 - \epsilon)}{2 + (\gamma - 1)(1 - \epsilon)}, \quad (\text{E.1}) $$

$$ \frac{c_{\text{sh}}^2 r_{\text{sh}}}{GM_{\text{ns}}} = \frac{1}{M_{\text{sh}}^2} \frac{v_{\text{sh}}^2 r_{\text{sh}}}{GM_{\text{ns}}}. \quad (\text{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 \ll r_{\text{sh}}$ and $M^2 \ll 1$:

$$ c^2 \left[ 1 + (\gamma - 1) \frac{M^2}{2} \right] = (\gamma - 1) \frac{GM_{\text{ns}}}{r} + c_{\text{sh}}^2 \left[ 1 + (\gamma - 1) \frac{M_{\text{sh}}^2}{2} \right] - (\gamma - 1) \frac{GM_{\text{ns}}}{r_{\text{sh}}}, \quad (\text{E.3}) $$

$$ c^2 \sim (\gamma - 1) \frac{GM_{\text{ns}}}{r}. \quad (\text{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 $\gamma = 4/3$ we approximate

$$ M^2 \sim M_{\text{sh}}^2 \left( \frac{r}{r_{\text{sh}}} \right)^3, \quad (\text{E.5}) $$

$$ v \sim v_{\text{sh}} \frac{r}{r_{\text{sh}}}, \quad (\text{E.6}) $$

$$ i\omega \int_{\text{sh}} \frac{M^2}{1-M^2} \frac{dr}{v} \sim -M_{\text{sh}}^2 \frac{i\omega r_{\text{sh}}}{|v_{\text{sh}}|} \int_{\text{sh}} \frac{x^2 dx}{1 - M_{\text{sh}}^2 x^3}, \quad (\text{E.7}) $$

$$ \sim \frac{i\omega r_{\text{sh}}}{3|v_{\text{sh}}|} \log \left( \frac{1 - M^2}{1 - M_{\text{sh}}^2} \right). \quad (\text{E.8}) $$
The oscillatory phase associated to $\delta f$ in Eq. (63) is made of a product of $(i\omega)$ with two explicit contributions and a third contribution from the definition of $\delta Y$ in Eq. (59). The sum of these three contributions is:

$$ \int_{\text{sh}} \frac{M^2}{1 - M^2} \frac{dr}{v} + \int_{\text{sh}} \frac{M^2}{1 - M^2} \frac{dr}{v} + \int_{\text{sh}} \frac{dr}{v} = \int_{\text{sh}} \frac{1 + M^2}{1 - M^2} \frac{dr}{v}. \quad (E.9) $$

With $\alpha_v = 1$ and $\tilde{r} \equiv r/r_{\text{sh}}$,

$$ \frac{1 + M^2 \frac{v_{\text{sh}}}{v}}{1 - M^2} \sim \frac{2 M_{\text{sh}}^2 \tilde{r}^2}{1 - M_{\text{sh}}^2 \tilde{r}^3} + \frac{1}{\tilde{r}}. \quad (E.10) $$

The integral can be estimated as follows:

$$ \int_{\text{sh}} \frac{1 + M^2}{1 - M^2} \frac{dr}{v} \sim \frac{r_{\text{sh}}}{|v_{\text{sh}}|} \left( \log \frac{r_{\text{sh}}}{r_{\text{ns}}} - \frac{2}{3} \log \frac{1 - M_{\text{sh}}^2}{1 - M_{\text{ns}}^2} \right), \quad (E.11) $$

$$ \sim \frac{r_{\text{sh}}}{|v_{\text{sh}}|} \left( \log \frac{r_{\text{sh}}}{r_{\text{ns}}} + \frac{2M_{\text{sh}}^2}{3} \right) \quad (E.12) $$

With a strong adiabatic shock $M_{\text{sh}}^2 \sim 1/8$, the correction to the advection timescale $\tau_{\text{adv}}^{\text{ns}}$ is thus of order 12% for $r_{\text{sh}}/r_{\text{ns}} = 2$ and 5% for $r_{\text{sh}}/r_{\text{ns}} = 5$. This correction is smaller if dissociation diminishes the value of $M_{\text{sh}}^2$.

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

$$ \int_{\text{sh}}^{r} \frac{\omega}{c^2} \frac{dX}{c^2} = \omega \int_{\text{sh}}^{r} \frac{M^2}{1 - M^2} \frac{dr}{v}, \quad (E.13) $$

$$ \sim \frac{\omega r_{\text{sh}}}{|v_{\text{sh}}|} \frac{M_{\text{sh}}^2}{\alpha_v + 2} \left[ 1 - \left( \frac{r}{r_{\text{sh}}} \right)^{\alpha_v + 2} \right]. \quad (E.14) $$

With $\omega_r \sim 2\pi|v_{\text{sh}}|/r_{\text{sh}}$ for a small shock radius this phase shift is not negligible since it reaches $0.74 \sim \pi/4$ at the inner boundary and it is linear in $\omega$.

We integrate Eq. (49) using Eq. (71) and neglecting $M^2 \ll 1$:

$$ X_{\text{sh}} \equiv \frac{r_{\text{sh}} v_{\text{sh}}}{\alpha_v + 1}, \quad (E.15) $$

$$ \frac{X}{X_{\text{sh}}} \sim \left( \frac{r}{r_{\text{sh}}} \right)^{\alpha_v + 1}. \quad (E.16) $$

The two solutions $Y_0^{\pm}$ of the homogeneous equation (60) are approximated as power laws with exponents $\alpha_{\pm}$.

$$ Y_0^{\pm}(r) \equiv \left( \frac{r}{r_{\text{ns}}} \right)^{\alpha_{\pm}}. \quad (E.17) $$

Injecting $Y_0^{\pm}(r)$ into Eq. (60),

$$ \alpha_{\pm} (\alpha_{\pm} - \alpha_v - 1) \left( \frac{r}{r_{\text{sh}}} \right)^{-2\alpha_v - 2} = - \left[ \frac{\omega^2 r_{\text{sh}}^2}{c^2} - \ell(\ell + 1) \frac{r_{\text{sh}}^2}{r^2} \right] \left( \frac{v_{\text{sh}}}{v} \right)^2. \quad (E.18) $$

For $\ell \ge 1$, the restoring force in Eq. (60) is independent of the frequency in the region where $\omega \ll \omega_{\text{Lamb}}$. Thus, using Eq. (E.4) for $\ell \ge 1$,

$$ \alpha_{\pm} (\alpha_{\pm} - \alpha_v - 1) = \ell(\ell + 1) - M_{\text{sh}}^2 \left( \frac{\omega r_{\text{sh}}}{v_{\text{sh}}} \right)^2 \left( \frac{r}{r_{\text{sh}}} \right)^3, \quad (E.19) $$

$$ \sim \ell(\ell + 1) \quad \text{if } \omega \ll \omega_{\text{Lamb}}. \quad (E.20) $$

In this region the approximate solution is $\alpha_{\pm} \equiv a \pm b$ with $a, b$ defined by Eqs. (75) and (76). The homogeneous solution $Y_0$ for $\ell \ge 1$ is a linear combination of power laws $Y_0^{\pm}$ satisfying the lower boundary condition (61):

$$ Y_0(r) = \left(\frac{r}{r_{\text{ns}}}\right)^{a-b} + \mathcal{R}_{\text{ns}} \left(\frac{r}{r_{\text{ns}}}\right)^{a+b}, \quad (E.21) $$

$$ \frac{1 - M_{\text{ns}}^2}{v_{\text{ns}}} \left( \frac{\partial Y_0}{\partial r} \right)_{\text{ns}} = \frac{i\omega}{c_{\text{ns}}^2} Y_0^{\text{ns}}. \quad (E.22) $$

With $M_{\text{ns}} \ll 1$ and using Eq. (72),

$$ a - b + (a + b)\mathcal{R}_{\text{ns}} = M_{\text{sh}}^2 \left(\frac{r_{\text{ns}}}{r_{\text{sh}}}\right)^{2+\alpha_v} \frac{i\omega r_{\text{sh}}}{v_{\text{sh}}} (1 + \mathcal{R}_{\text{ns}}), \quad (E.23) $$

$$ \mathcal{R}_{\text{ns}} = \frac{a - b - M_{\text{sh}}^{\frac{4}{\gamma+1}} \left( \frac{r_{\text{ns}}}{r_{\text{sh}}} \right)^{\frac{\gamma+3+\alpha_v(3-\gamma)}{\gamma+1}} \frac{i\omega r_{\text{sh}}}{v_{\text{sh}}}}{a + b + M_{\text{sh}}^{\frac{4}{\gamma+1}} \left( \frac{r_{\text{ns}}}{r_{\text{sh}}} \right)^{\frac{\gamma+3+\alpha_v(3-\gamma)}{\gamma+1}} \frac{i\omega r_{\text{sh}}}{v_{\text{sh}}}}. \quad (E.24) $$

From Eq. (E.24) we conclude that the coefficient $\mathcal{R}_{\text{ns}}$ is asymptotically independent of $\omega$ when $r_{\text{ns}} \ll r_{\text{sh}}$.

$$ \mathcal{R}_{\text{ns}} \sim \frac{a - b}{a + b}. \quad (E.25) $$

The acoustic function is independent of the frequency in the asymptotic limit $r_{\text{ns}} \ll r_{\text{sh}}$ for low frequency perturbations $\ell \ge 1$ driven by advection such that $\omega r_{\text{sh}}/|v_{\text{sh}}| \le 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 $\delta S = 0$ and Eq. (59) that the relation between $\partial Y_0/\partial r$ and $\delta v_r^0$ is

$$ \delta v_r^0 = - \frac{1}{\ell(\ell + 1)} \frac{\partial \delta A_0}{\partial r}, \quad (E.26) $$

$$ = \left( \frac{\partial Y_0}{\partial r} - \frac{i\omega M^2 Y_0}{1 - M^2 v} \right) e^{-\int_{\text{sh}}^{r} \frac{i\omega M^2}{1-M^2} \frac{dr}{v}}. \quad (E.27) $$

The fact that Eq. (74) imposes $\partial Y_0/\partial r = 0$ is approximately compatible with $\delta v_r^0 \sim 0$ to the extent that $M_{\text{ns}} \ll 1$.

## Appendix F: Approximation of the eigenfrequency equation for $\gamma = 4/3$

The advection time for $\alpha_v = 1$ is approximated using Eq. (79):

$$ e^{\int_{\text{sh}}^{r} \frac{i\omega}{v_r} dr} = \left( \frac{x}{x_{\text{sh}}} \right)^{\frac{i\omega r_{\text{sh}}}{|v_{\text{sh}}|}}. \quad (F.1) $$

Defining $\delta z \equiv i\omega r_{\text{sh}}/|v_{\text{sh}}| + 2 - b$ and noting that $a = 1$, the leading terms in Eq. (81) for $\ell \ge 1$ are:

$$ N(\delta z - 2 + b) + \frac{2b M_{\text{sh}}^2}{\ell(\ell + 1)} x_{\text{sh}}^{\delta z - 3} = \int_{1}^{x_{\text{sh}}} \left( \frac{x_{\text{sh}}}{x} \right)^{\delta z} \left( \frac{1}{x^{2b}} - 1 \right) \frac{dx}{x}, \quad (F.2) $$

The integral on the right hand side can be calculated explicitly. If $\delta z \ne 0$,

$$ \int_{1}^{x_{\text{sh}}} \left( \frac{x_{\text{sh}}}{x} \right)^{\delta z} \left( \frac{1}{x^{2b}} - 1 \right) \frac{dx}{x} = \frac{1}{\delta z} \left( 1 - \frac{\delta z x_{\text{sh}}^{-2b}}{2b + \delta z} - \frac{2b x_{\text{sh}}^{\delta z}}{2b + \delta z} \right), \quad (F.3) $$

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

$$ [1 - \delta z N(b - 2 + \delta z)] \left(1 + \frac{\delta z}{2b}\right) - \frac{\delta z}{2b x_{\text{sh}}^{2b}} = x_{\text{sh}}^{\delta z} \left[ 1 + \delta z \frac{2b M_{\text{sh}}^2}{\ell(\ell + 1) x_{\text{sh}}^3} \left(1 + \frac{\delta z}{2b}\right) \right]. \quad (F.4) $$
Reintroducing the normalised eigenfrequency $Z = i\omega r_{sh}/|v_{sh}|$,

$$ \frac{x_{sh}^{b-2}}{2b} \left\{ [1 - (Z + 2 - b) N(Z)] (Z + 2 + b) - \frac{Z + 2 - b}{x_{sh}^{2b}} \right\} = x_{sh}^{Z} \left\{ 1 + \left[ (Z + 2)^2 - b^2 \right] \frac{M_{sh}^2}{\ell(\ell + 1)x_{sh}^3} \right\}. \quad (F.5) $$

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