Corotational Instability of Inertial-Acoustic Modes in Black Hole Accretion Discs and Quasi-Periodic Oscillations

We study the global stability of non-axisymmetric p-modes (also called inertial-acoustic modes) trapped in the inner-most regions of accretion discs around black holes. We show that the lowest-order (highest-frequency) p-modes, with frequencies $\omega=(0.5-0.7) m\Omega_{\rm ISCO}$, can be overstable due to general relativistic effects, according to which the radial epicyclic frequency is a non-monotonic function of radius near the black hole. The mode is trapped inside the corotation resonance radius and carries a negative energy. The mode growth arises primarily from wave absorption at the corotation resonance, and the sign of the wave absorption depends on the gradient of the disc vortensity. When the mode frequency is sufficiently high, such that the slope of the vortensity is positive at corotation positive wave energy is absorbed at the resonance, leading to the growth of mode amplitude. We also study how the rapid radial inflow at the inner edge of the disc affects the mode trapping and growth. Our analysis of the behavior of the fluid perturbations in the transonic flow near the ISCO indicates that, while the inflow tends to damp the mode, the damping effect is sufficiently small under some conditions so that net mode growth can still be achieved. We further clarify the role of the Rossby wave instability and show that it does not operate for black hole accretion discs with smooth-varying vortensity profiles. Overstable non-axisymmetric p-modes driven by the corotational instability provide a plausible explanation for the high-frequency (>100 Hz) quasi-periodic oscillations (HFQPOs) observed from a number of black-hole X-ray binaries in the very high state. The absence of HFQPOs in the soft (thermal) state may result from mode damping due to the radial infall at the ISCO.

