A&A 392, 353-368 (2002)
DOI: 10.1051/0004-6361:20020912

Non-radial instabilities of isothermal Bondi accretion
with a shock: Vortical-acoustic cycle vs. post-shock acceleration

T. Foglizzo

Service d'Astrophysique, CEA/DSM/DAPNIA, CEA-Saclay, 91191 Gif-sur-Yvette, France

Received 17 May 2002 / Accepted 12 June 2002

Abstract
The linear stability of isothermal Bondi accretion with a shock is studied analytically in the asymptotic limit of high incident Mach number ${\cal M}_{1}$. The flow is unstable with respect to radial perturbations as expected by Nakayama (1993), due to post-shock acceleration. Its growth-time scales like the advection time from the shock ${r}_{\rm sh}$ to the sonic point $r_{\rm son}$. The growth rate of non-radial perturbations l=1 is higher by a factor ${\cal M}_{1}^{2/3}$, and is therefore intermediate between the advection and acoustic frequencies. Besides these instabilities based on post-shock acceleration, our study revealed another generic mechanism based on the cycle of acoustic and vortical perturbations between the shock and the sonic radius, independently of the sign of post-shock acceleration. The vortical-acoustic instability is fundamentally non-radial. It is fed by the efficient excitation of vorticity waves by the isothermal shock perturbed by acoustic waves. The growth rate exceeds the advection frequency by a factor $\log{\cal M}_{1}$. Unstable modes cover a wide range of frequencies from the fundamental acoustic frequency $\sim$ $c/{r}_{\rm sh}$ up to a cut-off $\sim$ $c/r_{\rm son}$ associated with the sonic radius. The highest growth rate is reached for l=1 modes near the cut-off. The additional cycle of acoustic waves between the shock and the sonic radius is responsible for variations of the growth rate by a factor up to 3 depending on its phase relative to the vortical-acoustic cycle. The instability also exists, with a similar growth rate, below the fundamental acoustic frequency down to the advection frequency, as vorticity waves are efficiently coupled to the region of pseudosound. These results open new perspectives to address the stability of shocked accretion flows.

Key words: accretion, accretion disks - hydrodynamics - instabilities - shock waves - stars binaries: close - X-rays: stars

1 Introduction

Hydrodynamic instabilities in accretion flows can help our understanding of the variability observed in the luminosity of X-ray binaries. Numerical simulations have revealed the existence of such a hydrodynamical instability in the accretion flow of a gas on a compact accretor moving at supersonic velocity (Bondi-Hoyle-Lyttleton accretion). A first step towards understanding the physical mechanism underlying this instability was made by Foglizzo & Tagger (2000, hereafter FT00) who recognized the unstable cycle of entropic and acoustic perturbations between the shock and the sonic surface. This cycle is unstable if there is a large enough temperature difference between the shock and the sonic surface. The academic case of shocked Bondi accretion was studied by Foglizzo (2001, hereafter F01) who revealed the importance of non radial perturbations, and vorticity in particular. Both vorticity and entropy perturbations are advected from the shock to the accretor, and both are coupled to the acoustic perturbations. This coupling was formulated in a compact way by Howe (1975). If the adiabatic index is in the range $1<\gamma\le5/3$, entropy and vorticity perturbations are intimately related in the shocked Bondi flow (Foglizzo 2002, in preparation), so that it is difficult to identify their respective roles in the instability mechanism. By contrast, their roles are well separated in the isothermal limit ($\gamma=1$), where entropy perturbations are absent from the problem. The present paper is therefore dedicated to the study of the linear stability of shocked accretion in the isothermal Bondi flow, where the incident Mach number is the only parameter. The stability of shocked isothermal flows was studied by Nakayama (1992, 1993) in the more general context of flows with small angular momentum. The study of Nakayama was restricted to axisymmetric perturbations, thus precluding any possible vortical acoustic cycle. In this approximation, Nakayama analytically obtained the result that the flow is unstable if the flow accelerates immediately after the shock surface. In this respect the shocked Bondi flow should be unstable. We are therefore interested here in an extension of Nakayama's results to the case of non radial perturbations, in order to include the effect of vorticity. It should be noted that isothermal flows are very particular as far as Bondi-Hoyle-Lyttleton (hereafter BHL) accretion is concerned, with unsettled issues concerning the influence of the numerical resolution. The instability observed in numerical simulations seems very weak in 3-D according to Ruffert (1996), whereas it is violent in 2-D (Ishii et al.1993; Shima et al.1998). Pogorelov et al.(2000) discussed the possible responsability of the numerical procedure in producing the instability. With a different approach, the instability of the shocked Bondi flow described in this paper could contribute to guide our physical understanding of more complex flows involving isothermal shocks.

The paper is organized as follows. Perturbed equations are described in Sect. 2 and eigenfrequencies are determined numerically in Sect. 3. Analytical methods are used to disentangle the effect of advection from the effect of the boundary. In the spirit of F01, the coupling between the vorticity and acoustic perturbations in the classical isothermal Bondi flow, without a shock, is described in Sect. 4. Boundary conditions are taken into account to build a vortical acoustic-cycle in Sects. 5 and 6. The instability due to post-shock acceleration is analysed in Sect. 7. The relationship between the vortical-acoustic instability and existing numerical simulations of BHL accretion and shocked discs is discussed in Sect. 8.

  
2 Linearized equations of the shocked Bondi flow

  
2.1 Properties of the unperturbed shocked flow

	Figure 1: Radial profile of the Mach number in the unperturbed shocked flow, for M1=5 (full line). The subsonic cavity stands between the sonic point and the shock radius. The dashed lines show the analytical continuation of the solutions beyond the shock radius.
Open with DEXTER

The Mach number in the unperturbed transonic flow satisfies the following equation deduced from the conservation of the mass flux and Bernoulli constant B and the regularity at the sonic radius:

$${r\over r_{\rm son}}{\cal M}^{1\over2}\exp\left({r_{\rm son}\over r}-{{\cal M}^{2}\over 4}\right)= {\rm e}^{3\over4}.$$ (1)

The sonic radius is half of the Bondi radius GM/c². The Bernoulli constant,

$$B\equiv {v^{2}\over2} + c^{2}\log\rho -{2r_{\rm son}\over r},$$ (2)

is conserved along the flow lines. By contrast with the adiabatic case $\gamma>1$, this quantity is not conserved through a shock. The shock radius ${r}_{\rm sh}$ corresponding to an incident Mach number ${\cal M}_{1}$ is determined from Eq. (1) together with the Rankine Hugoniot jump condition ${\cal M}_{2}=1/{\cal M}_{1}$. The Mach number profile is shown in Fig. 1. It should be noted that the presence of a radial shock is not guaranteed a priori, since the supersonic preshock flow ( $r>{r}_{\rm sh}$) could be continued without a shock down to the accretor (dashed line in Fig. 1). Several physical mechanisms could trigger the formation of a shock, such as the heating of protons to temperatures at which the fluid becomes collisionless (Mezaros & Ostriker 1983), the trapping of relativistic particules (Protheroe & Kazanas 1983) or the dissipation of magnetic fields (McCrea 1956; Scharlemann 1981). In the context of Bondi-Hoyle-Lyttleton accretion, the existence of a stagnation point behind the accretor implies that a fraction of the supersonic flow decelerates to subsonic velocities, and therefore naturally produces a shock. The shock radius deduced from Eq. (1) increases with ${\cal M}_{1}$ for spherical accretion:

$${{r}_{\rm sh}\over r_{\rm son}}\sim {\rm e}^{3\over4}{\cal M}_{1}^{1\over2}.$$ (3)

Compared to BHL accretion, the presence of an accretion shock at a distance exceeding the accretion radius $2GM/v_{\infty}^2$ is rather artificial and due simply to the assumption of purely radial velocity. In the highly supersonic limit ${\cal M}_{1}\gg1$ initially considered by Hoyle & Lyttleton (1939), any trajectory with an angular momentum larger than $\sim$2 $GM/v_{\infty}$ would of course miss the accretor. Bearing this feature in mind, the case of strong shocks studied in this paper is still very useful in order to understand the mechanisms involved, because of the separation of timescales it enables.

  
2.2 Scaling of timescales in the shocked Bondi flow

If the shock is strong, the timescale $\tau_{\rm adv}$ associated with advection is much longer than the sound crossing time ${r}_{\rm sh}/c$:

$$\tau_{\rm adv}\equiv\int^{{r}_{\rm sh}}_{r_{\rm son}} {{\rm d... ...}\over 3}{{r}_{\rm sh}\over c}\gg {2{r}_{\rm sh}\over c} \cdot$$ (4)

The third fundamental timescale of the flow is related to the presence of the sonic radius $\sim$ $r_{\rm son}/c$. Together with Eq. (3), the scaling of the three special timescales of the problem is thus the following for strong shocks:

$$\tau_{\rm adv}\gg{{r}_{\rm sh}\over c}\gg {r_{\rm son}\over c}\cdot$$ (5)

2.3 Differential system for the perturbations of the Bondi flow

Figure 2: Schematic structure of the perturbed acoustic (full line) and vorticity (dotted line) fields, depending on the frequency of the perturbation. The curve beyond the shock radius shows the acousticperturbation in the Bondi flow without a shock. The turning point $r_{\rm t}$ separates the region of pseudosound from the region of propagation of acoustic pertrbations. The vortical-acoustic coupling is most efficient in the region between $r_{\rm son}$ and $r_{\rm eff}$.

Open with DEXTER

As in the case $\gamma>1$ studied by F01, the vorticity equation can be directly integrated, so that the three quantities ( $r^{2}\delta w_{\rm r},rv\delta w_{\theta},rv\delta w_{\varphi}$) are conserved when advected. The perturbations of the radial velocity $\delta v_{\rm r}$ and density $\delta \rho$ in the isothermal Bondi flow are conveniently described by the functions f, g defined by:

  

$$\equiv v\delta v_{\rm r}+c^{2}{\delta\rho\over\rho},$$ (6)

$$\equiv {\delta v_{\rm r}\over v}+{\delta\rho\over\rho}\cdot$$ (7)

They satisfy the same differential equations as the functions f, g in the Appendix B of F01, by suppressing the entropy terms:

  

$$v{\partial f\over\partial r}+{i\omega {\cal M}^2f\over 1-{\cal M}^2} {i\omega v^2 g\over 1-{\cal M}^2},$$ = (8)

$$v{\partial g\over\partial r}+{i\omega {\cal M}^2g\over 1-{\cal M}^2} {i\omega \mu^{2}f\over c^2(1-{\cal M}^2)} + {i\delta K_{\rm sh}\over r^2\omega }{\rm e}^{i\omega\int_{{r}_{\rm sh}}^r {{\rm d} r\over v}},$$ = (9)

where the function $\mu(r,\omega,l)$ is defined by:

  

$$\mu^{2} \equiv 1 - {\omega_{l}^{2}\over \omega^{2}},$$ (10)

$$\omega_{l}^2 \equiv {l(l+1)\over r^2}(c^{2}-v^{2}).$$ (11)

The constant $\delta K$ is related to vorticity as follows

$$\delta K\equiv r^{2}v\cdot(\nabla\times w).$$ (12)

By contrast with the case of adiabatic flows (F01), $\delta K$ is the unique source term for the excitation of acoustic waves. For radial perturbations ( $l=0, \delta K=0$), the equations analysed by Nakayama (1992, 1993), with $\Psi\equiv f/i\omega$, are recovered.

The position of the turning point $r_{\rm t}$ of non-radial acoustic waves is defined by $\mu=0$ in Eq. (10). This turning point lies inside the subsonic cavity if $\omega>\omega_{\rm sh}$, with

