A&A 409, 1-7 (2003)
DOI: 10.1051/0004-6361:20031080

Non-axisymmetric instabilities in shocked accretion flows with differential rotation

Wei-Min Gu^1,2 - T. Foglizzo¹

1 - Service d'Astrophysique, CEA/DSM/DAPNIA, CEA-Saclay, 91191 Gif-sur-Yvette, France
2 - Department of Physics, Xiamen University, Xiamen, Fujian 361005, PR China

Received 10 April 2003 / Accepted 30 June 2003

Abstract
The linear stability of a shocked isothermal accretion flow onto a black hole is investigated in the inviscid limit. The outer shock solution, which was previously found to be stable with respect to axisymmetric perturbations, is, however, generally unstable to non-axisymmetric ones. Eigenmodes and growth rates are obtained by numerical integration of the linearized equations. The mechanism of this instability is based on the cycle of acoustic waves between their corotation radius and the shock. It is a form of the Papaloizou-Pringle instability, modified by advection and the presence of the shock. As such it can be generalized to non isothermal shocked accretion flows. Blobs and vortices are generated by the shock as a by-product of the instability.

Key words: accretion, accretion disks - black hole physics - hydrodynamics - instabilities - shock waves

1 Introduction

Hydrodynamic instabilities of shocked accretion flows may explain some of the properties of X-ray binaries, such as their time variability. The structure of stationary accretion flows involving shocks was described by Fukue (1987) and Chakarabarti (1989a,b). Even with the simple inviscid hypothesis, the structure of shocked accretion flows is complex, and their stability is not yet fully understood. As noted by Nakayama (1992), the calculations of Chakrabarti (1989a,b) are a study of the forced oscillations of the flow rather than an analysis of its intrinsic stability. Nakayama (1992, 1993) introduced a new type of global instability between the inner sonic point and the shock. He found that, of the two possible shock positions, the inner one is unstable due to post-shock acceleration, while the outer one is stable due to post-shock deceleration. His conclusion was confirmed by Nobuta & Hanawa (1994), whose numerical simulations showed that the inner shock is completely destroyed by perturbations, while the outer one is stable. All the above works, however, were based on axisymmetric calculations. Molteni et al. (1999, hereafter MTK) performed 2D simulations of an adiabatic flow with an outer shock and found a non-axisymmetric instability. They showed that the instability saturates at a low level, and a new asymmetric configuration develops, with a deformed shock rotating steadily. MTK pointed out that this effect may have relevant observational consequences, such as quasi-periodic oscillations (QPO). The mechanism of the instability was not explained by MTK, who briefly mentioned a possible link with the numerical simulations of Blaes & Hawley (1988). The Papaloizou-Pringle instability (1984, hereafter PPI) simulated by Blaes & Hawley is known to take place in discs or tori, in which the radial velocity is initially zero, whereas the flow simulated by MTK involves radial advection, an inner sonic point and an outer shock. The effect of advection on the PPI was investigated by Blaes (1987), who found that the PPI is strongly stabilized by advection at the inner boundary. The interpretation of the results of MTK in terms of the PPI is thus not obvious a priori. What is more, the flow studied by MTK is also potentially unstable by the advective-acoustic mechanism (Foglizzo & Tagger 2000; Foglizzo 2001, 2002), based on the cycle of entropy/vorticity perturbations and acoustic waves in the subsonic region between a stationary shock and a sonic surface.

The aim of this study is thus to understand the instability mechanism at work in shocked accretion flows. For the sake of simplicity, the present linear stability analysis is focused on isothermal accretion flows with constant angular momentum. The paper is organised as follows. Linearized equations and boundary conditions are described in Sect. 2 and solved numerically in Sect. 3. The instability mechanism is analysed in Sect. 4, and compared to the simulations by MTK in Sect. 5.

2 Equations

	Figure 1: Radial profile of the Mach number in the unperturbed flow, for which l=1.87 and $c_{\rm s}=0.1$. The two circles denote the inner and the outer sonic points, respectively. The dashed arrow shows the shock at $r_{\rm sh}= 6.8$.
Open with DEXTER

An inviscid, isothermal flow around a black hole is considered, in the pseudo-Newtonian potential introduced by Paczynski & Wiita (1980), $\Phi \equiv -GM/(r-r_{\rm g})$. Equations are made dimensionless by using the Schwarzschild radius and the speed of light as reference units, i.e., $r_{\rm g} \equiv 1$and $c \equiv 1$. In these units, the Keplerian frequency is denoted $\Omega_{\rm K}=1/(r-1)(2r)^{1\over2}$. In order to overcome the technical difficulty of treating a realistic vertical structure, three types of simplifying assumptions can be used; constant thickness, constant angle (conical flows), or vertical equilibrium. Chakrabarti & Das (2001) proved that there is no essential difference among these three assumptions from the point of view of the existence of stationary solutions. The thickness of the flow is approximated as a constant for the sake of simplicity, as in Nakayama (1992), Nobuta & Hanawa (1994), Blaes (1987) or MTK, although the approximate thickness $H\sim c_{\rm s}/\Omega_{\rm K}$ deduced from the balance of the vertical pressure force and the vertical gravity is not constant.

The stationary flow is described by the conservation of mass and the Bernoulli equation:

 

$$\rho r v_r= {\rm const.},$$ (1)

 

$${v_r^2\over 2}+{l^2\over 2r^2}+c_{\rm s}^2\log{\rho}-{1\over 2(r-1)}={\rm const.},$$ (2)