called "very high state"; see Done, Gierlinski & Kubota 2007), and they typically have high amplitudes and high coherence (Q > 10), and can vary in frequency on short timescales (minutes). However, it is the weaker, transient High-Frequency QPOs (HFQPOs, that have attracted more attention, since their frequencies do not vary significantly in response to sizable (factors of 3-4) luminosity changes and are comparable to the orbital frequencies at the Innermost Stable Circular Orbit (ISCO) of black holes with mass M ∼ 10M⊙. As such, HFQPOs potentially provide a probe to study the effects of strong gravity. HFQPOs are usually observed in the very high state of the X-ray binaries, and have low amplitudes (0.5 − 2% rms at 2-60 keV) and low coherence (Q ∼ 2 − 10). Out of the seven black-hole binaries from which HFQPOs have been reported, four show pairs of QPOs (first discovered in GRO J1655-40;Strohmayer 2001) with frequency ratios close to 2 : 3 (300 and 450 Hz in GRO J1655-40, 184 and 276 Hz in XTE J1550-564, 113 and 168 Hz in GRS 1915+105, 165 and 240 Hz in H1743-322; note that GRS 1915+105 also has a second pair of QPOs with f = 41 and 67 Hz).
It is worth noting that QPO (with period of ∼ 1 hour) in X-ray emission has recently been detected in the active galaxy RE J1034+396 (Gierlinski et al. 2008). This could be the "supermassive" analog of the HFQPOs detected in black-hole X-ray binaries.
Despite the observational progress, the origin of the HFQPOs remain unclear. A number of possibilities/models have been suggested or studied to various degrees of sophistication. We comment on some of these below: -Stella, Vietri & Morsink (1999) and others (see Schnittman & Bertschinger 2004, Schnittman 2005 suggested that orbiting hot spots (blobs) in the disc oscillating with epicyclic frequencies may provide variability in the X-ray emission. However the radial positions of such blobs are free parameters, which must be tuned to match the observed QPO frequencies, and it is also not clear that the blobs can survive the differential rotation of the disc.
- Abramowicz & Kluzniak (2001) suggested that HFQPOs involve certain nonlinear resonant phenomenon in the disc (e.g., coupling between the radial and vertical epicyclic oscillations of the disc fluid element; Kluzniak & Abramowicz 2002). This was motivated by the observed stability of the QPO frequencies and the commensurate frequency ratio. However, so far analysis has been done based only on toy models involving coupled harmonic oscillators (e.g. Rebusco 2004;Horak & Karas 2006) and no fluid dynamical model producing these resonances has been developed (see Abramowicz et al 2007 andRebusco 2008 for recent reviews). Petri (2008) considered the resonant oscillation of a test mass in the presence of a spiral density wave, but the origin of the wave is unclear.
-Acoustic oscillation modes in pressure-supported accretion tori have been suggested as a possible source of the observed QPOs (Rezzolla et al. 2003;Lee, Abramowicz & Kluziniak 2004; see also Blaes, Arras & Fragile 2006;Schnittman & Rezzolla 2006, Blaes et al 2007, Sramkova et al 2007. In this model, the commensurate mode frequencies arise from matching the radial wavelength to the size of the torus. Note that the QPO frequencies are determined mainly by the radial boundaries of the torus, which must be tuned to match the observed QPO frequencies. It is also not clear that the accretion flow in the very high state (in which HFQPOs are observed) is well represented by such a torus (e.g. Done et al 2007).
- Li & Narayan (2004) considered the dynamics of the interface between the accretion disc and the magnetosphere of a central compact object (see also Lovelace & Romanova 2007). The interface is generally Rayleigh-Tayor unstable and may also be Kelvin-Helmholtz unstable. While such an interface is clearly relevant to accreting magnetic neutron stars, Li & Narayan suggested that it may also be relevant to accreting black holes and that the strongly unstable interface modes may give rise to QPOs with commensurate frequencies.
-Perhaps the theoretically most appealing is the relativistic diskoseismic oscillation model, according to which general relativistic (GR) effects produce trapped oscillation modes in the inner region of the disc (Kato & Fukue 1980;Okazaki et al. 1987;Nowak & Wagoner 1991;see Wagoner 1999;Kato 2001 for reviews; see also Tassev & Bertschinger 2007 for the kinematic description of some of these wave modes). A large majority of previous studies have focused on disc g-modes (also called inertial modes or inertial-gravity modes, whose wavefunctions -such as the pressure perturbation, contain nodes in the vertical direction), because the trapping of the g-mode does not require a reflective inner/outer disc boundary. Kato (2003a) and Li, Goodman & Narayan (2003) showed that the g-mode that contains a corotation resonance (where the wave patten frequency equals the rotation rate of the background flow) in the wave zone is heavily damped. Thus the only nonaxisymmetric (m = 0) g-modes of interest are those trapped around the maximum of Ω + κ/m (where Ω is the rotational frequency, κ is the radial epicyclic frequency and m is the azimuthal mode number; see Fig. 1 below). Unfortunately, the frequencies of such modes, ω ≃ mΩISCO, are too high (by a factor of 2-3) compared to the observed values, given the measured mass and the estimated spin parameter of the black hole (Silbergleit & Wagoner 2007; see also Tassev & Bertschinger 2007). Axisymmetric g-modes (m = 0) may still be viable in the respect, and recent studies showed that they can be resonantly excited by global disc deformations through nonlinear effects (Kato 2003a(Kato ,2008Ferreira & Ogilvie 2008). Numerical simulations (Arras, Blaes & Turner 2006;Reynolds & Miller 2008), however, indicated that while axisymmetric g-mode oscillations are present in the hydrodynamic disc with no magnetic field, they disappear in the magnetic disc where MHD turbulence develops. Also, Fu & Lai (2008) carried out an analytic study of the effect of magnetic fields on diskoseismic modes and showed that even a weak (sub-thermal) magnetic field can "destroy" the self-trapping zone of disc g-modes, and this may (at least partly) explain the disappearance of the g-modes in the MHD simulations.
- Tagger  for a review) developed the theory of accretion-ejection instability for discs threaded by strong (of order or stronger than equipartition), large-scale poloidal magnetic fields. They showed that such magnetic field provides a strong coupling between spiral density waves and Rossby waves at the corotation, leading to the growth of the waves and energy ejection to disc corona.  suggested that normal modes trapped in the inner region of the disc become strongly unstable by a combination of accretion-ejection instability and an MHD form of the Rossby wave instability (see Lovelace et al. 1999;Li et al. 2000; see section 6 below). The Tagger model has the appealing feature that the instability leads to energy ejection into the disc corona, and thus explains why HFQPOs manifest mainly as the variations of the nonthermal (power-law) radiation from the systems.

This Paper
In this paper we study the global corotational instability of nonaxisymmetric p-modes (also called inertial-acoustic modes) trapped in the inner-most region of the accretion disc around a black hole. The p-modes do not have vertical structure (i.e., the wavefunctions have no node in the vertical direction). We focus on these modes because their basic wave properties (e.g. propagation diagram) are not affected qualitatively by disc magnetic fields (Fu & Lai 2008) and they are probably robust under hydromagnetic effects and disc turbulence (see Reynolds & Miller 2008).
The corotational instability of p-modes studied in this paper relies on the well-known GR effect of test-mass orbit around a black hole: Near the black hole, the radial epicyclic frequency κ reaches a maximum (at r = 8GM/c 2 for a Schwarzschild black hole) and goes to zero at the ISCO (rISCO = 6GM/c 2 ). This causes non-monotonic behavior in the fluid vortensity, ζ = κ 2 /(2ΣΩ) (assuming the surface density Σ is relatively smooth), such that dζ/dr > 0 for r < r peak (where r peak is the radius where ζ peaks) and dζ/dr < 0 for r > r peak . The vortensity gradient dζ/dr plays an important role in wave absorption at the corotation resonance (Tsang & Lai 2008a; see also Goldreich & Tremaine 1979 for corotational wave absorption due to external forcing). We show that the p-modes with frequencies such that the corotation radii lie inside the vortensity peak can grow in amplitude due to absorption at the corotation resonance, and that the overstability can be achieved for several modes with frequencies closely commensurate with the azimuthal wavenumber m.  have studied similar modes in discs threaded by strong magnetic fields, but in our analysis magnetic fields play no role.
The trapping of the p-modes requires the existence of a (partially) reflecting boundary at the disc inner edge, close to the ISCO. One may suspect that the rapid radial inflow at the ISCO will diminish any potential instabilities in the inner accretion disc (see Blaes 1987 for the case of thick accreting tori). Our analysis of the wave perturbations in the transonic accretion flow (see section 5) suggests that waves are partially reflected at the sonic point, and global overstable p-modes may still be produced under certain conditions (e.g., when the surface density of the flow varies on sufficiently small length scale around the sonic point). Even better mode trapping (and therefore larger mode growth) may be achieved when the system is an accretion state such that the inner disc edge does not behave as a zero-torque boundary (see section 7 for discussion and references).
Our paper is organized as follows. After summarizing the basic fluid equations for our problem in section 2, we give a physical discussion of the origin of the corotational instability of disc p-modes in section 3. We present in section 4 our calculations of the growing p-modes with simple reflective inner disc boundary conditions. Section 5 contains our analysis of the effect of the transonic radial inflow at the ISCO on the p-mode growth rate. In section 6, we discuss the role of the Rossby wave instability and show that it is not effective in typical accretion discs under consideration. In section 7 we discuss the application of our results to HFQPOs in black hole X-ray binaries.

SETUP AND BASIC EQUATIONS
We consider a geometrically thin disc and adopt cylindrical coordinate system (r, φ, z). The flow is assumed to be barotropic, so that the vertically integrated pressure, P = R p dz, depends only on the surface density, Σ = R ρ dz. We use the pseudo-Newtonian potential of Paczynski & Wiita (1980) with rS = 2GM/c 2 the Schwarzschild radius. The free-particle (Keplerian) orbital and radial epicyclic (angular) frequencies are ΩK = " 1 r dΦ dr The function κ peaks at r = (2 + √ 3)rS and declines to zero at rISCO = 3rS (while for a Schwarzschild black hole in GR, κ peaks at r = 4rS). The unperturbed flow has velocity u0 = (ur, rΩ, 0). Since pressure is negligible for thin discs, we have Ω ≃ ΩK .
Neglecting the self-gravity of the disc we have the linear perturbation equations: where δΣ, δu and δh = δP/Σ are the (Eulerian) perturbations of surface density, velocity and enthalpy, respectively. For barotropic flow, δh and δΣ are related by where cs is the sound speed, with c 2 s = dP/dΣ. We assume all perturbed quantities to be of the form e imφ−iωt , where m is a positive integer, and ω is the wave (angular) frequency. The perturbation equations then become is the wave frequency in the frame corotating with the unperturbed fluid. Except very near the inner edge of the disc, rin ≃ rISCO, the unperturbed radial velocity is small, |ur| ≪ rΩ. In our calculations of the disc modes, we will neglect ur and set the last terms on the left-hand sides of equations (6)-(8) to zero (However, ur plays an important role in determining the inner boundary condition of the fluid perturbations at rin; see section 5). Eliminating the velocity perturbations in favor of the enthalpy, we obtain our master equation where For concreteness we assume the surface density to have a the power-law form where p is the density index. The above equations adequately describe disc p-modes (also called inertial-acoustic modes), which do not have vertical structure. Other disc modes (g-modes and c-modes) involve the vertical degree of freedom [see Kato 2001 for a review; also see Fig. 1 of Fu & Lai (2008) for a quick summary], and their stability properties are studied by Kato (2003a), Li et al. (2003) and Tsang & Lai (2008b).
To determine the global modes of the disc, appropriate boundary conditions must be specified. These are discussed in sections 4 and 5.

P-MODES AND THEIR GROWTH DUE TO COROTATION RESONANCE: A PHYSICAL DISCUSSION
A WKB analysis of the wave equation (10) yields the dispersion relation for the local plane wave δh ∝ expˆi R r k(s)ds˜: Far from the singularity (ω = 0) at the corotation radius rc, this reduces to the well-known dispersion relation of spiral density wave with no self-gravity (e.g., Shu 1992), Density waves (p-modes) can propagate inside the inner Lindblad resonance radius rIL (defined byω = −κ), and outside the outer Lindblad resonance radius rOL (defined byω = κ), i.e., in the region where ω/m < Ω − κ/m and ω/m > Ω + κ/m, respectively (see Fig. 1). Between rIL and rOL, waves are evanescent except that a very narrow Rossby wave zone exists around the corotation radius. Indeed, in the vicinity ofω = 0, equation (13) reduces tõ Figure 1. Wave propagation diagram for non-axisymmetric p-modes in thin accretion discs around black holes. In the upper panel, the three solid curves depict the disc rotation profile Ω(r) and Ω ± κ/m (where κ is the radia epicyclic frequency); note that the three curves join each other at the disc inner radius r in = r ISCO (the inner-most circular orbit) since κ(r ISCO ) = 0. The wavy lines (of height ω/m) indicate the propagation zones for inertial-acoustic waves. Disc p-modes are trapped between r in and the inner Lindblad resonance radius (where ω/m = Ω − κ/m), but can tunnel through the corotation barrier. The lower panel depicts the disc vortensity profile, ζ = κ 2 /(2ΩΣ), which has a maximum at the radius r peak (as long as the surface density Σ does not vary too strongly with r). P-modes with ω/m > Ω peak = Ω(r peak ) (the upper wavy line) are overstable due to wave absorption at the corotation resonance radius rc (where ω/m = Ω) since (dζ/dr)r c > 0. P-modes with ω/m < Ω peak (the lower wavy line) tend to be damped by wave absorption at rc since (dζ/dr)r c < 0. Note that a narrow Rossby wave zone (labeled by thick horizontal bars) exists just outside or inside the corotation radius -the location of this Rossby zone determines the sign of the corotational wave absorption. where is the vortensity of the (unperturbed) flow. For (dζ/dr)r c > 0, the Rossby wave zone lies between rc and rc + ∆rR, where ∆rR = (2cs/κ)|ν| and R rc+∆r R rc k dr = πν, with the number of wavelengths in the Rossby zone given by (Tsang & where q ≡ − (d ln Ω/d ln r) rc , and the second equality assumes Σ ∝ r −p . For (dζ/dr)r c < 0, the Rossby wave zone lies inside rc, between rc − |∆rR| and rc (see Fig. 1). Note that since ν ∼ cs/(κr) ∼ H/r ≪ 1, no standing Rossby wave can exist in the Rossby zone (see also section 6).
Assuming that there exists a reflecting boundary at the inner disc radius rin ≃ rISCO (see sections 4.3 and 5), normal modes can be produced, with the waves partially trapped between rin and rIL -these are the p-modes that we will focus on in this paper. The mode eigen-frequency ω = ωr + iωi is generally complex, with the real part ωr determined approximately by the Sommerfeld "quantization" condition whereωr = ωr − mΩ, n is an integer and ϕ (of order unity) is a phase factor depending on the details of the (inner and outer) boundary conditions. The overstability of the p-mode is directly related to the reflectivity of the corotation barrier between rIL and rOL. In the WKB approximation, the imaginary part of the mode frequency is given by (Tsang & Lai 2008a; see also Narayan et al. 1987, who considered shearing-sheet model) where R is the reflectivity (see below). Thus the mode becomes overstable (ωi > 0) for |R| > 1 (termed "super-reflection") and stable (ωi < 0) for |R| < 1. Super-reflection in fluid discs arises because the waves inside the corotation radius and those outside carry energy or angular momentum of different signs: Since the wave inside rc has pattern speed ωr/m less than the fluid rotation rate Ω(r), it carries negative energy; outside rc, we have ωr/m > Ω(r), the wave carries positive energy. Consider an incident wave δh ∝ exp(−i R r k dr), carrying energy of the amount (−1), propagating from small radii toward the corotation barrier 1 . The wave reflected at rIL takes the form δh ∝ R exp(i R r k dr), and the transmitted wave in the region r > rOL is δh ∝ T exp(i R r k dr). Because of the corotation singularity, the wave energy can also be transferred to the background flow and dissipated at the corotation radius. Energy conservation then gives −1 = (−1)|R| 2 + |T | 2 + Dc, or where Dc is the wave energy dissipated at the corotation. Tsang & Lai (2008a) derived the analytical expressions (in the WKB approximation) for T , R and Dc. Two effects determine the reflectivity. (i) The transmitted wave (corresponding to the |T | 2 term) always carries away positive energy and thus increases |R| 2 . (ii) Wave absorption at the corotation can have both signs, depending on ν: For ν > 0, the Rossby wave zone lies outside rc, positive wave energy is dissipated and we have Dc > 0; for ν < 0, the Rossby zone lies inside rc and we have Dc < 0 (see Fig. 1). Tsang & Lai (2008a) showed explicitly that under most conditions, |Dc| ≫ |T | 2 (except when ν ≃ 0, for which Dc ≃ 0). In the limit of |ν| ≪ 1, we have where ΘII ≡ [these expressions are valid for ΘII, ΘIIa ≫ 1; see Tsang & Lai (2008a) for more general expressions]. Thus super-reflectivity (|R| 2 > 1) and growing modes (ωi > 0) are achieved when Note that typically |νcrit| ≪ 1; if |T | 2 is neglected compared to Dc, then νcrit = 0. As mentioned before, since κ is non-monotonic near the black hole, the vortensity ζ is also non-monotonic, attaining a peak value at r = r peak before dropping to zero at the ISCO. Therefore, p-modes with frequencies such that the corotation radius rc lies inside r peak are expected to be overstable by the corotational instability discussed above. In other words, when ωr/m > Ω peak ≡ Ω(r peak ), the corotation resonance acts to grow the mode. Note that Ω peak depends on the surface density profile as well as the spacetime curvature around the black hole (see Fig. 2). On the other hand, when ωr/m < Ω peak (ν < 0), the corotational wave absorption acts to damp the mode. However, when ωr/m is only slightly smaller than Ω peak (νcrit < ν < 0) mode growth can still be obtained due to wave leakage beyond the outer Lindblad resonance, though the growth rate will be small (see section 4.4 for examples).

CALCULATIONS OF TRAPPED, OVERSTABLE P-MODES
To determine the eigenvalues ωr and ωi of the trapped modes, we solve equations (6)-(8) (with ur = 0) or equation (10) subjected to appropriate boundary conditions at rin and rout. Figure 2. Critical mode frequency for corotational instability as a function of the disc surface density index p (with Σ ∝ r −p ). Wave absorption at the corotation resonance acts to grow the mode only if the corotation occurs in the region of positive vortensity gradient, i.e., if the mode pattern frequency ω/m > Ω peak (see Fig. 1). P-mode trapping also requires ω/m < Ω ISCO = Ω(r ISCO ). Note that weak mode growth can still occur when ω/m is slightly below Ω peak due to wave leakage beyond the outer Lindblad resonance radius. See text for detail.

"Landau" Integration Contour
When solving eigenvalue problem using the standard method (e.g. the shooting method as described in Press et al 1998), we encountered a conundrum: For ν > 0, we could find both a growing mode and a decaying mode, with almost the same ωr but opposite ωi. This appears to contradict our discussion in section 3. This conundrum arises because our numerical integration is confined to the real r axis. However, analogous to Landau's analysis of wave damping in a plasma (e.g., Lifshitz & Pitaevskii 1981), care must be taken in defining appropriate contour of integration across the corotation resonance. Indeed, at corotation, equation (10) contains a singular term, proportional to where Rc ≡ rc − ircωi/(qωr) is the complex pole, rc is determined by ωr = mΩ(rc) and q = −(d ln Ω/d ln r)c > 0. As discussed in Lin (1955) in the context of hydrodynamical shear flows, to obtain physically relevant solutions of the fluid system, it is necessary that the integration contour lies above the pole. This is the Landau contour. In essence, only by adopting such a Landau contour can one obtain the correct wave absorption (dissipation) at the corotation. For growing modes (ωi > 0), Im(Rc) < 0, our numerical integration along the real r axis constitutes the correct Landau contour. On the other hand, for decaying modes (ωi < 0), the real r axis is not the correct Landau contour as Im(Rc) > 0. Instead, to obtain physical solutions for these decaying modes, the integration contour must be deformed so that Rc lies below it (see Fig. 3). As we are primarily interested in over-stable modes in this paper, it is adequate to integrate along the real r axis in our calculation.

Re(r)
c i

Re(r)
c i Im(r) Im(r) Figure 3. "Landau" contour for integration across the corotation resonance. To calculate the growing mode (ω i > 0), it is adequate to integrate the fluid perturbation equations along the real r axis (upper panel). To obtain the physical solution for the shrinking mode (ω i < 0), the integration contour must deformed so that the corotational pole Rc lies below the contour.

Outer Boundary Condition
As we are interested in self-excited modes in the inner region of the disc, we adopt the radiative outer boundary condition. Specifically, far from the outer Lindblad resonance (r > rOL) we demand that only an outgoing wave exists: where k = p −D/c 2 s (see Tsang & Lai 2008a). This gives the boundary condition at some rout > rOL: In practice, we find that rout ∼ 2rOL would yield sufficiently accurate results.

Inner Boundary Conditions
To obtain global trapped modes, at least partial wave reflection must occur at rin. To focus on the effect of corotational instability discussed in section 3, in this section we consider two simple inner boundary conditions. We defer our analysis of the effect of radial inflow on the p-modes to section 5. (i) At the ISCO, the flow plunges into the black hole, we expect a sudden decrease in the surface density of the disc. Thus, it is reasonable to consider the free surface boundary condition, i.e., the Lagrangian pressure perturbation ∆P = 0. Using δur = −iωξr (where ξr is the radial Lagrangian displacement), and ∆P = δP + ξrdP/dr, we have where we have used dP/dr = (dP/dΣ)(dΣ/dr) = −pc 2 s Σ/r for barotropic, power-law discs (Σ ∝ r −p ). (ii) We assume that the radial velocity perturbation vanishes at the inner boundary, i.e., δur = 0. This was adopted by  in their calculations of overstable global modes due to accretion-ejection instability.
Both of these boundary conditions correspond to zero loss of wave energy at the inner boundary: If a wave from large radii impinges toward rin, the reflected wave will have the same amplitude. However, the phase shifts due to reflection differ in the two cases, and the resulting mode frequencies ωr are different. Since the corotational wave amplification depends on ωr, the mode growth rate ωi will also be different. The radius r is in units of GM/c 2 . The disc sound speed is cs = 0.1rΩ, and the m = 2 modes are obtained using the inner boundary condition ∆P (r ISCO ) = 0. The left panels show the p-mode for the disc with a density profile Σ ∝ r −1 , with the eigenvalues ωr = 0.467mΩ ISCO , ω i /ωr = 0.0029 [where Ω ISCO = Ω(r ISCO )]; the right panels show the mode for the disc with constant surface density profile (p = 0), with eigenvalues ωr = 0.464mΩ ISCO , ω i /ωr = 0.00073. Note that the model shown on the left panels has rc < r peak (the radius of peak vortensity) and thus F (rc+) < F (rc−), while the model shown on the right panels has rc > r peak and F (rc+) > F (rc−). In both models there is a positive flux for r > rc due to the outward propagating wave. The inserts on the left panels show the blowups of the real wavefunctions near the corotation radius. Figure 5. The real and imaginary frequencies of disc p-modes (with azimuthal wave numbers m = 1, 2, 3) as a function of the surface density index p (where Σ ∝ r −p ). The modes are calculated assuming the inner boundary condition ∆P (r ISCO ) = 0. The left panels are for discs with cs = 0.1rΩ and the right panels for cs = 0.2rΩ. The dotted lines denote the lower bound ωr/m = Ω peak for which the corotational wave absorption acts to enhance mode growth.

Numerical Results
We solve for the complex eigen-frequency ω = ωr + iωi using the shooting method, with a fifth-order Runge-Kutta integrator (Press et al. 1992). As discussed in section 4.1 we only calculate the growing modes (ωi > 0). We consider disc models with different surface density profile (characterized by the index p), sound speed cs, and inner boundary conditions. For a given set of disc parameters and the azimuthal mode wavenumber m, the lowest order (highest frequency) mode has the best chance of being overstable. This is easily understood from our discussion in section 3 (see Fig. 1): a low-frequency wave has to penetrate a wider evanescent barrier for the corotational amplifier to be effective, and when ωr < mΩ peak the corotation resonance acts to damp the mode. For most disc models we have considered, the lowest order mode (of a given m) is the only mode that has ωi > 0.
Our numerical results are presented in Figures 4-7. Figure 4 gives two examples of the eigenfunctions of overstable trapped p-modes, obtained with the inner boundary condition ∆P = 0. In addition to δh and δur, we also plot the angular momentum flux carried by the wave across the disc (e.g., Goldreich & Tremaine 1979;Zhang & Lai 2006) where the second equality follows from equations (6)-(8) (with ur = 0). We see from Fig. 4 that outside the corotation radius (rc), F is nearly constant since only the outgoing wave exists in this region and the wave action is conserved in the limit of ωi ≪ ωr. Inside rc, the interference between the ingoing and outgoing waves gives rise to the variation of F . At rISCO, F approaches zero since no wave action is lost through the disc inner boundary when ∆P = 0. 2 Figure 4 also shows a flux jump across the corotation. In the limit of ωi ≪ ωr, dδh/dr, δur and δu φ are discontinuous across rc (although δh is continuous), giving rise to the flux discontinuity (see Tsang & Lai 2008a): This discontinuity signifies the corotational wave absorption, the sign of which depends on ν ∝ (dζ/dr)r c , as discussed in section 3. Thus, the model shown on the left panels of Fig. 4 has rc < r peak and ν > 0, and the mode growth is primarily driven by wave absorption at the corotation. The model shown on the right panels of Fig. 4, on the other hand, has rc > r peak and ν < 0, thus the corotational wave absorption acts to damp the mode, and the overall growth of the mode is due to the outgoing wave beyond rc, and the growth rate is much smaller than the model shown on the left panels. Figure 7. The real and imaginary frequencies of disc p-modes (with azimuthal wave numbers m = 1, 2, 3) as a function of the surface density index p (where Σ ∝ r −p ). The modes are calculated assuming the inner boundary condition δur(r ISCO ) = 0. The left panels are for discs with cs = 0.1rΩ and the right panels for cs = 0.2rΩ. The dotted lines denote the lower bound ω/m = Ω peak for which the corotational wave absorption acts to enhance mode growth.
Figures 5-6 show the frequencies of the fundamental (no node/highest frequency) growing p-modes (with m = 1, 2, 3) for various disc parameters, again obtained with the inner boundary condition ∆P = 0. We consider p in the range between −1.5 and 1.5, and cs up to 0.3rΩ. For a given sound speed, the real mode frequency ωr depends very weakly on p, but the growth rate ωi increases with p (see Fig. 5) since a larger value of p leads to a larger ν and enhanced wave absorption at the corotation [see equation (17)]. In general, as the sound speed increases, the effective wavelength of the mode increases, and ωr decreases in order "fit in" the trapping zone between rin and rIL (see Fig. 6). The mode growth rates ωi depends on cs in a non-monotonic way because of two competing effects: As cs increases, less attenuation occurs in the evanescent zone, and more wave energy can be absorbed at the corotation and propagate to the outer edge of the disc; these tend to increase ωi. On the other hand, increasing cs also leads to smaller ωr, which shifts the corotation resonance to a larger radius and leads to decreasing ν and ωi.
Note that the growing modes shown in Fig. 5 extend below the ωr/m = Ω peak boundary due to the propagation of waves beyond the corotation radius, as discussed in section 3 [see eqs. (23)]. Such modes (with νcrit < ν < 0) grow significantly slower than the modes with ν > 0 as the flux is attenuated by the entire barrier between rIL and rOL.
For comparison, Figure 7 shows the disc mode frequencies and growth rates when the inner boundary condition δur = 0 is adopted. The different boundary condition leads to a different phase shift ϕ and higher mode frequency, but the results are similar to those illustrated in Fig. 5. In particular, as p increases, ωr remains approximately constant while ωi increases.

EFFECT OF RADIAL INFLOW ON THE P-MODE GROWTH RATE
Our mode calculations presented in section 4 neglect the radial velocity of the accretion flow and assume a loss-less inner disc boundary condition (either ∆P = 0 or δur = 0 at rISCO). In real discs, the radial inflow velocity ur is not negligible as r approaches rISCO, and the flow goes through a transonic point (where ur = −cs) at a radius very close to rISCO. We expect that part of the fluid perturbations may be advected into the black hole and the inner disc boundary will not be completely loss-less. Here we study the effect of the transonic flow on the p-mode growth rate.
We note that just as the accretion disc is not laminar but turbulent, the accretion flow around rISCO is complicated. General relativistic MHD simulations in 3D are only beginning to shed light on the property of the black hole accretion flow (e.g., Beckwith, Hawley & Krolik 2008;Shafee et al. 2008;Noble, Krolik & Hawley 2008), and many uncertainties remain unresolved. Here, to make analytic progress, we adopt a simple viscous transonic flow model, which qualitatively describes the inner accretion flow of the black hole as long as the flow remains geometrically thin (see Afshordi & Paczynski 2003).

Boundary Condition at the Sonic Point
We rewrite equations (6)- (7) as where ′ stands for d/dr and we have used rΣur = constant for the background flow. Solving for δh ′ and δu ′ r we have Clearly, in order for the perturbation to be regular at the sonic point rs, where ur = −cs, we require A2 + csA1 = 0 at r = rs.
For definiteness, we characterize the variations of Σ and cs at the sonic point rs ≃ rISCO by the two length scales: From rΣur = constant, we also have u ′ r = cs(r −1 s + L −1 Σ ) at r = rs. Then equation (34) becomes This is the boundary condition for the fluid perturbations at the sonic point.

Properties of the Transonic Flow
Before exploring the effect the radial inflow on the disc modes, we first estimate the length scale for the surface density variation, LΣ, using the viscous slim disc model (e.g., Muchotrzeb & Paczynski 1982;Matsumoto et al. 1984;Abramowicz et al. 1988) The basic steady-state slim disc equations arė where l = r 2 Ω is the specific angular momentum of the flow, Ω is the actual rotation rate, ΩK is given by equation (2), νvis is the kinetic viscosity, and l0 is the eigenvalue that must be solved so that flow pass through the sonic point smoothly. We shall use the α-disc model, so that νvis = αHcs, with H ≃ cs/ΩK . To estimate LΣ, we assume l(r) ≃ r 2 ΩK (r) for r > ∼ rISCO and l0 ≃ lK(rISCO). Equation (39) gives valid for r > ∼ rISCO. At the radius r = rISCO + ∆r, we have ur(∆r/rISCO) 2 ≃ −4αHcs/rISCO. The sonic point (ur = −cs) is at ∆r ≃ 2 √ αHrISCO, and ur = −cs/2 at ∆r ≃ 2 √ 2αHrISCO. Thus u ′ r (rs) ∼ cs/LΣ, with where β = cs/(rΩ). This should be compared to the disc thickness H = βr: depending on the value of α, both LΣ < H and LΣ > H are possible. The value and sign of Lc depend on the thermodynamical and radiative properties of the flow, and cannot be estimated in a simple way. It is reasonable to expect |Lc| ∼ rs.

Reflectivity at the Sonic Point
We can understand qualitatively the effect of the transonic boundary condition on the p-mode by calculating the reflectivity Rs of the inner boundary.
Consider a density wave δh ∝ exp(i R r k dr) in the wave zone rin = rs < r < rIL, traveling toward the inner disc boundary.
To apply the boundary condition (36) to equation (42), we neglect ur in equations (6)-(8) at r = rs + ε, with ε ≪ rs and rs ≃ rISCO. Implicit in this procedure is the assumption that the fluid perturbations do not vary significantly between rs and rs + ε. We then obtain where and When considering the damping of the p-mode due to the transonic flow, the quantity |Rs| 2 − 1 is the most relevant (see section 5.4 below). Let we have Using β = cs/(rΩ),ω = −ωmΩISCO (where 0 <ω < 1), we find from equation (46) that (49) Figure 8 shows how the reflectivity depends on various parameters of the disc inner edge. In particular, for small LΣ/H = LΣ/(βrs), i.e., when the surface density of the disc decreases rapidly at the sonic point, |Rs| 2 is only slightly smaller than unity and the wave loss at the inner edge of the disc is small.

Mode Growth Rate in the WKB Approximation
Consider the p-mode trapped between rin = rs ≃ rISCO and rIL. With the reflectivity at rIL given by R (see section 3), we can write the wave amplitude for r < rIL as 4 δh ∝ " D rΣk On the other hand, with the reflectivity at rin = rs given by Rs, the wave can also be expressed as (42). For stationary waves we therefore require where Θr and Θi are real. The real eigen-frequency ωr is given by where RRs = |RRs| exp(iϕ), and n is an integer. The mode growth rate ωi is determined by |RRs| = exp(−2Θi), or For Θi = R r IL r in ki dr ≪ 1 and ki ≃ ωiωr/(cs √ω 2 r − κ 2 ), we obtain where we have assumedωr < 0. Equation (54) is to be compared with (19), where perfect reflection at rin is assumed. Clearly, to obtain growing modes we require |RRs| > 1. For a given |R| > 1, growing modes are possible only when the loss at the sonic point is sufficiently small (i.e., |Rs| is sufficiently close to unity).

Numerical Results
We solve equations (6)-(8) (with ur = 0) subjected to the radiative outer boundary condition (26) and the transonic inner boundary condition (36). Figure 9 depicts an example of the p-mode wavefunctions. Again, the discontinuity in the angular momentum flux F at rc signifies wave absorption; since rc < r peak , this leads to mode growth. Comparing with Fig. 4, here the angular momentum flux at rin is significantly nonzero, indicating wave loss through the sonic point. Nevertheless, the corotational instability is sufficiently strong to overcome the loss and makes the mode grow. Figures 10-11 show the fundamental p-mode frequencies and growth rates as a function of the disc parameters. Consistent with the result of section 5.3 (see Fig. 8), growing modes are obtained for sufficiently small LΣ. Negative Lc also tends to reduce wave loss at rs and make the growing modes possible. Such values of LΣ and Lc are not unreasonable for black hole accretion discs.
It is important to note that while the mode growth rates ωi depend sensitively on the inner disc parameters, particularly the physical property of the transonic flow near the ISCO, the real mode frequencies ωr show only weak dependence on the inner disc parameters (e.g., ωr decreases with increasing sound speed; see Fig. 6). Thus we may expect that kHz QPOs appear only in certain accretion states of the black hole, and the frequencies do not vary much as the accretion rate changes. Lovelace et al. (1999) (see also Li et al. 2000) have shown that when the vortensity ζ = κ 2 /(2ΩΣ) has an extremum at a certain radius (r peak ) in the disc 5 , it is possible to form normal Rossby modes around r peak . If the trapped Rossby waves can propagate on both sides of the corotation, a standing pattern of waves of opposite energies are formed, making the mode unstable -This is the "Rossby wave instability".  have considered the MHD version of the instability and suggested that it played a role in the diskoseismic modes around black holes (see also Tagger 2006).

THE ROLE OF ROSSBY WAVE INSTABILITY
We do not find any trapped Rossby modes in our calculation. To clarify the issue in light of works by Lovelace et al. and by Tagger   where we have used d ln ζ dr = r − rmin L 2 ζ , (for |r − rmin| ≪ rmin).
In this case, for a mode with ωr = mΩ(rmin) (or rc = rmin), we find V eff (rc) ≃ −2/(qL 2 ζ ) + 1/H 2 (for ωi ≪ ωr). When V eff (rc) < 0, or when Rossby waves can propagate on both sides of the corotation, leading to mode growth -this is the Rossby wave instability. Thus, the Rossby wave instability would operate if there existed a "sharp" vortensity minimum in the disc (with ζ varying on the lengthscale comparable or less than the disc thickness) -this is not the case for typical black hole accretion discs considered in this paper.

DISCUSSION
High-frequency QPOs (HFQPOs) in black-hole X-ray binaries have been studied observationally for more than a decade now and they provide a potentially important tool for studying the strong gravitational fields of black holes (see Remillard & McClintock 2006). Despite much theoretical effort, the physical mechanisms that generate these QPOs remain unclear (see section 1.1 for a brief review of existing theoretical models). Ultimately, numerical simulations of realistic accretion discs around black holes may provide the answer. However, such simulations are still at their early stage of development and have their own limitations (e.g., De Villiers & Hawley 2003;Machida & Matsumoto 2003Arras et al. 2006;Fragile et al. 2007;Reynolds & Miller 2008;Beckwith et al. 2008;Shafee et al. 2008;Noble et al. 2008), semi-analytical study remains a useful, complementary approach in order to identify the key physics involved.
In this paper, we have studied the global instability of the non-axisymmetric p-modes in black-hole accretion discs. These modes have frequencies ω ∼ (0.5 − 0.7)mΩISCO (where m is the azimuthal wave number, ΩISCO is the disc rotation frequency at the inner-most stable circular orbit), where the pre-factor (0.5-0.7) depends on the inner disc structure. Recent works Reynolds & Miller 2008;Fu & Lai 2008) suggested that, unlike other diskoseismic modes (g-modes and c-modes), the p-modes may be robust in the presence of disc magnetic fields and turbulence. Our linear analysis showed that due to GR effects, the p-modes may grow in amplitude due to wave absorptions at the corotation resonance. For a given m, only the lowest-order p-mode has sufficiently high frequency (ω > mΩ peak ; see Fig. 1) to be driven overstable by the corotational instability, while high-order (lower frequency) modes are damped by the corotational wave absorption.
The greatest uncertainty of our calculation of the p-mode growth rate concerns the boundary condition at the inner disc edge near the ISCO. In particular, the rapid radial inflow at the ISCO has the tendency to damp the mode (see Blaes 1986). While our analysis in section 5 indicates that this damping does not completely suppress the mode growth under certain disc conditions, it suggests that mode growth may not always be achieved in real black-hole accretion discs. Observationally, it is of interest to note that HFQPOs are observed only when the X-ray binaries are in the steep power-law state, while they do not appear in other spectral states (Remillard & McClintock 2006). In particular, HFQPOs are absent in the thermal (soft-high) state, believed to correspond to geometrically thin discs extending down to the ISCO. It is reasonable to expect that in this state p-modes are damped due to the rapid radial inflow.
Our current understanding of the steep power-law state (also called very high state) of black-hole X-ray binaries is rather limited. A thermal-radiation-emitting disc is suggested by spectral modelings, but it is not clear whether the disc is truncated at the ISCO or slightly larger radius (see Done et al. 2007). The observed power-law radiation component requires a significant corona that Compton up-scatters the disc thermal radiation. It is possible that in the steep power-law state, the inner disc behaves as a more reflective boundary (modeled in section 4) than a transonic flow (modeled in section 5), and thus more robust p-mode growth can be achieved. One possibility is that a significant magnetic field flux can accumulate in the inner disc when the disc accretion rate is sufficiently high (see Bisnovtyi-Kogan & Lovelace 2007;Rothstein & Lovelace 2008 and references therein). Such a magnetic field may also enhance the corotational instability and induce variability in the power-law radiation flux (see . Although the p-mode growth rates depend sensitively on a number of (uncertain) disc parameters (particularly those related to the inner disc boundary), the mode frequencies are more robust (see [10][11]. More precisely, the real mode frequency can be written as ωr =ωmΩISCO, whereω < 1 depends weakly on m and has only modest dependence on disc parameters (e.g. sound speed). This implies a commensurate frequency ratio as observed in HFQPOs (note that in some of our models, the m = 2, 3 modes have the largest growth rates; see Fig. 5). The fact thatω < 1 would also make the numerical values of the p-mode frequencies more compatible with the measurements of the QPO frequencies and black hole masses.
We note that our calculations in this work are done with a pseudo-Newtonian potential. For direct comparison with observations a fully general relativistic calculation 6 is needed including a careful treatment of the corotation singularity. Including the effect of black hole spin would likely increase the value of ω by modifying rISCO and ΩISCO, whileω will likely remain similar to the non-spinning case discussed above. We plan to study these effects in future work.