$$\omega_{\rm sh}\equiv l^{1\over2}(l+1)^{1\over2}(1-{\cal M}_{... ...over 2} {c\over {r}_{\rm sh}}\propto {c\over{r}_{\rm sh}}\cdot$$ (13)

Below $\omega _{\rm sh}$ (or inside the turning point $r_{\rm t}$ if it exists), the acoustic perturbation is not propagating. It was named "pseudosound'' by Ffowcs Williams (1969).

The threshold between absorbed and trapped sound was introduced in F01 through the cut-off frequency $\omega_l^{\rm cut}$, which is independent of the shock strength:

$$\omega_l^{\rm cut}\sim {l^{1\over2}(l+1)^{1\over2}\over2}{c\over r_{\rm son}}\propto {c\over r_{\rm son}}\cdot$$ (14)

This cut-off frequency used to be refered to as a "refraction'' cut-off in F01 for $\gamma>1$, due to the non-homogeneity of the sound speed which bends trajectories outwards. By contrast, in the isothermal flow considered here, acoustic trajectories are only bent inwards, due to the sole effect of flow acceleration. The existence of a cut-off for non-radial acoustic waves can be understood by considering the frequency dependence of the direction $\psi$ of propagation of the incoming acoustic flux relative to the radial direction. This angle is approximated in Appendix E in the WKB approximation (Eq. (E.11)). The higher the frequency, the smaller $\psi$, so that high frequency waves with a given order l are pointing towards the accretor and are absorbed inside the sonic sphere.

The scaling of frequencies stressed in Eq. (5) for strong shocks consequently separates three ranges of eigenmodes schematized in Fig. 2:
(i) absorbed sound $\omega_l^{\rm cut}<\omega$,
(ii) trapped sound $\omega_{\rm sh}<\omega<\omega_l^{\rm cut}$,
(iii) pseudosound $1/\tau_{\rm adv}<\omega<\omega_{\rm sh}$. The corresponding structure of the vorticity wave is also displayed in Fig. 2, in each range of frequencies.

2.4 Boundary condition at the shock

The boundary conditions at the shock associated with the variables f,g can be computed with the same method as Nakayama (1992), extended to the case of non radial perturbations. Let us denote by $\Delta\zeta$ the radial displacement of the shock produced by a sound wave propagating against the stream in the subsonic region, and $\Delta v=-i\omega\Delta\zeta$ its velocity. The incident Mach number ${\cal M}_{1}'$ in the frame of the shock is altered by the perturbations of velocity and displacement as folllows:

 

$${\cal M}_{1}'({r}_{\rm sh}+\Delta\zeta) {\cal M}_{1}({r}_{\rm sh})+{\Delta v\over c}+ \Delta\zeta{\partial{\cal M}_{1}\over\partial r},$$ = (15)

$${\cal M}_{1}({r}_{\rm sh})+\left(1-{i\eta c{\cal M}_{2}\over\omega{r}_{\rm sh}}\right) {\Delta v\over c},$$ = (16)

where the parameter $\eta$ measures the strength of the local gradient of the Mach number immediately after the shock:

$$\eta\equiv {\partial\log {\cal M}_{2}\over\partial \log r}=-{... ...cal M}_{2}^{2}} \left(1-{r_{\rm son}\over{r}_{\rm sh}}\right),$$ (17)

The boundary conditions $f_{\rm sh},g_{\rm sh}$, written as functions of $\Delta v$, are deduced from the conservation of the mass flux and impulsion across the shock:

  

$$f_{\rm sh} -c^{2}{\cal M}_{2}(1-{\cal M}_{2}^{2}) \left\lbrack(1-{\cal M}_{2... ...\eta c\over\omega{r}_{\rm sh}}+{\cal M}_{2}\right\rbrack {\Delta v\over v_{2}},$$ (18)

$$g_{\rm sh} (1-{\cal M}_{2}^{2}){\Delta v\over v_{2}}\cdot$$ (19)

The non radial perturbation of the velocity, and the perturbed vorticity computed in Appendix A enable us to relate $\delta K_{\rm sh}$ to $\Delta v$:

$$\delta K_{\rm sh}=l(l+1)c^{2}(1-{\cal M}_{2}^{2})^{2} \left( ... ..._{2}\over\omega{r}_{\rm sh}}\right){\Delta v\over v_{2}}\cdot$$ (20)

Equation (20) can be used together with Eqs. (18)-(19) in order to express $f_{\rm sh},g_{\rm sh}$ as functions of $\delta K_{\rm sh}$ for non-radial perturbations. In Appendix C the boundary value problem is reduced to a single equation incorporating the boundary conditions both at the shock and at the sonic point. The eigenfrequencies satisfy the following equation, where f₀ is the unique homogeneous solution which is regular at the sonic point:

 

$${\partial^{2} f_{0}\over\partial r^{2}}+\left\lbrack {\partial\lo... ...al r} +{{\omega^{2}\over c^{2}}-{L^{2}\over r^{2}}\over 1-{\cal M}^{2}}f_{0}=0,$$ (21)

 

$${\left\lbrack (1-{\cal M}_{2}^{2})\eta - {\cal M}_{2}{i\omega {r}_{\rm sh}\over c}\right\rbrack {\partial f_{0}\over\partial r}({r}_{\rm sh})}$$

$${+{i\omega {r}_{\rm sh}\over c}\left\lbrack {1+{\cal M}_{2}^{2}\o... ...ver c} -{\cal M}_{2}\eta\right\rbrack {f_{0}({r}_{\rm sh})\over {r}_{\rm sh}}=}$$

$${l(l+1) \left({i\omega {r}_{\rm sh}\over c}+{\cal M}_{2}\eta\righ... ..._{\rm sh}}^r{1+{\cal M}^{2}\over 1-{\cal M}^{2}}{{\rm d} r\over v} }{\rm d} r.}$$ (22)

Equation (22) is independent of the normalization of f₀. It should be noted that the integral on the right hand side is well defined despite the singularity of the phase near the sonic radius. This equation describes the perturbation of the shock by the interplay of the acoustic and vortical perturbations. The same calculation in any other potential (e.g.the Paczynski-Wiita potential) would lead to the same system of Eqs. (21)-(22), only the shape of ${\cal M}$ (and its derivative $\eta$) described by Eq. (1) would be affected. The left hand side of Eq. (22) involves only the acoustic perturbation f₀, and is independent of the vortical perturbation. By contrast, the integral on the right hand side describes the acoustic feed back of the vortical perturbation coupled to the acoustic field.

   
3 Spectrum of eigenfrequencies

3.1 Numerical procedure to determine the spectrum of eigenfrequencies

The regularity of the solution at the sonic radius has already been discussed in F01: the sonic point is a regular singularity if $\gamma<5/3$. For a given incident Mach number ${\cal M}_{1}$, numerical intergration is performed from the sonic point towards the shock in order to determine the eigenfrequencies of perturbations with a latitudinal number l. For a given value of $\omega$, the unique solution which is regular at the sonic point is expanded in a Frobenius series in order to start the numerical integration away from the singularity. A Runge-Kutta algorithm is then used to simultaneously integrate four functions, namely $f_{0}({r}_{\rm sh})$, $g_{0}({r}_{\rm sh})$, the integral on the right hand side of Eq. (22) and the integral phase inside it. $\omega$ is varied and shooting is repeated until Eq. (22) is satisfied.

3.2 Purely growing radial instability

Figure 3: Growth rate of the purely growing radial instability, multiplied by the advection time (thick dashed line). The thin dashed lines correspond to the lower and upper bounds determined by Nakayama (1993). Anticipating on the calculations of the rest of the paper, the four other curves are analytical estimates of the growth rate of l=1 perturbations at high Mach number. The growth rate of the entropic-acoustic instability increases along the dotted line at low frequency(Eq. (74)), and is bounded by the two thin lines at high frequency (Eqs. (67)-(68)). The thick full line describes the most unstable branch (Eq. (85)).

Open with DEXTER

The radial instability found numerically in Fig. 3 is consistent with the result of Nakayama (1992, 1993), who found that the local acceleration of the flow immediately after the shock is a source of non-oscillating radial instability. The growth rate is shown as a function of the incident Mach number, in units of the advection time $\tau_{\rm adv}$. The lower and upper bounds found analytically by Nakayama (1993, Eq. (27)) are also displayed in Fig. 3:

$${\vert\nu\vert\over1+{\cal M}_{1}^{2}}\le\omega_{i}\le{\vert\nu\vert\over1+{\cal M}_{1}},$$ (23)

where the parameter $\nu$ describes the acceleration of the flow and is proportionnal to $\eta$:

  

$$\nu \equiv -{1\over v_{2}}\left({\partial\Phi\over\partial r}-{2c^{2}\over{r}_{\rm sh}}\right),$$ (24)

$$({\cal M}_{1}^{2}-1){\eta c{\cal M}_{2}\over{r}_{\rm sh}}\cdot$$ (25)

Note that Eq. (24) is corrected for a factor of 2 for spherical flows, as noticed by Nakayama (1993). The lower and upper bounds of Eq. (23) cover a wide range of frequencies from the advection frequency $\sim$ $2/(3\tau_{\rm adv})$ up to the acoustic frequency $\sim$ $2c/{r}_{\rm sh}$. In the present case of isothermal radial accretion, the growth rate computed numerically converges precisely towards the value $2/\tau_{\rm adv}$.

3.3 Non-radial oscillatory instabilities

Figure 4: Real and imaginary parts of the eigenfrequencies of the modes l=0,1,2,3 in the isothermal flow, for ${\cal M}_1=100$ ( $c\tau _{\rm adv}/4\pi {r}_{\rm sh}\sim 2.8$). The real part is measured in units of the cut-off frequency $\omega _1^{\rm cut}$ (lower axis) and advection frequency (upper axis). The imaginary part is normalized to the advection time. Analytical estimates of the growth rates of the modes l=0,1,2 at low freqency are indicated using big symbols (Eqs. (78), (83) and (85) in Sects. 6 and 7).

Open with DEXTER

Figure 5: Real and imaginary parts of the eigenfrequencies of the modes l=1 in the isothermal flow, for different Mach numbers. Frequencies are multplied by the advection time. The strong l=1 instability due to postshock acceleration corresponds to the thick line in both plots. In the upper plot, the growth rate of the fastest of the mode l=0 is also displayed (thick dashed line). The analytical approximation of the fastest growth rates $\omega _{i}^{+},\omega _{i}^{\rm max}$ of the l=1instabilities due to the entropic-acoustic cycle and to post-shock acceleration are shown as the thin dashed line (Eqs. (68)) and the thin dotted line Eq. (85). In the bottom plot, the cut-off frequency $\omega _1^{\rm cut}$, the minimum acoustic frequency $\omega _{\rm sh}$, and the frequency $\omega _{r}^{\rm max}$ of the most unstable mode (Eq. (84)) are indicated as references.

Open with DEXTER

The effect of non radial perturbations ($l\ge1$) involves an integral on the right hand side of Eq. (22), which reflects the cumulative excitation of acoustic waves by the vorticity perturbations advected from the shock to the sonic radius. Figure 4 shows the eigenspectrum of the isothermal Bondi flow for M₁=100, l=0,1,2,3. This spectrum is characterized as follows:

(i) The most unstable mode is asymmetric (l=1) and well separated from the other eigenmodes. The growth rate of this low frequency mode is significantly larger than the advection time, nearly four times faster than the unstable radial mode. The analytical estimates of the growth rates at low frequency correspond to Eqs. (77)-(78), (81)-(83) and (84)-(85) in Sects. 6 and 7.

(ii) Non radial unstable modes are very numerous and cover a wide range of frequencies, from the advection frequency ($\sim$ $v_{\rm sh}/{r}_{\rm sh}$), up to the cut-off frequency ( $\propto c/r_{\rm son}$).

(iii) A striking feature of the eigenspectrum in Fig. 4 is the apparent oscillation of the imaginary part as a function of the real part of the frequency. This behaviour is explained in Sect. 5 on the basis of the vortical-acoustic cycle.

The variation of the growth rate with the incident Mach number is shown in Fig. 5. The higher the Mach number, the more unstable the non radial modes, while the radial instability saturates. The growth rate of the most unstable l=1 mode increases with the incident Mach number much faster than in the other modes. The real part of its complex frequency is intermediate between the advection and acoustic frequencies. It is computed analytically in Sect. 7.

The rest of the paper is dedicated to understanding this apparently complicated spectrum, and to questioning the role of post-shock acceleration.

  
4 The sound of vorticity

4.1 Frequency dependence of the vortical-acoustic coupling

Figure 6: Efficiency $\vert{{\cal Q}}_{\rm adv}\vert$ of the excitation of acoustic waves propagating against the stream by the advection of a vortex in an isothermal flow, for l=1,2,3. The frequency is in units of the cut-off frequency $\omega _1^{\rm cut}$. The dotted lines correspond to the asymptotic behaviour ${{\cal Q}}_{\rm adv}\propto \omega ^{(2l-1)/3}$ at low frequency determined analytically in Eq. (30).

Open with DEXTER