where $\rho$ is the density, $\upsilon_r$ is the radial velocity, $c_{\rm s}$ is the sound speed, and l is the specific angular momentum. An example solution of unperturbed flow is shown in Fig. 1, for which l=1.87 and $c_{\rm s}=0.1$. As a consequence of the inviscid hypothesis, stationary accretion flows onto a black hole are very sub-Keplerian in their outer parts. The question of the origin of this sub-Keplerian flow is still uncertain (Chakrabarti & Titarchuk 1995; Molteni et al. 2001). In Fig. 1, the unperturbed flow is supersonic between the outer sonic point ( $R_{\rm s2} = 44.5$) and the shock ( $r_{\rm sh}= 6.8$), becomes subsonic between the shock and the inner sonic point ( $R_{\rm s} = 2.36$) and goes into the central black hole supersonically. This example solution is precisely the one chosen by Nobuta & Hanawa (1994), in which they showed that the outer shock is stable to axisymmetric perturbations. Our numerical results in Sect. 3, however, indicate that it is unstable to non-axisymmetric perturbations.

The mass conservation equation and the Euler equation are written as follows:

 

$$\frac{\partial\rho}{\partial t}+\nabla\cdot (\rho\upsilon) = 0,$$ (3)

 

$$\frac{\partial\upsilon}{\partial t}+w\times \upsilon +\nabla\left... ...\frac{\upsilon ^2}{2} +c_{\rm s}^2\log \rho -\frac{1}{2(r-1)}\right\rbrack = 0,$$ (4)

where w is the vorticity. In order to write the linearized equations in the simplest form, the two functions f,gare defined as follows:

 

$$f \equiv \frac{\delta\rho}{\rho}+\frac{\upsilon_r\delta\upsilon_r}{c_{\rm s}^2} +\frac{\upsilon_\varphi\delta\upsilon_\varphi}{c_{\rm s}^2} ,$$ (5)

 

$$g \equiv \frac{\delta\rho}{\rho}+\frac{\delta\upsilon_r}{\upsilon_r} ,$$ (6)

where f is the perturbation of the Bernoulli constant and g is the perturbation of the mass accretion rate. The frequency $\omega'$ of the perturbation measured in the rotating frame is defined as:

$$\omega'\equiv \omega- m\Omega,$$ (7)

where m is the azimuthal wave number and $\Omega \equiv l/r^2$ is the angular velocity. With the standard method of linear stability analysis, (perturbations proportional to e $^{-i\omega t+im\varphi}$, see Appendix B), the following differential system is obtained:

  