Let us assume that a vorticity wave is advected from infinity towards the accretor: what is the sound produced by the advection of this vorticity? This question is solved in Appendix D by computing the coupling coefficient ${\cal Q}_{\rm adv}$ measuring the outgoing acoustic flux F_- associated with a vorticity perturbation $\delta K$:

$$F_{-}={{\dot M}_{0}\over c^2}\vert{{\cal Q}}_{\rm adv}\vert^{2} \vert\delta K\vert^{2},$$ (26)

The relationship between the acoustic flux $F_{\pm}$ and the acoustic perturbation $f_{\pm}$ is deduced from the WKB approximation in Appendix B:

$$F_{\pm}\sim {{\dot M}_{0}\over c^{2}}{\mu\over{\cal M}}\left\vert f_{\pm}\right\vert^{2}.$$ (27)

Thus the complex efficiency ${\cal Q}_{\rm adv}$ is defined at a radius $R\gg c/\omega$, in the region of propagation of acoustic waves, as follows:

  

$${{\cal Q}}_{\rm adv}(R) \equiv \left({\mu\over{\cal M}}\right)^{1\over2} {f_{-}\over\delta K },$$ (28)

$$-{1\over 2i\omega c} {\rm e}^{{i\omega\over c}\int_{r_{\rm t}}^{R... ...t_{R}^r{1+{\cal M}^{2}\over 1-{\cal M}^{2}}{{\rm d} r\over{\cal M}}} {\rm d} r.$$ = (29)

${\cal Q}_{\rm adv}$ is approximately independent of R in the WKB region of propagation of acoustic waves, i.e.away from their turning point. This coefficient is approximated in Appendix D in the low frequency limit $c/R\ll\omega\ll c/r_{\rm son}$, i.e.below the cut-off frequency, using the function ${\cal H}(l)$ defined from an integral of the Spherical Bessel function j_l (Eq. (D.7)):

$$\vert{\cal Q}_{\rm adv}\vert\sim {3^{l+1\over3}{\rm e}^{2l-1\... ...)\over2}\left({\omega r_{\rm son}\over c}\right)^{2l-1\over3},$$ (30)

with ${\cal H}(1)\sim 0.45$, ${\cal H}(2)\sim 0.33$ and ${\cal H}(3)\sim 0.28$. The direct integration of Eq. (29) is compared to the approximation obtained in Eq. (30) for l=1,2,3 in Fig. 6. The asymptotic behaviour $\vert{\cal Q}_{\rm adv}\vert\propto \omega^{2l-1/3}$is correctly reproduced. The scaling factor of this power law obtained analytically is too low by a factor 2 for l=1 perturbations, which can be attributed to the roughness of the approximation. This calculation clearly indicates that the vortical-acoustic coupling is most efficient near the cut-off frequency. The same conclusion was reached in F01 concerning the efficiency of the entropic-acoustic coupling.

  
4.2 Region of coupling

This approximation enables us to determine the region where the coupling between the vorticity and the acoustic perturbations is most efficient. It occurs mainly in the subsonic region $r_{\rm son}<r<r_{\rm eff}$ where the wavelength of the vorticity perturbation is largest (see Fig. 2):

$$\omega \int_{r_{\rm son}}^{r_{\rm eff}} {{\rm d} r\over v}\sim 2\pi.$$ (31)

A similar result was obtained in F01 concerning entropy perturbations for $\gamma=5/3$. Due to the differences in velocity profiles in Eq. (31), $r_{\rm eff}\propto \omega^{-2/3}$ for $\gamma=5/3$ at high frequency (Eqs. (E.6), (E.11) and (E.12) of F01), whereas

$$r_{\rm eff}\propto \left({\omega r_{\rm son}\over c}\right)^{-{1\over3}}r_{\rm son}$$ (32)

in the isothermal flow at low frequency. It is therefore important to note that the most efficient coupling at low frequency lies in the region of pseudosound well within the turning point ( $r_{\rm t}\sim c/\omega$) of acoustic waves:

$$r_{\rm son}\ll r_{\rm eff}\ll r_{\rm t}.$$ (33)

According to Eq. (31), the upper bound $r_{\rm eff}$ of the coupling region for perturbations near the advection frequency coincides with the shock radius.

  
5 Vortical-acoustic instability at high frequency

5.1 The vortical-acoustic cycle for $\omega _{\rm sh}\ll \omega \ll \omega _l^{\rm cut}$

The formalism developed in FT00 for the entropic-acoustic cycle within the WKB approximation can be transposed to the case of the vortical-acoustic cycle. This formalism requires the frequency to be high enough so that acoustic waves can be identified and their propagation time measured. An initial perturbation of vorticity $\delta K(t_{0},{r}_{\rm sh})$ at the shock, advected towards the accretor, triggers the excitation of an acoustic flux F_- propagating outwards. As it reaches the shock surface, this acoustic flux induces new vorticity perturbations $\delta K(t_{0}+\tau_{\cal Q},{r}_{\rm sh})$, where  $\tau_{\cal Q}$ is the duration of this cycle. The global efficiency ${\cal Q}$ of the cycle is naturally

$${\cal Q}\equiv {\delta K(t_{0}+\tau_{\cal Q},{r}_{\rm sh})\over\delta K(t_{0},{r}_{\rm sh})}\cdot$$ (34)

${\cal Q}$ and $\tau_{\cal Q}$ depend a priori on the frequency $\omega_{r}$ and spatial structure l of the perturbation considered. The growth (or damping) rate of the cycle is identified with the imaginary part of the eigenfrequency, given by:

$$\omega_{i}\equiv {1\over\tau_{\cal Q}}\log \vert{\cal Q}\vert.$$ (35)

The timescale $\tau_{\cal Q}$ is dominated by the advection timescale $\tau_{\rm adv}$ if the shock is strong ( ${\cal M}_{2}\ll1$, ${r}_{\rm sh}\gg r_{\rm t}$):

$$\tau_{\cal Q}\sim\int_{r_{\rm t}}^{{r}_{\rm sh}} {1\over1-{\cal M}}{{\rm d} r\over \vert v\vert} \sim \tau_{\rm adv} .$$ (36)

The global efficiency ${\cal Q}$ can be decomposed in terms of the efficiencies of the coupling between F_- and $\delta K$ during advection ( ${\cal Q}_{\rm adv}$) and at the shock ( ${\cal Q}_{\rm sh}$), defined immediately after the shock according to Eq. (28) and:

  

$${{{\cal Q}}_{\rm sh}\equiv \left({{\cal M}_{2}\over\mu}\right)^{1\over2} {\delta K_{\rm sh} \over f_{-}},}$$ (37)

$${{\cal Q}= {{\cal Q}}_{\rm adv}{{\cal Q}}_{\rm sh}{\rm e}^{-i\omega\tau_{\cal Q}}.}$$ (38)

The modulus of ${\cal Q}_{\rm sh}$is directly interpreted as a coupling efficiency between F_- and $\vert\delta K\vert^{2}$:

$$\vert\delta K_{\rm sh}\vert^{2}= {c^{2}\over{\dot M}_{0}} \vert{{\cal Q}}_{\rm sh}\vert^{2}F_{-}.$$ (39)

The efficiency ${\cal Q}_{\rm sh}$ is computed in Appendix E, following a method introduced by D'Iakov (1958) and Kontorovich (1958, 1959). ${\cal Q}_{\rm sh}$ can be approximated by a WKB analysis when the wavelength of the perturbation is small compared to the lengthscale of the local gradients of the flow, by keeping only the first order terms in ${\cal M}_{2}$ and $c/\omega r$, at high frequency and for strong shocks:

$${{\cal Q}}_{\rm sh}\sim-{2l(l+1)\over{\cal M}_{2}^{1\over2}} \left(1-{\cal M}_{2}- i{\eta c\over\omega r}\right)\cdot$$ (40)

This calculation proves that the flow reacceleration does not affect the vortical-acoustic coupling at the shock for

$${\omega r\over c}\gg \left\vert\eta\right\vert\sim 2,$$ (41)

and that ${\cal Q}_{\rm sh}$ is approximately independent of frequency in the range $\omega_{\rm sh}\ll\omega<\omega_l^{\rm cut}$. As a consequence of Eqs. (37) and (40), the vorticity is produced at the shock with an efficiency proportional to the shock strength:

$$\delta K_{\rm sh}\propto {\cal M}_{1} f_{-}\propto {\cal M}_{1}{\delta p_{-}\over p}\cdot$$ (42)

The global efficiency ${\cal Q}$ should consequently peak near the cut-off frequency, and drop abruptly above it. The frequency dependence of ${\cal Q}$ appears in Fig. 8 (top picture, full line), as deduced from Fig. 6 and Eq. (40) in the WKB approximation. Thus the instability should be strongest near the cut-off frequency with the following scaling, obtained from Eqs. (14), (30) and (40) for $\omega _{\rm sh}\ll \omega \ll \omega _l^{\rm cut}$:

$$\vert{\cal Q}\vert\sim \left({\omega\over\omega_l^{\rm cut}}\... ...er 3} \left({{\cal M}_{1}\over{\cal M}_{0l}}\right)^{1\over2},$$ (43)

where the scaling factor ${\cal M}_{0l}$ depends on l through both ${\cal Q}_{\rm adv}$ and ${\cal Q}_{\rm sh}$:

$${\cal M}_{0l}\sim \left({4\over3{\rm e}^{3\over2}}\right)^{2l... ...(l+1)\right\rbrack^{2l+5\over3}} {1\over3{\cal H}^{2}(l)}\cdot$$ (44)

${\cal M}_{0l}\propto l^{2l/3}$ diverges for $l\gg10$, thus favouring the instability of low degree modes. According to Fig. 6, the extrapolation of Eq. (43) to $\omega\sim\omega_l^{\rm cut}$ gives a very rough upper bound of $\vert{\cal Q}\vert(\omega_l^{\rm cut})$, particularly overestimlating it for high degree perturbations. Although ${\cal M}_{01}\sim1.96$ exceeds ${\cal M}_{02}\sim 0.95$ and ${\cal M}_{03}\sim0.69$, the actual efficiency $\vert{\cal Q}\vert$ is maximal for the mode l=1, as checked in Fig. 8 by multiplying the curve $\vert{{\cal Q}}_{\rm adv}\vert$ (Fig. 6) by $\vert{\cal Q}_{\rm sh}\vert$ (Eq. (40)). Using Eqs. (35) and (36), the growth rate $\omega_{i}$ at high frequency is deduced from Eq. (43):

$$\omega_{\rm sh}\ll\omega_{r}\ll\omega_l^{\rm cut},$$ (45)

  

$$\omega_{i} \sim {1\over2\tau_{\rm adv}}\log \left({\omega\over\omega_l^{\rm cut}}\right)^{2(2l-1)\over 3} {{\cal M}_{1}\over{\cal M}_{0l}},$$ (46)

$$\sim {1\over2\tau_{\rm adv}}\log \left({2\omega\over{\rm e}^{3\over4}\... ...}}\right)^{2(2l-1)\over 3} {{\cal M}_{1}^{2(2-l)\over3}\over{\cal M}_{0l}}\cdot$$ (47)

Equation (47) indicates that the range of unstable frequencies below the cut-off frequency gets narrower for $l\ge 2$. This is also confirmed by the eigenspectrum obtained numerically for l=1,2,3 in Fig. 4. Note that the efficiency $\vert{\cal Q}\vert$ was estimated using $\vert{{\cal Q}}_{\rm adv}\vert$ which was computed for perturbations with a purely real frequency. Thus the estimate of the growth rate in Eq. (46) implicitly assumes that the imaginary part is small compared to the real part. This is true since $\log{\cal M}_{1}/\tau_{\rm adv}\ll \omega_l^{\rm cut}$. This calculation proves that the vortical-acoustic cycle is unstable for strong shocks, with a growth rate exceeding the growth rate of the radial mode, and a mechanism which is independent of the sign of the local flow acceleration immediately after the shock. A more refined description of the entropic-acoustic cycle is developped in the next section in order to understand the apparent oscillations in the spectrum obtained numerically in Figs. 4 and 5.

5.2 Additonnal contribution of the acoustic cycle

  
5.2.1 Dispersion relation for the double cycle

The role of the purely acoustic cycle had been anticipated in FT00. It is characterized by a time scale $\tau_{\cal R}$ and a global efficiency ${\cal R}$. Let us show how the simultaneous existence of the two cycles can explain the apparent oscillation in the eigenspectrum of Fig. 4. The following extension of the analysis of FT00 is valid for both entropic-acoustic and vortical-acoustic cycles, by replacing "vortical'' by "entropic''.

A perturbation f is influenced by the two cycles as follows:

$$f(t)={\cal Q}f(t-\tau_{Q}) + {\cal R}f(t-\tau_{R}).$$ (48)

A solution of the form $f(t)\propto \exp(-i\omega t)$ satisfies Eq. (48) if the complex eigenfrequency ( $\omega_{r},\omega_{i}$) is a solution of the following dispersion equation (Eq. (25) of FT00):

$${\cal Q}{\rm e}^{i\omega\tau_{\cal Q}}+{\cal R}{\rm e}^{i\omega\tau_{\cal R}}=1.$$ (49)

This dispersion relation is recovered in Appendix F as a WKB approximation of the exact dispersion relation (22). The analysis of the dispersion relation (49) in Appendix G enables us to extract physical information from the complicated eigenspectrum in Fig. 4, such as the ratio of timescales $\tau _{\cal Q}/\tau _{\cal R}$, and the dimensionless efficiencies $\vert{\cal Q}\vert$ and $\vert{\cal R}\vert$.

5.2.2 Ratio of the two cycles timescales $\tau _{\cal Q}/\tau _{\cal R}$

	Figure 7: Number of eigenmodes in each oscillation of the eigenspectrum in Fig. 4, identified as the ratio of the timescales of the two cycles. For strong shocks, $\tau _{\cal Q}/\tau _{\cal R}\sim {\cal M}_{1}/6$.
Open with DEXTER

The timescale $\tau_{\cal R}$ of the acoustic cycle for strong shocks can be written approximately:

$$\tau_{\cal R}\sim \int_{r_{\rm t}}^{{r}_{\rm sh}} {2\over 1-{... ... c}\sim 2{{r}_{\rm sh}\over c}\propto {\cal M}_{1}^{1\over2}.$$ (50)

According to Eqs. (36) and (50), the ratio $\tau _{\cal Q}/\tau _{\cal R}$ for a strong shock is simply:

$${\tau_{\cal Q}\over\tau_{\cal R}}\sim{{\cal M}_{1}\over6}\cdot$$ (51)

As shown in Appendix G, this ratio is the average number of eigenmodes in each oscillation of the eigenspectrum. As an illustration, this number is reported in Fig. 7 for the modes l=1,2,3 in the case ${\cal M}_1=100$ corresponding to Fig. 4, in good agreement with Eq. (51).

5.2.3 Efficiencies $\vert{\cal Q}\vert$ and $\vert{\cal R}\vert$

	Figure 8: Global efficiencies $\vert{\cal Q}\vert$, $\vert{\cal R}\vert$ asociated to non radial perturbations l=1,2,3 deduced from the eigenspectrum in Fig. 4 in the framework of the double cycle model (dotted line), compared with the efficiencies deduced from the calculation of ${{\cal Q}}_{\rm adv},{{\cal Q}}_{\rm sh}$ and ${{\cal R}}_{\rm adv},{{\cal R}}_{\rm sh}$.
Open with DEXTER

As noticed in FT00, the acoustic cycle can be neglected ( $\omega_{i} \sim\log \vert{\cal Q}\vert/\tau_{\cal Q}$) if the efficiencies and timescales of the two cycles are such that $\vert\alpha\vert\ll1$, with

$$\alpha\equiv {{\cal R}\over{\cal Q}^{\tau_{\cal R}\over\tau_{\cal Q}}}\cdot$$ (52)

If $\alpha$ is not negligible, the acoustic cycle can contribute to either stabilize or destabilize the vortical-acoustic cycle. The most stabilizing effect of the acoustic cycle is the effective reduction of ${\cal Q}$ by a factor 2 (see Appendix G). By contrast, its destabilizing contribution can be much larger comparatively, if the acoustic time is short compared to the advection time. This is the case for strong shocks. Depending on the relative phases of the two cycles, the growth rate can then cover the following range:

$${1\over\tau_{\cal Q}}\log {\vert{\cal Q}\vert\over 1+\vert\al... ...\cal Q}}\log {\vert{\cal Q}\vert\over 1-\vert\alpha\vert}\cdot$$ (53)

Conversely, the values of the dimensionless physical parameters $\vert{\cal Q}\vert,\vert{\cal R}\vert$ can be extracted from the eigenspectrum, as in Fig. 4, by measuring
(i) the extremal values $\omega_{i}^{+},\omega_{i}^{-}$,
(ii) the period $\Delta\omega_{r}$ of the oscillations
(iii) the average number n of modes per period.

As demonstrated in Appendix G, the timescale of the acoustic cycle is simply $\tau_{\cal R}=2\pi/\Delta\omega_{r}$. Thus the oscillations are well resolved in the eigenspectrum only if the advection time is significantly longer than the acoustic time (i.e. ${\cal M}_{1}>6$).

The values of $\vert{\cal Q}\vert,\vert{\cal R}\vert$ are determined using the following relations:

  

$$\vert{\cal Q}\vert {\cosh \pi{\Delta\omega_{i}\over\Delta\omega_{r}} \over \cosh (n-... ...ga_{i}\over\Delta\omega_{r}}} \exp2n\pi{{\bar\omega}_{i}\over\Delta\omega_{r}},$$ (54)

$$\vert{\cal R}\vert {\sinh n\pi{\Delta\omega_{i}\over\Delta\omega_{r}} \over \cosh (n... ...{i}\over\Delta\omega_{r}}} \exp2\pi{{\bar\omega}_{i}\over\Delta\omega_{r}}\cdot$$ (55)

where

$${{\bar \omega}_{i}\equiv {\omega_{i}^{+}+\omega_{i}^{-}\over2}}$$ (56)

$${\Delta\omega_{i}\equiv \omega_{i}^{+}-\omega_{i}^{-},}$$ (57)

As an illustration, Fig. 8 is obtained from the eigenspectrum of Fig. 4 using Eqs. (54) and (55). The accuracy of the measurements of $\Delta \omega_{r}, \Delta\omega_{i}$ and $\omega_{i}$ is of course of the order 1/n. As a check of consistency, the value of $\vert{\cal Q}\vert$ obtained from the product of $\vert{{\cal Q}}_{\rm adv}\vert$ (Fig. 6) and $\vert{\cal Q}_{\rm sh}\vert$ (Eq. (40)) is also displayed, showing an excellent agreement except for the lowest frequency modes where WKB approximation breaks down. It is remarkable that $\vert{\cal Q}\vert$ seems to be maximized at low frequency for l=1 whereas it is maximized near the cut-off frequency $\omega_l^{\rm cut}$ for $l\ge 2$. The efficiency $\vert{\cal R}\vert$ of the acoustic cycle is bounded by one and decreases to zero above the cut-off frequency, as could be expected from FT00 and F01.

5.2.4 WKB approximation of ${\cal R}$ and growth rate

The global efficiency ${\cal R}$ can be decomposed in terms of the efficiencies ${{\cal R}}_{\rm adv}$ and ${{\cal R}_{\rm sh}}$ of the coupling between the acoustic fluxes $F_{\pm}$, such that the global efficiency ${\cal R}$ is

$${\cal R}\equiv{{\cal R}_{\rm adv}}{{\cal R}_{\rm sh}} {\rm e}^{-i\omega\tau_{\cal R}}.$$ (58)

The efficiency ${{\cal R}}_{\rm adv}$ measures the outgoing acoustic flux F_- produced by the deviation of an ingoing acoustic flux F₊:

$${{\cal R}_{\rm adv}} \equiv {f_{-}\over f_{+}},$$ (59)

$${\cal R}^{*}{\rm e}^{{i\omega\over c}\int_{r_{\rm t}}^{{r}_{\rm sh}} {2\mu\over1-{\cal M}^{2}}{\rm d} r},$$ (60)

$$F_{-}=\vert{{\cal R}_{\rm adv}}\vert^{2}F_{+}.$$ (61)

${{\cal R}}_{\rm adv}$ is simply the limit when $\gamma\to1$ of the efficiency computed in F01: it is close to unity below the cut-off frequency $\omega_l^{\rm cut}$, and decreases exponentially above this cut-off. Conversely, an acoustic flux F_- reaching a shock produces a reflected ingoing acoustic flux F₊ with the efficiency ${{\cal R}_{\rm sh}}$:

 

$${{\cal R}_{\rm sh}}\equiv{f_{+}\over f_{-}},$$ (62)

 

$$F_{+}= \vert{{\cal R}_{\rm sh}}\vert^{2} F_{-}.$$ (63)

Keeping only the first order terms in ${\cal M}_{2}$ and $c/\omega r$, the complex efficiency is computed in Appendix E:

$${{\cal R}_{\rm sh}}\sim-1+2{\cal M}_{2}+ 2i{\eta c\over\omega r}\cdot$$ (64)

Thus the global efficiency $\vert{\cal R}\vert$ should be close to unity below the cut-off frequency. The product $\vert{{\cal R}}_{\rm sh}{{\cal R}}_{\rm adv}\vert$ is displayed in Fig. 8 (bottom picture), in good agreement with the value of $\vert{\cal R}\vert$ deduced from Fig. 4 and Eq. (55).

The extremal values $\omega_{i}^{\pm}$ of the growth rate of the vortical-acoustic cycle are computed in Appendix G (Eq. (G.7)) for strong shocks:

$$\omega_{i}^{\pm}\sim{1\over\tau_{\cal Q}}\log {\vert{\cal Q}\vert\over 1\mp\vert\alpha\vert}\cdot$$ (65)

The asymptotic scaling of $\alpha$ is deduced from Eqs. (43), (51) and (64) :

$$\alpha\sim 1-{3\over{\cal M}_{1}}\log\left\lbrack{{\cal M}_{1... ...\omega_l^{\rm cut}}\right)^{{2\over3}(2l-1)}\right\rbrack\cdot$$ (66)

Together with Eq. (36), and despite the roughness of the approximation in Eq. (30) near the cut-off frequency, the leading order of the asymptotic values of $\omega_{i}^{\pm}$ is comparable to:

  

$$\omega_{i}^{-} \sim {1\over2\tau_{\rm adv}} \log {{\cal M}_{1}\over4{\cal M}_{0l}},$$ (67)

$$\omega_{i}^{+} \sim {3\over2\tau_{\rm adv}} \log {{\cal M}_{1}\over (9{\cal M}_{0l})^{1\over3}\log^{2\over3}{{\cal M}_{1}\over{\cal M}_{0l}}}\cdot$$ (68)

The values of ${\cal Q},{\cal R}, \tau_{\cal R}/\tau_{\cal Q}$are thus responsible for a dispersion of the growth rate by a factor 3in Eqs. (67)-(68) near the cut-off frequency, depending on the relative phases of the two cycles. This factor 3 is consistent with both Figs. 4 and 5.

  
6 Vortical-acoustic instability at low frequency

The region of efficient coupling between the vorticity and acoustic perturbations lies in the region of pseudosound according to Eq. (33). Thus it seems natural to expect an instability at low frequency as long as the coupling region determined from Eq. (32) lies inside the shock radius, i.e.down to the advection frequency. This case is similar to the hole tone instability (e.g.whistling kettle) studied by Chanaud & Powell (1965), or in the oscillations of impinging shear layers (see a review by Rockwell 1983). In some of these cycles, the vorticity perturbations are coupled to the acoustic field in the region of pseudosound. As stressed by Chanaud & Powell (1965), this does not preclude the use of the term "acoustic feedback'' so that we may talk of a vortical-acoustic cycle even in this range of frequencies. Although the problem cannot be treated in the WKB limit, the homogeneous solution can be approximated by a Spherical Bessel function of the first kind in the region of pseudosound far from the sonic point. The calculation in Appendix H shows that the case l=1 and $l\ge 2$ must be treated separately. In the domain of pseudosound $1\ll \vert\omega\tau_{\rm adv}\vert\ll{\cal M}_{1}$, Eq. (22) is reduced to:

 

$${i\omega r\over {\cal M}_{2} c}r^{2}{\partial^{2} f_{0}\over\partial r^{2}} \sim -l(l+1)(l+\eta-2)f_{0}$$