$$\frac{\partial f}{\partial r} = \frac{i\omega{\cal M}^2}{\upsilon... ...al M}^2)} {\rm e}^{\int_{r_{\rm sh}}^r\frac{i\omega'}{\upsilon_r}{\rm d} r} \ ,$$ (8)

$$\frac{\partial g}{\partial r} = \frac{i\omega'}{\upsilon_r(1-{\ca... ...ht\rbrack {\rm e}^{\int_{r_{\rm sh}}^r\frac{i\omega'}{\upsilon_r}{\rm d} r} \ ,$$ (9)

where $B \equiv r\upsilon_rw_z$, w_z is the vorticity along the rotation axis, ${\cal M}\equiv -\upsilon_r/c_{\rm s}$ is the radial Mach number, and the subscript "sh'' denotes the shock position. The boundary conditions corresponding to a perturbed shock velocity $\Delta v_r$ are obtained in Appendix C:

   

$$f_{\rm sh} \left( \frac{\omega'}{\omega}-\frac{i{\cal M}c_{\rm s}\eta}{\omeg... ...\over1-{\cal M}^2}\right) (1-{\cal M}^2)^2\frac{\Delta \upsilon_r}{\upsilon_r},$$ = (10)

$$g_{\rm sh} \frac{\omega'}{\omega}(1-{\cal M}^2)\frac{\Delta \upsilon_r}{\upsilon_r},$$ = (11)

$$B_{\rm sh} -imc_{\rm s}^2\left(\frac{\omega'}{\omega}-\frac{i{\cal M}c_{\rm s}\eta}{\omega r}\right) (1-{\cal M}^2)^2\frac{\Delta \upsilon_r}{\upsilon_r},$$ = (12)

where $\eta \equiv ({\rm d}\log {\cal M}/ {\rm d}\log r)\vert _{\rm sh}$ (Eq. (A.1)), ${\cal M}$ and $\upsilon_r$ are calculated at the post-shock side.

The vorticity $\delta w_z$ produced by the shock is:

 

$$\delta w_z= \frac{B_{\rm sh}}{r\upsilon_r}{\rm e}^{\int_{r_{\rm sh}}^r\frac{i\omega'}{\upsilon_r}{\rm d} r}.$$ (13)

There are two differential equations for f,g and one unknown parameter $\omega$ in our system, thus three boundary conditions are needed to solve the equations. In addition to the two boundary conditions Eqs. (10), (11) at the shock, a third equation is obtained from the critical condition at the sonic point,

 

$$\omega g_{\rm son}-\omega'f_{\rm son}-\frac{ilB_{\rm sh}}{r_{\rm ... ...e}^{\int_{r_{\rm sh}}^{r_{\rm son}}\frac{i\omega'}{\upsilon_r}{\rm d} r}= 0 \ .$$ (14)

These three boundary conditions are used to numerically solve the differential system Eqs. (8), (9) and to determine the eigenfrequencies $\omega$. A single equation corresponding to this boundary value problem is formulated in Appendix D.

3 Numerical results

	Figure 2: Eigenspectrum of the flow l=1.87 and $c_{\rm s}=0.1$, showing the instability of the modes $1\le m\le 5$. The empty triangles correspond to the stable axisymmetric eigenmodes found by Nobuta & Hanawa (1994).
Open with DEXTER

	Figure 3: The maximum growth rate of the non radial modes is reached for $m\sim 10$in the flow l=1.87 and $c_{\rm s}=0.1$. The pattern frequency $\omega _r/m$ of these unstable modes is concentrated in a limited band.
Open with DEXTER

The standard Runge-Kutta method is used to integrate differential equations from the sonic point to the shock. Figure 2 shows the eigenspectrum of the flow l=1.87 and $c_{\rm s}=0.1$ studied by Nobuta & Hanawa (1994), for perturbations $0\le m \le 5$. The stability of axisymmetric perturbations confirms the results of Nobuta & Hanawa (1994). In contrast, non-axisymmetric perturbations are unstable. The highest growth rate is reached for m=10 perturbations, as seen in Fig. 3. Determining numerically how the most unstable m depends on the value of $(l,c_{\rm s})$ is beyond the scope of this paper. In what follows, numerical calculations are restricted to m = 1perturbations for the sake of simplicity.

Figure 4: The two thick solid lines are the threshold for shock-included solutions. The dotted lines measure the shock strength by the value of ${\cal M}_{\rm sh}$indicated on the left. The thin solid lines correspond to the value of the minimum Mach number ${\cal M}_{\rm min}$ indicated on the right. The 19 filled circles correspond to unstable outer shock solutions, andthe 5 empty squares correspond to stable ones.

Open with DEXTER

Figure 4 shows the domain $(l,c_{\rm s})$ of angular momentum and sound speed for which an isothermal inviscid flow, subsonic far from the accretor, may accrete onto a black-hole through a stationary shock. This domain is limited by two solid lines. The solid line on the left corresponds to the solutions in which the inner shock and the outer shock are identical, whereas the solid line on the right corresponds to extremely weak shocks, ${\cal M}_{\rm sh}\sim 1$. A total of 24 shocked solutions were considered for numerical calculation, in order to get a broad view of their stability properties. Since the inner shock solution was already found unstable to axisymmetric perturbations, we concentrated on the stability of the outer shock solution. Solutions at the bottom right corner in Fig. 4, for which the advection timescale is much longer than the acoustic one ( ${\cal M}_{\rm min} \le 0.01$), were avoided for computing time reasons. Among this sample of 24 shocked flows, 19 were found to be unstable and 5 to be stable. The results shown in Fig. 4 suggest that the shock is generally unstable to non-axisymmetric perturbations, except for a very narrow region close to the upper right border. In other words, the shock might be stable only for ${\cal M}_{\rm sh}\sim 1$.

4 Instability mechanism

4.1 Comparison with the advective timescale

	Figure 5: Growth rates compared to three timescales: advection $\tau _{\rm adv}$ (squares), acoustic $\tau _{\rm ac}$ (circles), and rotation $\Omega (r_{\rm son})$ (crosses).
Open with DEXTER

In an isothermal flow, two types of cycles may exist between the sonic point and the shock, i.e., the purely acoustic cycle and the vortical-acoustic cycle. From the results of Foglizzo (2002), we would expect the vortical-acoustic cycle to be particularly unstable for very strong shocks, ${\cal M}_{\rm sh}\ll 1$, with a growth time comparable to the advection timescale $\tau _{\rm adv}$:

$$\tau_{\rm adv} \equiv \int_{r_{\rm son}}^{r_{\rm sh}}\frac{{\rm d} r}{\vert\upsilon_r\vert}\cdot$$ (15)

The numerical results seem to exclude an explanation based on a vortical-acoustic cycle:

(i) the flow is generally unstable to non-axisymmetric perturbations, even for mild shocks: the range of shock strengths among the 19unstable flows considered is very wide ( $0.06\le{\cal M}_{\rm sh}\le0.94$).

(ii) the growth time can be much shorter than the advection time. The growth time of the instability is plotted in Fig. 5, in units of the rotation frequency at the sonic radius $\Omega (r_{\rm son})$, the advection time $\tau _{\rm adv}$ and also in units of the acoustic time $\tau _{\rm ac}$ approximated by the following lower bound:

$$\tau_{\rm ac} \equiv {2\over1-{\cal M}_{\rm min}^2}{r_{\rm sh}-r_{\rm son}\over c_{\rm s}}\cdot$$ (16)

The dispersion of the points in Fig. 5 is smallest when the growth rate is measured in units of the acoustic time, suggesting a purely acoustic cycle. If this is the case, what is the mechanism?

4.2 The corotation region in the PPI mechanism

The corotation radius ${r_{\rm co}}$ of the perturbation is defined by $\omega= m\Omega_0$, where $\Omega_0\equiv \Omega({r_{\rm co}})$:

$${r_{\rm co}}\equiv \left({lm\over\omega _r}\right)^{1\over2}\cdot$$ (17)

According to Fig. 6, the corotation radius of the most unstable modes is always located between the sonic point and the shock. Figures 5 and 6 are hints in favour of the PPI. The mechanism of the PPI was formulated most simply by Goldreich & Narayan 1985 (hereafter GN85) in thin discs. The more subtle effect of a corotation resonnance (Narayan et al. 1987; Kato 1987) can be neglected in the present study, since the epicyclic frequency is zero ( $\Omega=l/r^2$). The PPI is based on the exchange of energy and angular momentum between acoustic waves propagating inside and those propagating outside the corotation radius (Mark 1976). The thin disc hypothesis in these studies is technically important in order to treat high frequency perturbations in the WKB approximation $\omega=m\Omega_0\gg c_{\rm s}/r$. According to Fig. 7, the high frequency approximation is applicable only in the limit $c_{\rm s}\to 0$, which coincides with strong shocks according to Fig. 4. Following GN85 and Kato (1987), we believe that the physical mechanism captured analytically in thin discs is also relevant in thicker discs ( $\Omega_{\rm K} r/c_{\rm s}\sim 1$), even if a quantitative estimate of the growth rate is precluded.

	Figure 6: The position of the corotation radius is indicated for 19 unstable flows (circles) and 5 stable ones (squares).
Open with DEXTER

	Figure 7: The value of the ratio $c_{\rm s}/\Omega _{\rm K}r$ at the sonic point (dotted lines) and at the shock (full lines) is indicated on the lines.
Open with DEXTER

The PPI occurs with either an inner or an outer reflecting boundary, or even more efficiently with both. Let us recall this mechanism in the case of a single reflecting boundary, noting ${\cal T}_{\rm co}$ the fractional amplitude of the transmitted wave through the corotation zone, and ${\cal R}$ the fractional amplitude of the acoustic wave reflected by the boundary. An incident wave with energy +1transmits an energy $-\vert{\cal T}_{\rm co}\vert^2$ on the other side of corotation, thus amplifying the energy of the reflected wave by a factor $1+\vert{\cal T}_{\rm co}\vert^2$. If $\tau$ is the duration of the acoustic cycle, the growth rate $\omega_i$ of this instability is:

 

$$\omega_i={1\over\tau}\log\left\lbrack\vert{\cal R}\vert(1+\vert{\cal T}_{\rm co}\vert^2)^{1\over 2}\right\rbrack .$$ (18)

The following estimate of ${\cal T}_{\rm co}$ was obtained by GN85 and confirmed by Kato (1987):

 

$${\cal T}_{\rm co}\sim\exp\left\lbrack-{\pi\over 2} {c_{\rm s}\ove... ...t({m^2\over{r_{\rm co}}^2}+{\kappa^2\over c_{\rm s}^2}\right)\right\rbrack\cdot$$ (19)

The instability is thus slow in thin discs if $\kappa>0$, and the most unstable modes correspond to $m\sim \kappa {r_{\rm co}}/c_{\rm s}$. Let us now investigate the effects of radial advection on this mechanism.

4.3 Effect of advection on the corotation region

The effect of advection on the corotation region can be handled analytically in the high frequency limit. The second order differential equation satisfied by the homogeneous solution f₀ is obtained in Appendix E:

 

$$\left\lbrace{\partial^2 \over\partial r^2}+k^2\right\rbrace \left... ...omega'{{\cal M}\over1-{\cal M}^2}{{\rm d} r\over c_{\rm s}}}\right\rbrack=0 \ ,$$ (20)

In the high frequency limit, acoustic waves can propagate in the region where the function k²(r) is positive, and are evanescent in the region where k²(r)<0. k² is approximated at high frequency as follows:

 

$$k^2\sim {1\over(1-{\cal M}^2)} \left\lbrack{\omega'^2\over(1-{\cal M}^2)c_{\rm s}^2}-{m^2\over r^2}\right\rbrack\cdot$$ (21)

k² is positive near the sonic point ( ${\cal M}\sim1$), negative near the corotation radius ( $\omega'\sim0$), and positive again at large radius. A calculation with $\kappa=0$ and ${\cal M}>0$ in Appendix F shows that ${\cal T}_{\rm co}$ is independent of the Mach number in the high frequency approximation, if the gradient of the Mach number is neglected:

 

$${\cal T}_{\rm co}\sim\exp\left\lbrack-{\pi\over 2} {mc_{\rm s}\over \vert{\dot\Omega}_0\vert{r_{\rm co}}^2}\right\rbrack \ ,$$ (22)

where ${\dot\Omega}_0$ is the derivative of the rotation frequency at the corotation radius. This estimate at high frequency thus favours modes with a low azimuthal number m. Beyond the high frequency approximation, this calculation shows that the amplification mechanism at corotation can operate even with radial advection.

4.4 Crucial role of the boundaries

In order to estimate the effect of radial advection at the inner boundary, Blaes (1987) considered an inviscid flow with constant thickness in a pseudo-Newtonian potential, and a uniform angular momentum. It differed from the present study by an adiabatic hypothesis and, most importantly, a leaking outer boundary. As stressed by GN85, the PPI requires at least one reflecting boundary. The stabilizing effect of advection found by Blaes (1987) is thus directly related to the impossibility, for acoustic waves moving away from the corotation, to be reflected towards it. By contrast, a stationary shock standing at the outer boundary of the flow is a good reflector of acoustic waves, if the shock is not too weak. A continuity argument would lead to expect the reflection coefficient ${\cal R}_{\rm sh}$ to decrease to zero in the limit ${\cal M}_{\rm sh}\to 1$. This is confirmed by a calculation of ${\cal R}_{\rm sh}$ in the high frequency limit (see Appendix E), which gives the same result as in isothermal flows without rotation (Foglizzo 2002):

 

$${\cal R}_{\rm sh}\sim- {1-{\cal M}_{\rm sh}\over1+{\cal M}_{\rm sh}}\cdot$$ (23)

According to Eqs. (18) and (23), a strong isothermal shock (i.e. $\vert{\cal R}_{\rm sh}\vert\sim 1$) is a sufficient condition for the instability of the acoustic cycle between the shock and corotation in thin discs. Acoustic energy in shocked flows is thus trapped between the shock and the corotation radius, and may leak inside through the sonic radius without damping the instability. On the basis of the results obtained by Blaes (1987), a strong decrease of the growth rate should be expected in the weak shock limit. This trend is clear in Fig. 5, when measuring the growth rate, as Blaes (1987), in units of the rotation frequency at a fixed radius. What is more, all the stable modes correspond to weak shocks ( ${\cal M}_{\rm sh}>0.8$) according to Fig. 6.

	Figure 8: Growth rate measured in units of the total acoustic time $\tau _{\rm ac}$ (circles) compared to the acoustic time $\tau _{\rm out}$between the corotation and the shock (stars).
Open with DEXTER

Figure 8 shows that the dispersion of the points is further decreased by measuring the growth rate in units of the acoustic timescale $\tau _{\rm out}$ between the corotation and the shock, approximated by:

 

$$\tau_{\rm out}\sim {2\over1-{\cal M}_{\rm min}^2}{r_{\rm sh}-{r_{\rm co}}\over c_{\rm s}}\cdot$$ (24)

The combination of Eqs. (18), (22)-(24) gives an analytical estimate of the growth rate which is justified only in the high frequency limit, i.e. for $c_{\rm s}\to 0$ and strong shocks. A direct application of this estimate outside its range of validity is unable to explain the strong dispersion of the points in Fig. 8 for weak shocks. As noted in Sect. 4.2 with Fig. 7, solutions with a weak shock coincides with those with the largest parameter ( $c_{\rm s}/\Omega _{\rm K}r$), i.e.the least adapted to a high frequency approximation. It should be noted that part of the dispersion of the growth rate in Fig. 8 may come from the contribution of the vorticity perturbations, when the advection time is comparable to the acoustic time: they trigger additionnal acoustic waves which may add or substract to the acoustic power present in the acoustic cycle. Could this effect be essentially stabilizing? According to Fig. 9, for a given shock strength (i.e. ${\cal M}_{\rm sh}=0.9$), the flow seems to be stabilized by a strong advection.

	Figure 9: Incidence of advection on the PPI for weak shocks ( ${\cal M}_{\rm sh}=0.9$). The growth rate of the instability is measured in units of the acoustic time $\tau _{\rm out}$ between the corotation and the shock.
Open with DEXTER

5 Comparison with adiabatic simulations

Although non-linear processes might redistribute the energy of the unstable modes, it could be tempting to identify the frequency $\omega _r^{\rm sim}$, observed in the simulations, with the real part of the most unstable eigenmode. Based on five regular oscillation examples, MTK found that the oscillation period P is always in the same range, 60 < P < 280 despite the large variation in the rotation rate at the shock distance, $100 < 2\pi /\Omega _{\rm sh} < 2000$. We find that all the five regular oscillation examples listed in MTK satisfy the condition, $\Omega_{\rm sh} < \omega _r^{\rm sim} \equiv 2\pi /P < \Omega_{\rm son}$, which is a requirement of the PPI mechanism.

A quantitative comparison with the results of MTK is hampered by our choice of an isothermal flow rather than an adiabatic one, and also by our restriction to m=1 perturbations which are not necessarily the most unstable ones. A global understanding of which azimuthal wavenumber corresponds to the most unstable perturbation has not been reached yet. The only isothermal example considered in Fig. 3 with m>1perturbations favours $8\le m\le12$, which seems higher than the low azimuthal wavenumbers present in the simulations of MTK.

Despite these uncertainties, the present isothermal calculation is simple enough to illustrate the possibility of a PPI in a shocked accretion flow, occuring on a shorter time scale than any possible advective-acoustic cycle. This instability mechanism should thus also apply to adiabatic flows, including those simulated by MTK, if the shock is able to reflect acoustic waves. The calculation of ${\cal R}_{\rm sh}$ in the adiabatic hypothesis (Foglizzo & Tagger 2000, Eq. (C13)) would be unchanged by rotation and non radial perturbations, in the high frequency limit, for a strong shock:

 

$${\cal R}_{\rm sh}\sim -{2^{1\over2}-\gamma^{1\over2}(\gamma-1)^{1\over2} \over 2^{1\over2}+\gamma^{1\over2}(\gamma-1)^{1\over2}}\cdot$$ (25)

Comparing Eq. (25) with Eq. (23), the reflection coefficient of strong adiabatic shocks is significantly smaller than strong isothermal ones in the high frequency limit. The adiabatic hypothesis in the work of MTK may thus marginally influence the strength of the instability ( $\vert{\cal R}_{\rm sh}\vert\sim 0.36$ for $\gamma=4/3$).

A refined comparison of the growth rate and oscillation frequency with the results of MTK would require a linear stability analysis in an adiabatic flow involving also $m\ge 2$perturbations. This will be the subject of a forthcoming paper.

6 Conclusion

The main results can be summarized as follows:
1. Despite post-shock deceleration, the outer shock solution is generally unstable with respect to non-axisymmetric perturbations.
2. The growth time is comparable to the acoustic time if the shock is not too weak. The flow might be stable only when ${\cal M}_{\rm sh}$ is close to 1.
3. A WKB analysis of thin discs with advection proves that the acoustic cycle between the corotation radius and the shock is unstable if the shock is strong, thus extending the mechanism of the Papaloizou-Pringle instability to shocked discs with advection.
4. Extrapolating the results established in thin discs, the instability of thicker discs could be interpreted as a form of the Papaloizou-Pringle instability modified by advection. This interpretation is supported by the existence of a corotation radius within the disc, a growth time comparable to the acoustic time, and the role of the shock strength. The influence of advection in the specific case of weak shocks, however, is still unexplained (Fig. 9).

One particularity of this instability should be pointed out: although the instability mechanism relies on acoustic waves, vortical waves are continuously generated by the perturbed shock and advected to the accretor (Eqs. (12), (13)). These vortical waves (or the entropy/vorticity waves in the case of non isothermal flows) may play the role of the "blobs'' invoked frequently to model the time variability of X-ray binaries (e.g.references in Belloni et al. 2002). The observed properties of these blobs could therefore be used to test the validity of a model of shocked accretion.

Acknowledgements

The authors would like to thank Dr M. Tagger for discussions. W. Gu is grateful to Prof. Lu for his encouragement.

References

 
• Belloni, T., Psaltis, D., & van Der Klis, M. 2002, ApJ, 572, 392 In the text NASA ADS  
• Blaes, O. M. 1987, MNRAS, 227, 975 In the text NASA ADS  
• Blaes, O. M., & Hawley, J. F. 1988, ApJ 326, 277 In the text NASA ADS  
• Chakrabarti, S. K. 1989a, PASJ, 41, 1145 In the text NASA ADS  
• Chakrabarti, S. K. 1989b, MNRAS, 240, 7 NASA ADS  
• Chakrabarti, S. K., & Das, S., 2001, MNRAS, 327, 808 In the text NASA ADS  
• Chakrabarti, S. K., & Titarchuk, L., 1995, ApJ, 455, 623 In the text NASA ADS  
• Foglizzo, T., & Tagger M. 2000, A&A, 363, 174 In the text NASA ADS  
• Foglizzo, T. 2001, A&A, 368, 311 In the text NASA ADS  
• Foglizzo, T. 2002, A&A, 392, 353 In the text NASA ADS  
• Fukue, J. 1987, PASJ, 39, 309 In the text NASA ADS  
• Goldreich, P., & Narayan, R. 1985, MNRAS, 213, 7 (GN85) In the text  
• Kato, S. 1987, PASJ, 39, 645 In the text NASA ADS  
• Mark, J. W. K. 1976, ApJ, 205, 363 In the text NASA ADS  
• Molteni, D., Kuznetsov, O. A., Bisikalo, D. V., & Boyarchuk, A. A. 2001, MNRAS, 327, 1103 In the text NASA ADS  
• Molteni, D., Tóth, G., & Kuznetsov, O. A. 1999, ApJ, 516, 411 (MTK) 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  
• Narayan, R., Goldreich, P., & Goodman, J. 1987, MNRAS, 228, 1 In the text NASA ADS  
• Nobuta, K., & Hanawa, T. 1994, PASJ, 46, 257 In the text NASA ADS  
• Paczynski, B., & Wiita, P. J. 1980, A&A, 88, 23 In the text NASA ADS  
• Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 NASA ADS  
• Papaloizou, J. C. B., & Pringle, J. E. 1985, MNRAS, 213, 799 NASA ADS

  
Online Material

Appendix A: Description of the unperturbed flow

We consider a simple one-dimensional steady flow under the assumption of a constant thickness. The derivative of the Mach number is obtained by derivating the stationary flow Eqs (1), (2):

 

$$\eta\equiv \frac{{\rm d}\log {\cal M}}{{\rm d}\log r} = -\frac{1+(l^2-l_{\rm K}^2)/r^2c_{\rm s}^2} {1-{\cal M}^2} \ ,$$ (A.1)

where $l_{\rm K}$ is the Keplerian angular momentum, $l_{\rm K}^2 \equiv r^3/2(r-1)^2$. The continuity of the flow at the sonic point requires

$$r_{\rm son}^2c_{\rm s}^2+l^2-\frac{r_{\rm son}^3}{2(r_{\rm son}-1)^2} = 0.$$ (A.2)

The $(l,c_{\rm s})$ parameter space for the shocked accretion flows is shown in Fig. 4.

Appendix B: Linearized equations for perturbations

The mass conservation and Euler equations are linearized in order to obtain the following three equations:

$$\upsilon_r\frac{\partial g}{\partial r}+\frac{im}{r}\delta\upsilon_\varphi-i\omega'\frac{\delta\rho}{\rho} = 0 \ ,$$ (B.1)

$$c_s^2\frac{\partial f}{\partial r}-i\omega\delta\upsilon_r-\frac{l}{r}\delta w_z= 0 \ ,$$ (B.2)

$$\frac{imc_{\rm s}^2}{r}f-i\omega\delta\upsilon_\varphi+\upsilon_r\delta w_z= 0 \ .$$ (B.3)

The evolution of vorticity is deduced from the curl of Euler equation:

$$\frac{\partial w}{\partial t}+w(\nabla\cdot v)+(v\cdot \nabla)w-(w\cdot \nabla)v = 0 \ .$$ (B.4)

The vorticity can be integrated when linearized, leading to Eq. (13). The differential system of Eqs. (8), (9) is obtained from the above equations. The perturbation $\delta\upsilon_\varphi$ is related to f and $B_{\rm sh}$(Eqs. (B.3), (13)),

 

$$\delta\upsilon_\varphi= \frac{1}{i\omega r}\left(imc_{\rm s}^2f +... ...m sh}{\rm e}^{\int_{r_{\rm sh}}^r\frac{i\omega'}{\upsilon_r}{\rm d} r}\right) .$$ (B.5)

Equations (8), (9), (B.5) enable us to get $(\delta\upsilon_r,\delta\upsilon_\varphi,\delta\rho)$ from (f,g).

Appendix C: Boundary conditions at the shock

The isothermal shock conditions are written as follows:

$$\rho_-\upsilon_{r-} = \rho_+\upsilon_{r+} ,$$ (C.1)

$$\rho_-(\upsilon_{r-}^2+c_{\rm s}^2) = \rho_+(\upsilon_{r+}^2+c_{\rm s}^2) ,$$ (C.2)

$$\upsilon_{\varphi -} = \upsilon_{\varphi +} ,$$ (C.3)

where the subscripts "-'' and "+'' denote pre-shock and post-shock quantities, respectively. Let the shock be perturbed by $\Delta\xi$ in the radial direction:

$$\Delta\xi\propto {\rm e}^{-i\omega t+im\varphi}.$$ (C.4)

The perturbed velocity $\Delta \upsilon_r$ of the shock and its angle $\Delta\theta$ with the azimuthal direction are

$$\Delta \upsilon_r= \frac{\partial\Delta\xi}{\partial t} = -i\omega\Delta\xi,$$ (C.5)

$$\Delta\theta = \frac{1}{r}\frac{\partial\Delta\xi}{\partial\varphi} = \frac{im}{r}\Delta\xi= -\frac{m}{\omega r}\Delta \upsilon_r .$$ (C.6)

Considering the velocity component perpendicular to the shock, in the frame of the shock, the isothermal shock conditions become the following:

$$\left(\rho_- + \frac{\partial\rho_-}{\partial r}\Delta\xi\right) ... ...t) =\left(\rho_+ + \frac{\partial\rho_+}{\partial r}\Delta\xi+\delta\rho\right)$$

$$\qquad\qquad\times\left(\upsilon_{r+}+\frac{\partial\upsilon_{r+}... ...al r}\Delta\xi-\frac{\omega'}{\omega}\Delta \upsilon_r+\delta\upsilon_r\right),$$ (C.7)

$$\left(\rho_- + \frac{\partial\rho_-}{\partial r}\Delta\xi\right) ... ...\xi-\frac{\omega'}{\omega} \Delta \upsilon_r\right)^2 +c_{\rm s}^2\right\rbrack = \left(\rho_+ + \frac{\partial\rho_+}{\partial r}\Delta\xi+\delt... ...{\omega} \Delta \upsilon_r+\delta\upsilon_r\right)^2 +c_{\rm s}^2\right\rbrack,$$ (C.8)

$$\upsilon_{\varphi -} +\frac{\partial\upsilon_{\varphi -}}{\partia... ...phi +}}{\partial r}\Delta\xi+\upsilon_{r+}\Delta\theta +\delta\upsilon_\varphi.$$ (C.9)

The boundary conditions Eqs. (10), (12) are deduced from the above three equations.

Appendix D: Boundary value problem

With three boundary conditions Eqs. (10), (11), (14), the differential system could in principle be directly integrated using the Runge-Kutta method. The singularity at the sonic point requires to integrate from the sonic point towards the shock. The boundary value problem is conveniently reduced to a single equation involving the solution f₀ of the homogeneous system. The homogeneous differential equations are expressed as follows:

  

$$\frac{\partial f_0}{\partial r} \frac{i\omega{\cal M}^2}{\upsilon_r(1-{\cal M}^2)} \left(-\frac{\omega'}{\omega}f_0 + g_0\right),$$ = (D.1)

$$\frac{\partial g_0}{\partial r} \frac{i\omega'}{\upsilon_r(1-{\cal M}^2)} \left(\frac{\omega'}{\omega}f_0-{\cal M}^2g_0\right) - \frac{im^2c_{\rm s}^2}{\omega r^2\upsilon_r}f_0.$$ = (D.2)

At the sonic point, the following boundary conditions are obtained:

 

$$g_0\vert _{\rm son} =\frac{ \omega'}{\omega}f_0\vert _{\rm son} \ .$$ (D.3)

The derivation at the sonic point can be obtained from L'Hospital method,

   

$$\frac{\partial f_0}{\partial r}\vert _{\rm son} \frac{1}{2r_{\rm son}}\frac{\frac{\omega'^2r_{\rm son}^2} {c_{\rm... ...}}{\eta_{\rm son}+\frac{i\omega'r_{\rm son}}{c_{\rm s}}} f_0\vert _{\rm son}\ ,$$ = (D.4)

$$\frac{\partial g_0}{\partial r}\vert _{\rm son} \left(\frac{im^2c_{\rm s}}{\omega r^2}-\frac{i\omega'^2} {\omegac... ... son} -\frac{\omega'}{\omega}\frac{\partial f_0}{\partial r}\vert _{\rm son}\ ,$$ = (D.5)

$$\eta_{\rm son} -{1\over c_{\rm s}} \left\lbrack{r_{\rm son}(r_{\rm son}+1)\over ... ...\rm son}-1)^3}-\left({l\over r_{\rm son}}\right)^2\right \rbrack^{1\over2}\cdot$$ = (D.6)

With the method in "Appendix C'' of Foglizzo (2002), the following single equation for the boundary value problem is obtained:

$$\left(f\cdot\frac{\partial f_0}{\partial r}-\frac{\partial f}{\partial r}\cdot f_0\right)\vert _{\rm sh}$$ (D.7)

$$\quad= -\frac{\upsilon_r}{1-{\cal M}^2}B_{\rm sh}\int_{r_{\rm... ...silon_r}\frac{1+{\cal M}^2}{1-{\cal M}^2}{\rm d} r} {\rm d} r,$$ (D.8)

where the parameter $\lambda$ is:

$$\lambda\equiv \frac{1}{r^2c_{\rm s}^2\upsilon_r}\left\lbrack\frac... ...rac{2l}{\upsilon_r} \frac{{\rm d}\log(r\upsilon_r)}{{\rm d} r}\right\rbrack \ ,$$ (D.9)

and the value of $\partial f/\partial r$ at the shock is deduced from Eqs. (8) and (10)-(12).

Appendix E: Reflection at the shock

Perturbations $f_{\rm sh},g_{\rm sh}$ at the shock are decomposed on acoustic waves $f_0^{\pm},g_0^{\pm}$ and advected perturbations $f_{\rm B},g_{\rm B}$:

  

$$f_{\rm sh} f_0^++f_0^-+f_{\rm B},$$ = (E.1)

$$g_{\rm sh} g_0^++g_0^-+g_{\rm B}.$$ = (E.2)

The second order differential equation deduced from the differential system (D.1), (D.2) is

$${\partial^2 f_0\over\partial r^2}+\left\lbrace {\partial\log\over... ...{\rm s}}{2{\cal M}\over1-{\cal M}^2}\right\rbrace{\partial f_0\over \partial r}$$

$$+\left\lbrace\omega'^2-{m^2c_{\rm s}^2\over r^2}+iv_r{\partial\omega'\over\partial r}\right\rbrace {f_0\over(1-{\cal M}^2)c_{\rm s}^2}=0.$$ (E.3)

This equation is written in Eq. (20) in a canonical form as in Kato (1987), with

$$\equiv {\omega'^2\mu^2\over(1-{\cal M}^2)^2c_{\rm s}^2} -{1\over2}{\partial^2\over\partial r^2}\log{1-{\cal M}^2\over{\cal M}}$$                             k²

$$-\left({1\over2}{\partial\over\partial r}\log{1-{\cal M}^2\over{\cal M}}\right)^2 \cdot$$ (E.4)

$$\mu^2 \equiv 1-{m^2c_{\rm s}^2\over\omega'^2r^2}(1-{\cal M}^2) \ .$$ (E.5)

In the high frequency limit, the WKB approximations of the ingoing f₀⁺ and outgoing f₀^- acoustic solutions are:

 

$$f_0^{\pm}\propto\left({\omega{\cal M}\over\omega'\mu}\right)^{1\o... ...i\omega'{{\cal M}\mp\mu\over1-{\cal M}^2}{{\rm d} r\over c_{\rm s}}\right)\cdot$$ (E.6)

The definition of the acoustic flux used in Foglizzo (2002) can thus be extended to rotating flows as follows;

 

$$F_\pm\equiv {{\dot M}_0\over c^2}{\mu\omega'\over{\cal M}\omega}\vert f_\pm\vert^2 \ .$$ (E.7)

The variations of $f_0^{\pm}$ are dominated by the phase variations in the WKB limit:

$${\partial f_0^{\pm}\over\partial r}\sim {i\omega'\over c_{\rm s}}{{\cal M}\mp\mu\over1-{\cal M}^2}f_0^{\pm},$$ (E.8)

$$g_0^{\pm}\sim\pm{\omega'\mu\over\omega{\cal M}}f_0^{\pm}.$$ (E.9)

The dominant contribution of the advected vorticity to f and g, at high frequency, is deduced from the differential system (8), (9) where the derivative were replaced by a multiplication by $i\omega'/v_r$:

  

$$f_{\rm B} \equiv {mv_r^2-l\omega'\over\omega'^2r^2+m^2v^2}{iB_{\rm sh}\over c_{\rm s}^2},$$ (E.10)

$$g_{\rm B} \equiv {imB_{\rm sh}\over\omega'^2r^2+m^2v^2}\cdot$$ (E.11)

Using Eqs. (E.10), (E.11) and the boundary conditions (10)-(12), the linear system (E.1), (E.2) is solved in order to obtain $f_0^{\pm}$:

$$f_0^{\pm}={\Delta v_r\over 2v_r}{{\cal M}\over\mu}{(1-{\cal M}^2)... ...eta c_{\rm s}\over\omega' r}{1-{\cal M}^2\over{\cal M}\mp\mu}\right\rbrack\cdot$$ (E.12)

The reflection coefficient ${\cal R}_{\rm sh}$ is defined using Eq. (E.7), so that $\vert{\cal R}_{\rm sh}\vert^2$ is the fraction of the reflected acoustic flux reflected inward:

$${\cal R}_{\rm sh} \equiv {f_0^+\over f_0^-} \ ,$$ (E.13)

$$- {1+\mu{\cal M}\over1-\mu{\cal M}} \left({{\cal M}-\mu\over{\cal... ...} \over 1+{i\eta c_{\rm s}\over\omega' r}{1-{\cal M}^2\over{\cal M}+\mu} }\cdot$$ = (E.14)

In the high frequency limit, required for the validity of Eqs. (E.6) to (E.11), $\mu\sim 1$and the reflection coefficient is given by Eq. (23).

Appendix F: Approximation of the amplification at corotation

Defining a new variable X as in Foglizzo (2001),

$${{\rm d} X\over{\rm d} r}\equiv {{\cal M}\over1-{\cal M}^2},$$ (F.1)

the differential equation of acoustic perturbations is simply

$${\partial^2 {\tilde f}\over\partial X^2}+ \left\lbrack (\Omega-\O... ...er r^2}(1-{\cal M}^2)\right\rbrack {m^2{\tilde f}\over{\cal M}^2c_{\rm s}^2}=0.$$ (F.2)

This equation can be approximated as a parabolic cylinder differential equation in the vicinity of the corotation, by linearizing the rotation frequency and neglecting the variation of the Mach number in this region:

$$X\sim {{\cal M}_0\over 1-{\cal M}_0^2}(r-{r_{\rm co}}),$$ (F.3)

$${\partial^2 {\tilde f}\over\partial r^2}+ \left\lbrack {{\dot... ...\rm co}}^2}\right\rbrack {m^2{\tilde f}\over1-{\cal M}_0^2}=0,$$ (F.4)

where the subscript "0'' denotes the corotation radius. The effect of advection is obtained formally from GN85 with $\kappa=0$ by replacing the sound speed $c_{\rm s}$ by $c_{\rm s}(1-{\cal M}_0^2)^{1\over2}$, and the azimuthal wavenumber m by $m/(1-{\cal M}_0^2)^{1\over2}$. The fractional amplitude of the transmitted wave, deduced from Eqs. (5)-(8) of GN85 is thus independent of the Mach number in this approximation:

$${\cal T}_{\rm co}=\exp\left(-{\pi\over 2} {mc_{\rm s}\over \vert{\dot\Omega}_0\vert{r_{\rm co}}^2}\right)\cdot$$ (F.5)

An equivalent way to obtain the same result is the WKB approximation of the tunelling obtained by integrating the variation of the radial wavenumber, as in Kato (1987). Denoting by $r_{\rm IL},r_{\rm OL}$ the two zeros of kdeduced from Eq. (21),

 

$${\cal T}_{\rm co} \equiv \exp(-\Phi_{\rm LL}),$$ (F.6)

$$\Phi_{\rm LL} \equiv -i\int_{r_{\rm IL}}^{r_{\rm OL}}k{\rm d} r,$$ (F.7)

$$\sim \int_{r_{\rm IL}}^{r_{\rm OL}}\left\lbrack {m^2\over r^2}- {\omeg... ...c_{\rm s}^2}\right\rbrack^{1\over 2}{{\rm d} r\over (1-{\cal M}^2)^{1\over 2}},$$ (F.8)

$$\sim {\pi\over 2}{mc_{\rm s}\over \vert{\dot\Omega}_0\vert {r_{\rm co}}^2}\cdot$$ (F.9)