$$+9\;l\;{\rm e}^{i\omega \tau_{\rm adv}} \Gamma\left({l+4\over3}\right)(i\omega\tau_{\rm adv})^{5-l\over3}f_{0}.$$ (69)

6.1 Low frequency global cycle l=1

The l=1 acoustic perturbation is approximated in Appendix H as $f_{0}\propto j_{1}(\omega r/c)$, and Eq. (69) is transformed into

$${9\over5}{\cal M}_{2}^{2}\left(i\omega\tau_{\rm adv}\right)^{... ...}\Gamma\left({5\over3}\right) {\rm e}^{i\omega\tau_{\rm adv}}.$$ (70)

A branch of solutions with $\omega_{r}\gg\omega_{i}$corresponds to a vortical-acoustic cycle between the shock and the sonic point:

$$\omega_{i}\tau_{\rm adv}\sim \log\left\lbrack{3\over2}\Gamma\left({5\over3}\right)\right\rbrack +{4\over3}\log \vert\omega_{r}\tau_{\rm adv}\vert ,$$

$${\rm if}\;\;\vert\omega\tau_{\rm adv}\vert\ll{\cal M}_{1}^{2\over3},$$ (71)

$$\omega_{i}\tau_{\rm adv}\sim \log\left\lbrack{5\over9}\Gamma\left... ...ht\rbrack +2\log{\cal M}_{1}-{5\over3}\log \vert\omega_{r}\tau_{\rm adv}\vert ,$$

$${\rm if}\;\;\vert\omega\tau_{\rm adv}\vert\ll{\cal M}_{1}^{2\over3}.$$ (72)

The maximum growth rate of this branch of solution is reached for a real frequency similar to the frequency of the local instability.

  

$$\omega_{r} \sim {10^{1\over3}{\cal M}_{1}^{2\over3} \over3\tau_{\rm adv}},$$ (73)

$$\omega_{i} \sim {8\over 9\tau_{\rm adv}}\log{{\cal M}_{1}\over{\cal M}_{0}},$$ (74)

$${\cal M}_{0}\equiv{2^{5\over8}3^{3\over8}\over5^{1\over2} \Gamma^{9\over8}\left({5\over3}\right)}\sim 1.17.$$ (75)

This growth rate is compared to the results of numerical calculations in Fig. 5. This growth rate exceeds the growth rate of the vortical-acoustic instability at high frequency (Eq. (46)), but is asymptotically smaller than the maximum growth rate $\omega_{i}^{+}$ reached when the purely acoustic and vortical acoustic cycles are in phase (Eq. (68)).

6.2 Low frequency global cycle $l\ge 2$

If $l\ge 2$, the acoustic feedback due to an isolated vortical perturbation near the shock is negligible compared to the integral effect for strong shocks. Eq. (69) is reduced to:

 

$$3\Gamma\left({l+4\over3}\right)\left(i\omega\tau_{\rm adv}\right)^{2-l\over3} {\rm e}^{i\omega\tau_{\rm adv}}=l-1 .$$ (76)

The most unstable modes are therefore the modes l=2, with

  

$$\omega_{r} \sim {2n\pi\over\tau_{\rm adv}}, \;\;{\rm with}\;\;1\ll n \ll{\cal M}_{1},$$ (77)

$$\omega_{i} \sim {\log 3\over\tau_{\rm adv}}>0.$$ (78)

These approximations are in excellent agreement with the numerical result of Fig. 4.

By contrast with the case of high frequency perturbations, it seems difficult to separate completely the effect of postshock acceleration from the global vortical-acoustic cycle in this range of frequencies.

  
7 Effects of post-shock acceleration: A global instability with a local criterion

The boundary condition (18) can be rewritten for strong shocks as follows:

$${f_{\rm sh}\over c^{2}}\sim \eta{\Delta\zeta\over {r}_{\rm sh... ...\eta\vert\over{\cal M}_{2}^{1\over2}}{c\over r_{\rm son}}\cdot$$ (79)

This equation is similar to the argument of Nobuta & Hanawa (1994) concerning the balance total pressures (thermal and dynamical) on both sides of the shock. Rather than treating the shock surface as a material surface pushed by the local pressures on both sides of it, Eq. (79) is interpreted as follows: an excess of Bernoulli perturbation f (which we may call "energy'') on the subsonic side of the shock is associated with a displacement of the shock in the direction of increase of the local Mach number. If $\eta <0$ (Eq. (17)), an excess of f produces an outward displacement of the shock. This statement alone is not conclusive: instability occurs only if the flow is not able to evacuate this excess of energy. Although the displacement of the shock indeed liberates some potential energy locally, the instability depends on the leakage of energy through the sonic radius (Nakayama 1993).

7.1 Asymptotic growth rate of the radial instability

For radial perturbations (l=0), the only way of evacuating the excess of energy is through acoustic perturbations. Acoustic energy, however, is trapped in the subsonic region if the real frequency of the mode is low enough. Thus a low frequency instability is expected. Although the criterion of the instability is indeed local, its growth rate depends on the spatial structure of the acoustic perturbation from the sonic radius to the shock, as illustrated by Eq. (22) for strong shocks:

$${\partial\log f_{0}\over\partial\log r}({r}_{\rm sh})\sim{1\o... ...eft(1+{i\eta c{\cal M}_{2}\over\omega{r}_{\rm sh}}\right)\cdot$$ (80)

f₀ is approximated as the acoustic perturbation of a uniform medium in spherical coordinates (see Appendix B), using a Spherical Bessel function $f_{0}\propto j_{0}(\omega r/ c)$ for ${r}_{\rm sh}\gg r_{\rm son}$. The solution of the dispersion relation (80) at low frequency is a purely imaginary eigenfrequency:

   

$$\omega_{r} \sim$$ 0, (81)

$$\omega_{i} \sim {-\eta\over 3+\eta}{3{\cal M}_{2}c\over{r}_{\rm sh}},$$ (82)

$$\sim -6{\partial \vert v\vert\over\partial r}\sim {2\over\tau_{\rm adv}}\cdot$$ (83)

The fact that the growth time scales like the advection time is not obvious a priori.

7.2 Asymptotic growth rate of the l = 1 instability

Figure 9: Interaction of the vortex with the shock. The contribution of a vorticity perturbation $\delta K_{\rm sh}>0$ to the Bernoulli perturbation is positive ( $v\delta v_{r}>0$), and induces an increase of the local accretion rate (g>0).Depending on the sign of the postshock acceleration, a vorticity perturbation contributes to a local expansion ($\eta >0$) or collapse ($\eta <0$) of the shock.

Open with DEXTER

Non radial perturbations generate vorticity, which contribute to the energy balance in the subsonic part of the flow. Figure 9 illustrates the contribution of a vorticity perturbation to the Bernoulli perturbation through the radial component of velocity $v~\delta v_{\rm r}$. Thus the vorticity perturbation participates to increase the local accretion rate (i.e.g>0) in the regions where the shock moves inward ( $\Delta\zeta<0$). Comparing this statement with $g\sim i\omega \Delta \zeta/{\cal M}_{2}c$ deduced from Eq. (19), we conclude that the vorticity contributes to the instability if $\nu<0$. A quantitative calculation of the growth rate requires taking into account the acoustic perturbation in the region from the shock to the sonic point. The fastest instability of the flow, well isolated from all the other modes of instability, is obtained by neglecting the last term on the right hand side of Eq. (70):

  

$$\omega_{r}^{\rm max} {5^{1\over3}3^{1\over2}\over2^{2\over3}} {c\over{r}_{\rm sh}}{1\over{\cal M}_{1}^{1\over3}},$$ (84)

$$\omega_{i}^{\rm max} {5^{1\over3}\over2^{2\over3}} {c\over{r}_{\rm sh}}{1\over{\cal M}_{1}^{1\over3}}\cdot$$ (85)

It is in excellent agreement with the numerical calculations in Figs. 4 and 5. This growth rate is intermediate between the advection and the fundamental acoustic frequencies. The real and imaginary parts being comparable, the growth rate is the time taken by the vorticity perturbation to travel by one wavelength $2\pi v/\omega_{r}^{\rm max}\propto {\cal M}_{2}^{2\over 3}{r}_{\rm sh}\ll{r}_{\rm sh}$.

  
8 Discussion

8.1 Relationship with other shock instabilities

Most shock instabilities identified and studied in astrophysics are related to the acceleration or deceleration of the shock itself. The most famous is of course the Rayleigh-Taylor instability of accelerated shocks, analysed by Bernstein & Book (1978). Decelerated shocks are also unstable to a rippling instability, studied by Vishniac (1983), Bertshinger (1986) and Vishniac & Ryu (1989). By definition, stationary shocks are stable with respect to these mechanisms, but can still be unstable. Nakayama (1992, 1993) pointed out the radial instability of the shock if the flow is immediately accelerated after the shock, in isothermal flows. The validity of the postshock acceleration criterion for adiabatic flows is still uncertain, even for radial perturbations (Nakayama 1994). The vortical-acoustic cycle studied in the present paper resembles by many aspects to the entropic-acoustic cycle studied by FT00 and F01 in adiabatic flows. Two distinct mechanisms, however, are at work:

(i) the vortical-acoustic cycle is fed by the vorticity production by a perturbed shock, which is highest for isothermal flows with a strong shock,

(ii) the entropic-acoustic cycle is fed by the temperature increase from the shock to the sonic radius, which is highest if $\gamma=5/3$(FT00, F01).

A calculation similar to Appendix E for adiabatic flows would show that the vorticity production by a perturbed shock is large only in the isothermal limit:

$$\left\vert{\delta K_{\rm sh}\over f_{-}}\right\vert\propto {1... ...l M}_{2}}< \left({2\gamma\over\gamma-1}\right)^{1\over2}\cdot$$ (86)

Although shocked adiabatic flows with $1<\gamma<5/3$ are subject to both entropic-acoustic and vortical-acoustic cycles in principle, their stability cannot be determined without a specific calculation.

8.2 Constraints on the instability mechanism of BHL accretion

The stability of BHL accretion can be addressed with new tools, using the present results and the entropic-acoustic cycle of FT00 and F01. Although a detailed analysis of the existing numerical simulations of BHL accretion is postponed to a forthcoming paper, let us outline the possible consequences of the entropic-acoustic and vortical-acoustic mechanisms on this specific flow. Numerical simulations of BHL accretion in 3-D show a strong instability for adiabatic flows with $\gamma=5/3$ and small accretors (Ruffert & Arnett 1994; Ruffert 1994): this coincides nicely with the properties of the entropic-acoustic cycle, which is most unstable if the temperature gradient in the subsonic part of the flow is strongest (FT00, F01). By contrast, nearly isothermal flows are rather stable in 3-D numerical simulations (Ruffert 1996). This could seem puzzling given the strong vortical-acoustic instability described in the present paper: why would the entropic-acoustic cycle be relevant for the BHL instability with $\gamma=5/3$, and the vortical-acoustic cycle be irrelevant to the BHL stability for $\gamma=1$? This apparent contradiction can be partly solved by remembering the difference of topology between the spherical Bondi flow and the BHL flow, following the remark of Sect. 2.1. All the numerical simulations of isothermal BHL accretion seem to agree with the fact that the shock is attached to the accretor, whereas it is detached in high resolution simulations of adiabatic flows with $\gamma\ge4/3$. This can be understood qualitatively in terms of the weakness of pressure forces in isothermal flows compared to adiabatic flows ( $P\propto\rho^\gamma$). More quantitatively, Foglizzo & Ruffert (1997) proved that the shock in the isothermal BHL flow cannot be detached unless the sonic surface extends up to distances comparable to the Bondi radius GM/c², i.e.much larger than the accretion radius. In view of this topological difference between isothermal and adiabatic BHL flows, it seems easier to extend the results obtained for shocked Bondi accretion to the subsonic region ahead of the accretor, for detached shocks, than in the very non-radial region of accretion behind the accretor for attached shocks. Although true, however, this argument is not conclusive. Even in a non radial isothermal flow, a small pressure perturbation of the shock is able to generate vorticity perturbations very efficiently. The lack of instability of isothermal BHL accretion in 3-D must be sought for in the lack of acoustic feedback from the advected vorticity perturbation. Unfortunately we do not have quantitative arguments to explain why the acoustic feedback is so weak, apart from noticing that geometric compression of the vorticity perturbation from the shock to the sonic radius might be insufficient in the BHL flow. Any coherent explanation should of course also account for the strong instability observed in 2-D simulations of BHL accretion (Shima et al.1998).

8.3 Relationship with the numerical simulations of shocked inviscid flows with low angular momentum

Inviscid accretion flows with low angular momentum are much more complex than the Bondi flow in the sense that the shock position is not unique (Fukue 1987; Chakrabarti 1989). The numerical simulations of Nobuta & Hanawa (1994) did not detect the effects of vorticity perturbations since they were restricted to axisymmetric motion. These simulations showed the evolution of the radial instability due to postshock acceleration: the unstable shock is either absorbed by the accretor or inflates to reach the stable outer position.

Molteni et al. (1999) performed non axisymmetric numerical simulations of shocked inviscid accretion flows with low angular momentum, for $\gamma=4/3$. Unstable m=1 oscillations were observed, with no physical explanation. The entropic-acoustic and vortical-acoustic cycles might provide a physical basis to understand this instability.

8.4 Non linear evolution of the instability

If the existence of the shock is postulated a priori, as in the present study, the flow might simply converge towards another solution, as illustrated by the dashed lines in Fig. 1. The shock can either be absorbed by the accretor, leading to the fully supersonic solution, or could expand towards infinity to establish the Bondi solution if the outer boundary condition allows it.

In accretion flows stable with respect to the postshock acceleration criterion, but unstable through the vortical-acoustic (or entropic-acoustic) cycle, the instability might be saturated by the effect of the geometric dilution of the acoustic energy in the subsonic cavity. Let us assume that advected vorticity perturbations generate an acoustic flux propagating outward $F_{-}\propto \vert{\cal Q}_{\rm adv}\vert^{2}\vert\delta K_{\rm sh}\vert^{2}$ (Eq. (26)), and that the acoustic perturbation produces in turn a vorticity perturbation $\vert\delta K_{\rm sh}\vert\propto {\cal M}_{1}\vert\delta p_{-}/p\vert$ (Eq. (42)). The amplitude of the pressure perturbation $\vert\delta p_{-}/p\vert$ depends on the volume in which the acoustic flux F_- is diluted through Eq. (27): $\vert\delta p_{-}/p\vert\propto ({r}_{\rm sh}/r) F_{-}^{1/2}$. Thus the vortical-acoustic cycle is naturally stabilized when the shock reaches a distance $r_{\rm max}$ defined by

$$r_{\rm max}\sim \vert{\cal Q}\vert{r}_{\rm sh},$$ (87)

where $\vert{\cal Q}\vert$ is the global efficiency of the vortical-acoustic cycle in the linear regime. If the shock is not simply absorbed by the accretor, the non linear evolution of the instability could lead to quasi periodic oscillations of amplitude comparable to $r_{\rm max}$. However, given the number of unstable modes in a spectrum like Fig. 4, the vortical-acoustic instability might as well saturate into turbulence rather than be dominated by a single QPO. This issue can only be solved with numerical simulations.

   
9 Conclusions

The linear stability of shocked isothermal Bondi accretion has been studied by comparing the complex eigenfrequencies obtained through a direct numerical integration (Sect. 3) to the analytical results obtained for strong shocks by two different methods:

(i) an analytical estimate of the growth rate corresponding to a cycle of perturbations with a purely real frequency, obtained by separating the effects of advection from the boundary effects of the shock. This WKB approximation, valid in the range of acoustic waves ( $c/{r}_{\rm sh}\ll\omega\ll c/r_{\rm son}$), was used in Sects. 4 and 5,

(ii) an analytical estimate of the complex eigenfrequencies in the range of pseudosound ( $v_{\rm sh}/{r}_{\rm sh}\ll\omega\ll c/{r}_{\rm sh}$) using Spherical Bessel functions (Sects. 6 and 7).
The results obtained by these methods are summarized as follows:

- As expected by the postshock acceleration of Nakayama (1993), the isothermal Bondi accretion with a shock is unstable with respect to radial perturbations. Its growth rate is comparable to the advection time from the shock to the sonic point.

- The present analysis has revealed the existence of a new instability, based on the cycle of vortical and acoustic perturbations in the subsonic part of the flow. The analytical study of this instability at high frequency ( $c/{r}_{\rm sh}\ll\omega\ll c/r_{\rm son}$) proves that it is independent of the postshock acceleration criterion established by Nakayama (1992, 1993). It is fed by the efficient production of vorticity perturbations when the shock is perturbed non radially (Eq. (42)), and by the vortical-acoustic coupling in the region of the sonic radius which enables the acoustic feedback. In this sense this non radial instability is generic and can be expected in more complex situations such as shocked flows with a weak angular momentum accreting into a black hole, even if the flow is decelerated immediately after the shock.

In the shocked Bondi flow the vortical-acoustic instability is faster than the radial instability if the shock is strong, by a factor $\propto \log({\cal M}_{1})$.

The vortical-acoustic instability occurs for low degree perturbations on a wide range of frequencies below the cut-off frequency $\sim$ $c/r_{\rm son}$. A branch of unstable l=1 eigenmodes corresponds to a vortical-acoustic cycle in which the acoustic feedback is produced in the pseudosound domain ( $v_{\rm sh}/{r}_{\rm sh}<\omega<c/{r}_{\rm sh}$). The resulting growth rate is comparable to the growth rate of the vortical-acoustic cycle at high frequency.

The role of the purely acoustic cycle was pointed out at high frequency in order to explain the large variations of the growth rate from one mode to another. More generally, the formalism developped in Appendix G concerning the simultaneous acoustic and vortical-acoustic cycles applies to any context where the efficiencies ${\cal Q},{\cal R}$ and timescales $\tau_{\cal Q},\tau_{\cal R}$ can be defined.

- A strong l=1 oscillatory instability was found in the pseudosound domain, at a frequency which is intermediate between the acoustic and advection frequencies ( $\omega_{i}^{\rm max}\propto {\cal M}_{1}^{1/3}c/{r}_{\rm sh}$). With comparable real and imaginary parts of the eigenfrequency, vorticity perturbations are advected over a very short distance during one growth time. On the basis of the contribution of a vortex to the Bernoulli perturbation sketched in Fig. 9, this strong instability is a non radial consequence of post-shock acceleration.

- As outlined in Sect. 8, the vortical-acoustic mechanism can be used as a tool in order to analyse the instability observed in numerical simulations of more complicated accretion flows. The specific application to BHL accretion or shocked flows with low angular momentum will be developped in a forthcoming paper.

Appendix A: Vorticity perturbations produced by the perturbed shock

The non radial perturbation of velocity is deduced from the continuity of the velocity parallel to the shock, as in Landau & Lifshitz (1987, Chap. 90, p. 336):

  

$$\delta v_{\theta} {v_{1}-v_{2}\over r}{\partial\Delta\zeta\over\partial\theta},$$ (A.1)

$$\delta v_{\varphi} {v_{1}-v_{2}\over r\sin\theta} {\partial\Delta\zeta\over\partial\varphi}\cdot$$ (A.2)

The perturbed vorticity in the flow is deduced from the Eqs. (A.1)-(A.2) and the Euler equation:

  

w_r=0, (A.3)

$$w_{\theta} (1-{\cal M}_{2}^{2})^{2} \left(1-{i\eta c{\cal M}_{2}\over\omega{... ...ight){{\cal M}_{1}^{2}\over r\sin\theta}{\partial\Delta v\over\partial\varphi},$$ (A.4)

$$w_{\varphi} -(1-{\cal M}_{2}^{2})^{2} \left(1-{i\eta c{\cal M}_{2}\over\omega... ...sh}}\right){{\cal M}_{1}^{2}\over r} {\partial\Delta v\over\partial\theta}\cdot$$ (A.5)

The conserved quantity $\delta K_{\rm sh}$ defined by Eq. (12) is deduced from Eqs. (A.4)-(A.5), resulting in Eq. (20).

Appendix B: Approximations of the acoustic perturbation

B.0.1 High frequency limit: WKB approximation

If the wavelength of the perturbation is shorter than the lengthscale of the flow inhomogeneity ( $\omega r/c\gg1$), a WKB approximation enables us to describe the propagation of acoustic waves in the direction of the flow (f₀⁺) or against it (f₀⁺):

$$f_{0}^{\pm}\sim {{\cal M}^{1\over2}c^{2}\over \mu^{1\over2}}\... ...m t}}^r {{\cal M}\mp\mu\over1-{\cal M}^{2}}{{\rm d} r\over c},$$ (B.1)

where we have chosen the normalization of $f_{0}^{\pm}$ such that the lower bound of the integral is the turning point $r_{\rm t}$. The condition that $r_{\rm t}<{r}_{\rm sh}$ defines a minimum frequency $\omega_{l\ge1}^{\rm sh}$ using Eq. (11) and $L^{2}\equiv l(l+1)$:

$$\omega_{l\ge1}^{\rm sh}\equiv \omega_{l}({r}_{\rm sh})=L{c_{... ... {r}_{\rm sh}} \left(1-{\cal M}_{\rm sh}^{2}\right)^{1\over2}.$$ (B.2)

B.0.2 Low frequency limit: Uniform steady medium approximation

Far from the accretor, the flow velocity decreases like $\propto 1/r^{2}$ and the density of the gas is uniform. If the shock is strong ( ${r}_{\rm sh}\gg r_{\rm son}$), the homogeneous equation associated with Eq. (C.1) is approximated for $r\gg r_{\rm son}$ and $\omega\ll c/(r{\cal M})$ as follows:

$${\partial^{2} f\over\partial r^{2}}+{2\over r}{\partial f\ove... ... +\left({\omega^{2}\over c^{2}}-{L^{2}\over r^{2}}\right) f=0,$$ (B.3)

which is nothing more than the equation of acoustic waves in a uniform steady medium, in spherical coordinates. The solution f₀ can therefore be approximated with a spherical Bessel function of the first kind j_l, which is normalized here as in Eq. (B.1):

 

$$\sim {\rm e}^{3\over 4} c^{2}{\omega r_{\rm son}\over c} \left({\pi c\over 2\omega r}\right)^{1\over 2} J_{l+{1\over2}}\left({\omega r\over c}\right),$$ (B.4)

$$\equiv {\rm e}^{3\over 4}c^{2}{\omega r_{\rm son}\over c} j_{l}\left({\omega r\over c}\right).$$ (B.5)

Appendix C: Formulation of the boundary value problem

The differential system (8-9) is written as a single differential equation of second order:

 

$${\partial^{2} f\over\partial r^{2}}+\left\lbrack {\partial\log\ov... ...\omega {\cal M}\over (1-{\cal M}^{2})c}\right\rbrack{\partial f\over\partial r}$$

$$+{{\omega^{2}\over c^{2}}-{L^{2}\over r^{2}}\over 1-{\cal M}^{2}}... ...1-{\cal M}^{2})r^{2}} {\rm e}^{i\omega\int_{{r}_{\rm sh}}^r{{\rm d} r\over v}}.$$ (C.1)

Following the method used in F01, the general solution of Eq. (C.1) can be written using the solutions f₀,f₁ of the homogeneous equation, where f₀ is the unique solution which is regular at the sonic radius. Let us normalize f₀ using the following definition of the WKB solutions $f_{0}^{\pm}$ (Eq. (B.1)). The complex coefficient ${\cal R}^{*}$ is defined by:

$$f_{0}\equiv f_{0}^{+}+{\cal R}^{*}f_{0}^{-}.$$ (C.2)

The wronskien of the couple of solutions f₀⁺,f₀^- is deduced from Eqs. (B.1) and (C.1):

$${\partial f_{0}^{+}\over\partial r}f_{0}^{-}-{\partial f_{0}^... ...c}\int_{r_{\rm t}}^r{{\cal M}\over 1-{\cal M}^{2}}{\rm d} r}.$$ (C.3)

Let us normalize f₁ such that the wronskien of the couple of solutions f₀,f₁ is the same as in Eq. (C.3). The general solution is then:

 

$${\delta K_{\rm sh}\over 2i\omega c^{3}} {\rm e}^{-i\omega\int_{r_... ...{2}\over 1-{\cal M}^{2}}{{\rm d} r\over{\cal M}} }{\rm d} r\right\rbrack\right.$$

$$\qquad-\left.f_{1}\left\lbrack A_{1}+\int_{r_{\rm son}}^r {f_{0}\... ...{\cal M}^{2}}{{\rm d} r\over{\cal M}} }{\rm d} r\right\rbrack\right\rbrace\cdot$$ (C.4)

A Frobenius analysis of f₀,f₁ near the sonic points leads to:

$$\sim f_{0}({r}_{\rm sh}) + {\cal O}(r-r_{\rm son}),$$ (C.5)

$$\propto (r-r_{\rm son})^{-{i\omega\over {\dot{\cal M}}c}}.$$ (C.6)

thus the integrals in Eq. (C.4) are converging when $r\to r_{\rm son}$if $\omega_{i}>c(\partial{\cal M}/\partial r)(r_{\rm son})$. The regularity at the sonic radius therefore requires A₁=0. A combination of Eq. (C.4) and its derivative at the shock radius leads to eliminate A₀ as follows:

 

$$\left({\partial f_{0}\over\partial r}f-{\partial f\over\partial r... ... sh}}^r{1+{\cal M}^{2}\over 1-{\cal M}^{2}}{{\rm d} r\over{\cal M}}} {\rm d} r.$$ (C.7)

Equations (18)-(20) provide the boundary conditions at the shock in order to replace $f({r}_{\rm sh})$, $\partial f/\partial r ({r}_{\rm sh})$ and $\delta K_{\rm sh}$ in Eq. (C.7), resulting in Eq. (22).

Appendix D: Calculation of ${\cal Q}_\mathsf{adv}$

The analytical expression for ${\cal Q}_{\rm adv}$ can be determined by writting the solution corresponding to zero acoustic flux F₊ at an outer boundary R, using the couple of homogeneous solutions ( f₀⁺,f₀^-):

 

$${\delta K_{R}\over 2i\omega c^{3}} {\rm e}^{-i\omega\int_{r_{\rm ... ...{2}\over 1-{\cal M}^{2}}{{\rm d} r\over{\cal M}}} {\rm d} r\right\rbrack\right.$$

$$-\left.f_{0}^{-}\left\lbrack B_{1}+\int_{r_{\rm son}}^R {f_{0}^{+... ...{\cal M}^{2}}{{\rm d} r\over{\cal M}} }{\rm d} r\right\rbrack\right\rbrace\cdot$$ (D.1)

The regularity at the sonic radius requires

$$B_{1}=-{\cal R}^{*}B_{0}.$$ (D.2)

The condition of absence of an incoming acoustic flux at the outer boundary is obtained by canceling the coefficient of f₀⁺ at r=R, in the WKB limit of high frequency ( $\omega\gg c/R$):

$$B_{0}=-\int_{r_{\rm son}}^{R} {f_{0}^{-}\over{\cal M} r^{2}} ... ...{1+{\cal M}^{2}\over 1-{\cal M}^{2}}{{\rm d} r\over{\cal M}}}.$$ (D.3)

Together with Eqs. (D.2) and (D.3), Eq. (D.1) at the outer boundary becomes:

 

$$-{\delta K_{R}\over 2i\omega c^{3}} {\rm e}^{{i\omega\over c}\int_{r_{\rm t}}^{R}{{\rm d} r\over {\cal M}}}$$

$$f_{0}^{-}\int_{r_{\rm son}}^R {f_{0}\over{\cal M} r^{2}} {\rm e}^... ...m t}}^r{1+{\cal M}^{2}\over 1-{\cal M}^{2}}{{\rm d} r\over{\cal M}} }{\rm d} r.$$ (D.4)

${{\cal Q}}_{\rm adv}(R)$, defined by Eq. (28), is deduced from the asymptotic behaviour of f₀^- in the WKB approximation (Eq. (B.1)). This calculation is formally similar to the calculation of ${\cal Q}_{K}$ in F01, corrected for a phase shift. The efficiency ${\cal Q}_{\rm adv}$ involved in the vortical-acoustic cycle is deduced from Eq. (29), with $R={r}_{\rm sh}$. In the strong shock limit, below the cut-off frequency ( $c/r_{\rm son}\gg\omega\gg c/{r}_{\rm sh}$), the acoustic efficiency ${\cal Q}_{\rm adv}$ is approximated using the Spherical Bessel function j_l (B.5):

 

$$\vert{\cal Q}_{\rm adv}\vert \sim {{\rm e}^{-{3\over4}}\over2}{c\over \omega r_{\rm son}} \left\vert\int_{0}^\infty {\rm e}^{-i\lambda x^{3}}j_{l}(x){\rm d} x\right\vert,$$ (D.5)

$$\lambda \equiv {{\rm e}^{-{3\over2}}\over3}\left({c\over\omega r_{\rm son}}\right)^{2}\cdot$$ (D.6)

For $\lambda\gg1$, a function ${\cal H}(l)$ may be defined such that

$$\left\vert\int_{0}^\infty {\rm e}^{-i\lambda x^{3}}j_{l}(x){\... ...ot3\ldots(2l+1)} {{\cal H}(l)\over\lambda^{{l+1\over 3}}}\cdot$$ (D.7)

The main contribution to the integral in Eq. (D.7) comes from the region $x\sim \lambda^{-1/3}$. Equation (30) is obtained from Eqs. (D.5) and (D.7)

Appendix E: Decomposition of the perturbation onto vorticity waves and acoustic waves in the WKB approximation

The perturbations f,g immediately after the shock, defined by Eqs. (18) and (19) are decomposed as follows:

  

$$f({r}_{\rm sh})$$ = f_-+f₊+f_K, (E.1)

$$g({r}_{\rm sh})$$ = g_-+g₊+g_K, (E.2)

where f_K,g_K correspond to the vorticity wave associated with the vorticity perturbation $\delta K$, and $f_{\pm},g_{\pm}$ correspond to the purely acoustic waves propagating in the direction of the flow (index +) or against the flow (index -). An exact calculation can be made in the case of the reflexion of an acoustic wave $\delta p_{-}$ with wavevector ( $k_{\parallel},k_{\perp}$) on a plane shock in Cartesian coordinates, in the absence of a gradient of ${\cal M}$. The vorticity wave f_K,g_K is advected at the velocity of the fluid:

$${\partial f_K\over\partial r} {i\omega \over v}f_K,$$ = (E.3)

$${\partial g_K\over\partial r} {i\omega \over v}g_K.$$ = (E.4)

Replacing these derivatives in Eqs. (8), (9), we obtain:

  

$${{\cal M}_{2}^{2}c^{2}\delta K\over r^{2}\omega^{2}+v^{2}L^{2}},$$ (E.5)

$${\delta K\over r^{2}\omega^{2}+v^{2}L^{2}}$$ (E.6)

$$g_\pm=\pm{f_{\pm}\over{\cal M}_{2} c^{2}}\cdot$$ (E.7)

Equations (E.5)-(E.7) are used with Eqs. (18), (19) in order to obtain Eq. (20) and

 

$$f_{\pm} {{\cal M}_{2} c^{2}\over2\mu}{\Delta v\over v}(1-{\cal M}_{2}^{2}) {(\mu\mp{\cal M}_{2})^{2}\over1\mp\mu{\cal M}_{2}},$$ (E.8)

where $\mu$ in Cartesian coordinates is also defined by Eq. (10), but replacing l(l+1)/r² by the wavenumber $k_{\perp}^{2}$. According to the definitions of ${\cal Q}_{\rm sh}$ and ${{\cal R}_{\rm sh}}$ in Eqs. (37) and (62), together with Eqs. (20) and (E.8):

  

$${{\cal R}_{\rm sh}} -\left({\mu-{\cal M}_{2}\over\mu+{\cal M}_{2}}\right)^2 \left({1+\mu{\cal M}_{2}\over1-\mu{\cal M}_{2}}\right),$$ (E.9)

$${{\cal Q}}_{\rm sh} -2l(l+1)\left({\mu\over{\cal M}_{2}}\right)^{1\over2} {(1+\mu{\cal M}_{2})(1-{\cal M}_{2}^{2})\over(\mu+{\cal M}_{2})^{2}}\cdot$$ (E.10)

These equations show the decrease of the vortical-acoustic coupling for weak shocks ( ${\cal M}_{2}\sim1$), and the existence of maximal efficiency at low frequency. Indeed, the maximum ${{\cal Q}}_{\rm sh}\sim {3^{3\over2}\over8}l(l+1){\cal M}_{1}^{2}$ is reached for a frequency such that $\mu\sim {\cal M}_{2}/3$, with ${{\cal R}_{\rm sh}}\sim -1/ 4$.

The angle $\psi$ between the direction of propagation of the wave and the vector orthogonal to the shock surface is given by:

$$\tan\psi = \left\lbrack\left({\omega \over k_{\perp} c}\right)^{2}+{\cal M}^{2}-1 \right\rbrack^{-1}\cdot$$ (E.11)

As remarked by Kontorovich (1958), the reflected sound wave propagates away from the shock with the same angle as the incident wave ( $\psi_+=\psi_-$), as in a classical reflexion. The presence of a gradient of ${\cal M}_{2}$ in the Bondi flow precludes the use of these formulae at low frequency. In the following calculation we assume $\omega r/c\gg1$ and ${\cal M}_{2}\ll1$, and keep only the first order terms in ${\cal M}_{2}$ and $c/\omega r$. Thus $\mu\sim 1$. Neglecting the coupling between the vorticity and acoustic waves in the vicinity of the shock, the vorticity wave f_K,g_K is still advected at the velocity of the fluid. Eqs. (E.5) and (E.6) are now approximated at high frequency by

  

$$\sim {{\cal M}_{2}^{2}c^{2}\delta K\over r^{2}\omega^{2}},$$ f_K (E.12)

$$\sim {\delta K\over r^{2}\omega^{2}}\cdot$$ g_K (E.13)

Acoustic waves are described by Eqs. (8), (9) in the absence of vorticity perturbations, i.e.when $\delta K=0$. Using the WKB approximation of Eq. (B.1), the radial derivative of $f_{\pm}$ is approximated by:

$${\partial f_\pm\over\partial r}\sim{i\omega\over c} \left(\mp 1+{\cal M}_{2}-{i\eta c\over2\omega r}\right) f_\pm.$$ (E.14)

$g_{\pm}$ is deduced from Eqs. (8) and (E.14):

$$g_\pm\sim\left(\pm1 +{i\eta c\over2\omega r}\right) {f_{\pm}\over{\cal M}_{2} c^{2}}\cdot$$ (E.15)

Equations (E.12), (E.13) and (E.15) are used with Eqs. (18), (19) in order to obtain Eq. (20) and

 

$$f_{\pm} \sim {{\cal M}_{2} c^{2}\over2}{\Delta v\over v} \left(\pm1-{\cal M}_{2}- i{\eta c\over\omega r}\right),$$ (E.16)

Combining Eq. (20) with Eq. (E.16), we obtain the expressions for ${{\cal R}}_{\rm sh},{{\cal Q}}_{\rm sh}$ in Eqs. (40) and (64). The pressure perturbations  $\delta p_{\pm}$ associated with the acoustic waves $f_{\pm}$ are deduced from Eqs. (6)-(7) and Eq. (E.15):

 

$${\delta p_{\pm}\over p} {f_{\pm}-v^{2}g_{\pm}\over c^{2}-v^{2}},$$ = (E.17)

$$\sim (1\mp{\cal M}_{2}){f_{\pm}\over c^{2}}\cdot$$ (E.18)

Appendix F: WKB approximation of the dispersion relation

Using the integral expression of ${\cal Q}_{\rm adv}$ (Eq. (29)), Eq. (22) can be approximated as follows for strong shocks ( ${\cal M}_{2}\ll1$):

 

$$\left\lbrack \eta - {\cal M}_{2}{i\omega {r}_{\rm sh}\over c}\rig... ...ft\lbrack {i\omega {r}_{\rm sh}\over c}-{\cal M}_{2}\eta\right\rbrack f_{0}\sim$$

$$-{2l(l+1)\over{\cal M}_{2}^{1\over2}} \left\lbrack{i\omega {r}_{\rm sh}\over c}+{\cal M}_{2}\eta\right\rbrack {\cal Q}_{\rm adv}f_{0}^+ .$$ (F.1)

In the WKB limit $\omega {r}_{\rm sh}/c \gg 1$, Eq. (F.1) can be simplified into

 

$$f_{0}^+ +{\cal R}^{*}f_{0}^-\sim-2l(l+1){\cal Q}_{\rm adv}{f_{0}^+\over{\cal M}_{2}^{1\over2}}\cdot$$ (F.2)

Using the expressions for ${\cal Q}_{\rm sh}$ and ${{\cal R}_{\rm sh}}$ (Eqs. (40) and (64)), and the definition of ${\cal Q}$ and ${\cal R}$ (Eqs. (38) and (58)), the global dispersion relation (49) is recovered.

Appendix G: Analysis of the phase relation of the two cycles

Let us introduce the new complex variable

$$z\equiv {\cal Q}{\rm e}^{i\omega\tau_{\cal Q}}.$$ (G.1)

The resolution of the dispersion relation (49) is then equivalent to finding the complex number z satisfying:

$$z+\alpha z^\epsilon=1,$$ (G.2)

where the dimensionless parameter $\epsilon\equiv \tau_{\cal R}/\tau_{\cal{Q}}<1$ and $\vert\alpha\vert<1$ is defined by Eq. (52). The minimum and maximum values $z_{\pm}$ of |z| satisfy an equation similar to Eq. (E3) of FT00:

$$z_{\pm}\mp \vert\alpha\vert z_{\pm}^\epsilon=1.$$ (G.3)

The minimum and maximum effect of the acoustic cycle on the growth rate $\omega_{i}^{\pm}$ are directly related to $z_{\pm}$ through Eq. (G.1):

$$\omega_{i}^{\pm}\equiv {1\over\tau_{\cal Q}}\log{\vert{\cal Q}\vert\over z_{\mp}}\cdot$$ (G.4)

From Eq. (G.3),

z₊+z_-=2, (G.5)

we deduce that the maximum stabilizing effect of the acoustic cycle is to divide the efficiency ${\cal Q}$ by a factor 2:

$$1< z_{+}\le 2.$$ (G.6)

If the acoustic time is much shorter than the advection time ( $\epsilon\ll1$), we deduce from Eq. (G.3) the values of $x_{\pm}$:

 

$$z_{\pm}\sim 1 \pm \vert\alpha\vert.$$ (G.7)

Thus the acoustic cycle may participate efficiently to the instability if $\epsilon\ll1$ and $\vert\alpha\vert\sim 1$. In the more general case where $0<\epsilon<1$, the following bounds on x are obtained from Eq. (G.2):

$$1-\vert\alpha\vert\le \vert z\vert\le{1\over1-\vert\alpha\vert}\cdot$$ (G.8)

The resolution of Eq. (G.2) can be decomposed into a phase condition applied to the points of a continuous curve ${\cal C}$, defined by:

$$\vert z-1\vert=\vert\alpha \vert \vert z\vert^\epsilon.$$ (G.9)

${\cal C}$ is a closed curve containing the point (1,0) in its interior, and (0,0) in its exterior. It can be described in a univoque way by the angle $\varphi\equiv {\rm Arg}(z-1)\in[-\pi,\pi]$. Let us define the angle $\theta\equiv {\rm Arg}(z)\in[-\pi,\pi]$. The solutions of Eq. (G.2) are recovered by applying to the solutions of Eq. (G.9) the following phase condition:

$$\varphi+2k\pi=\epsilon(\theta+2k'\pi),$$ (G.10)

where k,k' are two integers. Since (0,0) is exterior to ${\cal C}$, the range of values covered by $\theta$ when $\varphi$ covers $[-\pi,\pi]$ is limited to $\vert\theta\vert\le\theta_{\rm max}<\pi$. The discrete solutions of the phase Eq. (G.10) can be seen in a graphic way as the intersection of the periodic curve $\theta(\varphi)$ with the straight line $\theta=\varphi/\epsilon$. The number of solutions in each phase $\varphi\in[0,2\pi]$ is therefore equal to $1/\epsilon$ on average. Comparing Eq. (G.1), Eq. (G.9) with the definitions of $\varphi,\theta$, one period of $\varphi$ corresponds to a variation of $\omega_{r}$ of $2\pi/\tau_{\cal R}$, and one period of $\theta$corresponds to $2\pi/\tau_{\cal Q}$. The values of $\vert{\cal Q}\vert$ and $\vert{\cal R}\vert$ in Eqs. (54) and (55) can be determined from the measurement of $\omega_{i}^{+}$, $\omega_{i}^-$ and the average number $n=\tau_{\cal Q}/\tau_{\cal R}$ of eigenmodes per period $\Delta\omega_{r}=2\pi/\tau_{\cal R}$, by eliminating $\vert\alpha\vert,x_{\pm}$from the set of Eqs. (G.3), (G.4), and (52).

Appendix H: Pseudosound approximation of the dispersion relation

The integral involved in Eq. (22) can be approximated in the low frequency limit by introducing the complex variable z:

$${{\rm d} z\over{\rm d} r}\equiv i{\omega\over {\cal M} c}{1+{\cal M}^{2}\over1-{\cal M}^{2}}\cdot$$ (H.1)

The contour of integration is deformed in the complex plane introducing a point on the real axis at $+\infty$, and performing two integrations by parts:

 

$${{r}_{\rm sh}\int_{r_{\rm son}}^{{r}_{\rm sh}}{f_{0}\over{\cal M}... ... d} r\over{\cal M}}}{\rm d} r={{r}_{\rm sh}c\over i\omega}{\rm e}^{z_{\rm sh}}}$$

$${\times\left\lbrace \int_{z_{\rm son}}^{+\infty} {f_{0}\over r^{2... ...0}\over r^{2}}{1-{\cal M}^{2}\over1+{\cal M}^{2}}\right){\rm d} z\right\rbrace}$$

$${-{c\over i\omega {r}_{\rm sh}}{1-{\cal M}^{2}\over1+{\cal M}^{2}... ...2}}{1-{\cal M}^{2}\over1+{\cal M}^{2}}\right)\right\rbrace_{{r}_{\rm sh}}\cdot}$$ (H.2)

The turning point in the Bondi flow would be at $r_{\rm t}\sim L c/\omega$. For $r\gg r_{\rm son}$, the homogeneous solution f₀ is approximated by a Spherical Bessel function of the first kind $j_{l}(\omega r/c)$. For $\omega r/c\ll 1$, it is approximated as follows:

$$f_{0}(r)\sim f_{0}({r}_{\rm sh}) \left({r\over {r}_{\rm sh}}\right)^l\cdot$$ (H.3)

The first integral on the right hand side of Eq. (H.2) is approximated using a Gamma function:

 

$${z\sim {i\omega r\over 3{\cal M} c},}$$ (H.4)

$${{{r}_{\rm sh}c\over i\omega} \int_{z_{\rm son}}^{+\infty} {f_{0}... ... z_{\rm sh}^{l+1\over3}} \int_{0}^{+\infty}z^{l-2\over3}{\rm e}^{-z}{\rm d} z,}$$ (H.5)

$${\sim{f_{0}({r}_{\rm sh})\over 3{\cal M} z_{\rm sh}^{l+1\over3}} \Gamma\left({l+1\over3}\right)\cdot}$$ (H.6)

The second integral on the right hand side of Eq. (H.2) is negligible if $r_{\rm t}\gg{r}_{\rm sh}\gg r_{\rm son}$. In view of the particular case l=1, the differential equation satisfied by f₀ is used to sum up terms of same order:

 

$$\eta r{\partial f_{0}\over\partial r} + l(l+1)f_{0} = r^{2}(1-{\c... ...{\partial^{2} f_{0}\over\partial r^{2}}+\left({\omega r\over c}\right)^{2}f_{0}$$

$$-{\cal M}^{2}\left(\eta+{2i\omega r\over c{\cal M}}\right) r{\partial f_{0}\over\partial r}\cdot$$ (H.7)

Using Eqs. (H.2), (H.6) and (H.7), the dispersion relation (22) is approximated in the pseudosound domain as follows:

 

$$r^{2}{\partial^{2}f_{0}\over\partial r^{2}} + {l(l+1)\over 3z_{\r... ...^{z_{\rm sh}}\over z_{\rm sh}^{l-2\over3}} \Gamma\left({l+4\over3}\right) \cdot$$ (H.8)

The left hand side of Eq. (H.8) is dominated by the second derivative of f₀ if $l\ge 2$:

 

$${l-1\over3}z_{\rm sh}^{l-2\over3}{\rm e}^{-z_{\rm sh}} \sim \Gamma\left({l+4\over3}\right)\cdot$$ (H.9)

If l=1, the second derivative of f₀ is approximated as follows:

$${j_{1}(x)\sim{x\over3},}$$ (H.10)

$${x^{2}{\partial^{2}j_{1}\over\partial x^{2}}\sim -{x^{3}\over 5} \cdot}$$ (H.11)

Thus Eq. (H.8) becomes:

 

$${9\over 5} {\cal M}^{2}z_{\rm sh}^{3} + {2\over 9} \left\lbrack 1... ...sim z_{\rm sh}^{4\over3} {\rm e}^{z_{\rm sh}} \Gamma\left({5\over3}\right)\cdot$$ (H.12)

References

 
• Bernstein, I. B., & Book, D. L. 1978, ApJ, 225, 633 In the text NASA ADS  
• Bertschinger, E. 1986, ApJ, 304, 154 In the text NASA ADS  
• Chakrabarti, S. K. 1989, MNRAS, 240, 7 In the text NASA ADS  
• Chanaud, R. C., & Powell, A. 1965, J. Acoust. Soc. Am., 37, 902 In the text  
• D'Iakov, S. P. 1958, Sov. Phys. JETP, 6, 729, 739  
• Ffowcs Williams, J. E. 1969, Ann. Rev. Fluid. Mech., 1, 197 In the text  
• Foglizzo, T., & Ruffert, M. 1997, A&A, 320, 342 In the text NASA ADS  
• Foglizzo, T., & Tagger, M. 2000, A&A, 363, 174 (FT00) In the text NASA ADS  
• Foglizzo, T. 2001, A&A, 368, 311 (F01) In the text NASA ADS  
• Fukue, J. 1987, PASJ, 39, 309 In the text NASA ADS  
• Howe, M. S. 1975, J. Fluid Mech., 71, 625 In the text NASA ADS  
• Hoyle, F., & Lyttleton, R. A. 1939, Proc. Cam. Phil. Soc., 35, 405 In the text  
• Ishii, T., Matsuda, T., Shima, E., et al. 1993, ApJ, 404, 706 In the text NASA ADS  
• Kontorovich, V. M. 1958, Sov. Phys. JETP, 6, 729, 1180 In the text  
• Kontorovich, V. M. 1959, Sov. Phys. Acoust., 5, 320 In the text  
• Landau, L. D., & Lifshitz, E. M. 1987, Fluid Mechanics, 6 (Pergamon Press) In the text  
• McCrea, W. H. 1956, ApJ, 124, 461 In the text NASA ADS  
• Meszaros, P., & Ostriker, J. P. 1983, ApJ, 273, L59 In the text NASA ADS  
• Molteni, D., Tóth, G., & Kuznetsov, O. A. 1999, ApJ, 516, 411 In the text NASA ADS  
• Nakayama, K. 1992, MNRAS, 259, 259 In the text NASA ADS  
• Nakayama, K. 1993, PASJ, 45, 167 In the text NASA ADS  
• Nakayama, K. 1994, MNRAS, 270, 871 In the text NASA ADS  
• Nobuta, K., & Hanawa, T. 1994, PASJ, 46, 257 In the text NASA ADS  
• Pogorelov, N. V., Ohsugi, Y., & Matsuda, T. 2000, MNRAS, 313, 198 In the text NASA ADS  
• Protheroe, R. J., & Kazanas, D. 1983, ApJ, 265, 620 In the text NASA ADS  
• Rockwell, D. 1983, AIAA J., 21, 645 In the text  
• Ruffert, M. 1994, A&AS, 106, 505 In the text NASA ADS  
• Ruffert, M. 1996, A&A, 311, 817 In the text NASA ADS  
• Ruffert, M., & Arnett, D. 1994, ApJ, 427, 351 NASA ADS  
• Scharlemann, E. T. 1981, ApJ, 246, L15 In the text NASA ADS  
• Shima, E., Matsuda, T., Anzer, U., Börner, G., & Boffin, H. M. J. 1998, A&A, 337, 311 In the text NASA ADS  
• Vishniac, E. T. 1983, ApJ, 274, 152 In the text NASA ADS  
• Vishniac, E. T., & Ryu, D. 1989, ApJ, 337, 917 In the text NASA ADS