# Interface Modes and Their Instabilities in Accretion Disc Boundary Layers

## Abstract

We study global non-axisymmetric oscillation modes trapped near the inner boundary of an accretion disc. Observations indicate that some of the quasi-periodic oscillations (QPOs) observed in the luminosities of accreting compact objects (neutron stars, black holes and white dwarfs) are produced in the inner-most regions of accretion discs or boundary layers. Two simple models are considered in this paper: The magnetosphere-disc model consists of a thin Keplerian disc in contact with a uniformly rotating magnetosphere with and low plasma density, while the star-disc model involves a Keplerian disc terminated at the stellar atomosphere with high density and small density scale height. We find that the interface modes at the magnetosphere-disc boundary are generally unstable due to Rayleigh-Taylor and/or Kelvin-Helmholtz instabilities. However, differential rotation of the disc tends to suppress Rayleigh-Taylor instability and a sufficiently high disc sound speed (or temperature) is needed to overcome this suppression and to attain net mode growth. On the other hand, Kelvin-Helmholtz instability may be active at low disc sound speeds. We also find that the interface modes trapped at the boundary between a thin disc and an unmagnetized star do not suffer Rayleigh-Taylor or Kelvin-Helmholtz instability, but can become unstable due to wave leakage to large disc radii and, for sufficiently steep disc density distributions, due to wave absorption at the corotation resonance in the disc. The non-axisymmetric interface modes studied in this paper may be relevant to the high-frequency QPOs observed in some X-ray binaries and in cataclysmic variables.

## 1 Introduction

Quasi-periodic variabilities have been observed in the timing data of various types of accreting objects. Several types of quasi-periodic oscillations (QPOs) are observed in X-ray binaries with accreting black holes (BHs) or neutron stars (NSs) (e.g., Remillard & McClintock 2006; Van der Klis 2006). Oscillations are also seen in the outbursts of accreting white dwarf (WD) systems (e.g., Patterson 1981; see Warner 2004 for a review).

In accreting NS and BH X-ray binaries the observed QPO frequencies ( Hz for the high-frequency QPOs in the BH systems and Hz for kHz QPOs in the NS systems) imply a source close to the central compact object where the Keplerian orbital frequencies are high. Since the BH systems lack a hard surface where oscillations may occur, it is likely that the source of the variability is in the inner regions of the disc itself or in some interface regions between the disc and the plunging flow. Gilfanov et al (2003), however, found that, based on spectral analysis of the disc emission components, the quasi-periodic variability in Low Mass NS X-ray Binary systems are most likely caused by variations in the disc boundary layer, rather than the disc itself.

In Cataclysmic Variables (CVs) the Dwarf Nova Oscillations (DNOs) seen during during outbursts have frequencies roughly corresponding to the Keplerian rotation rate at the WD surface (e.g., Patterson, 1981; Warner, 2004; Knigge et al, 1998), which imply an origin at or near the inner disk boundary.

Several models involving accretion disk boundary dynamics have been proposed in different contexts. Popham (1999) studied the effect of a non-axysymmetric bulge at the optically thick to optically thin transition radius as a model for DNOs. Piro & Bildsten (2004) examined the surface wave oscillations that would occur within the thin equatorial belt around a non-magnetized WD formed by the accretion spreading layer, while Warner & Woudt (2002) considered accretion onto a slipping belt. In the context of accreting magnetic (neutron) stars, Arons & Lea (1976) and Elsner & Lamb (1977) considered the interchange instability at the magnetosphere boundary. Spruit & Taam (1990) and Spruit, Stehle & Papaloizou (1995) investigated the stability of thin rotating magnetized discs. Of particular relevance to the present paper is the work of Li & Narayan (2004), who examined a simplified cylindrical model of the Rayleigh-Taylor and Kelvin-Helmholtz instabilities at the boundary between a magnetosphere and an incompressible rotating flow. There have also been a number of numerical simulations of the interface at the magnetosphere-disc boundary (see Romanova et al., 2008 and Kulkarni & Romanova, 2008 and references therein).

In this paper we study global non-axisymmetric oscillation modes confined near inner boundary of the accretion disc (interface modes). We consider two simple models. The first model involves the magnetosphere-disc boundary similar to the model of Li & Narayan (2004): we consider an uniformly rotating incompressible magnetosphere with low gas density (where magnetic pressure dominates), which truncates a thin barotropic accretion disc (where gas pressure dominates). This situation may arise from magnetic field build up due to accretion (e.g., Bisnovatyi-Kogan & Ruzmaikin, 1974, 1976; Igumenshchev et al., 2003; Rothstein & Lovelace, 2008) or by the magnetosphere of a central (neutron) star. Unlike Li & Narayan (2004), who restricted their model to incompressible fluid, our discs are compressible and we show that because of the differential rotation of the disc, finite disc sound speed plays an important role in the development of the instability of the interface modes.

In our second model we examine the interface modes for accretion onto a non-magnetic stellar surface. Though the structure of the boundary layer is non-trivial and may affect boundary modes (see, e.g., Carroll et al., 1985; Collins et al., 2000), we consider the instabilities for a thin disc truncated by a sharp transition to a dense uniformly rotating stellar atmosphere. This simplified model may provide insight for modes with characteristic radial length scale much greater than the radial length scale of the boundary layer.

In Section 2 we describe the basic setup for the magnetospheric boundary model, and in Section 3 we discuss the resulting interface mode instabilities. We describe the star-disc boundary and analyze its possible instabilities in Section 4. We then conclude in Section 5 with a discussion of possible applications of our findings.

## 2 Magnetosphere-Disc Setup

We begin by considering a simplified model of the magnetosphere-disc boundary similar to the one considered by Li & Narayan (2004). The magnetic field is assumed to be negligible in the disc region (), while the magnetosphere region () is assumed to be incompressible and have low density compared to the disc region, with purely vertical magnetic field. Unlike Li & Narayan (2004), who assumed infinite sound speed in the disc, our disc has sound speed much less than the disk rotation speed .

In terms of the vertically integrated density (), pressure (), magnetic field () and fluid velocity () the ideal MHD equations are:

(1) | |||||

(2) | |||||

(3) |

where is the total pressure, is the magnetic tension, and is the gravitational potential due to the central object (e.g. Fu & Lai, 2008). Using cylindrical coordinates (), we consider the case where the magnetic field is purely poloidal and in the disc plane, which gives . We assume an axisymmetric background flow with fluid velocity . The unperturbed flow satisfies the condition

(4) |

The linearized equations of (1) and (2) with perturbations of the form (assuming no vertical dependence) take the form:

(5) | |||||

(6) | |||||

(7) |

where is the Doppler shifted frequency, is the radial epicyclic frequency, and , and are the Eulerian perturbations of the fluid variables. Additionally, assuming a barotropic flow we have

(8) |

with the sound speed .

### 2.1 The Magnetosphere

In the inner, magnetically dominated region (), we assume the flow to be incompressible, and have uniform rotation ( const) and uniform surface density ( const). Equations (5) - (7) then reduce to

(9) | |||||

(10) | |||||

(11) |

As in Li & Narayan (2004) we define and find and

(12) |

For uniform rotation, , equation (12) has the solution . Since , we take the positive sign to be the physical solution so that the perturbation falls off away from the interface. Thus the exact solution for the region is

(13) |

### 2.2 The Disc

In the disc (), we take the magnetic field to be small, such that , and the angular velocity of the unperturbed flow to be nearly Keplerian, such that . Rewriting equations (5) - (7) we have

(14) | |||||

(15) | |||||

(16) |

where

(17) |

is the enthalpy perturbation. Eliminating the velocity perturbations in favor of the enthalpy, we obtain the second order ODE for the enthalpy perturbation in the disc,

(18) |

where . For concreteness we will assume a power-law disc surface density profile .

### 2.3 Matching Conditions Across the Interface

The matching conditions across the interface at between the magnetosphere and the disc region are given by demanding the continuity of the Lagrangian displacement in the radial direction , and the total Lagrangian pressure perturbation across the boundary. The former gives

(19) |

where the subscript “” implifes that the quantities are evaluated at . The total Lagrangian pressure perturbation for is given by

(20) | |||||

In the disc region, we have

(21) | |||||

The condition then gives

(22) |

## 3 Interface Modes at the Magnetosphere-Disc Boundary

Perturbations mainly confined to the magnetosphere-disc interface can become unstable due to Rayleigh-Taylor or Kelvin Helmholtz instability. In order to calculate the growth rates, we must solve the eigenvalue problem given by equation (18) with an outgoing wave boundary condition at some outer radius, and equation (22) at the interface radius .

### 3.1 Numerical Solution

We adopt the radiative outer boundary condition in the outer wave zone of the disc, such that far from the outer Lindblad resonance radius (where ) we have the solution of the form:

(23) |

with and (see Tsang & Lai 2008, Lai & Tsang, 2008). This gives the boundary condition at :

(24) |

We adopt (22) as the inner boundary condition for the disc and solve the eigenvalue problem using a standard shooting method (Press et al 1998). For the numerical solutions below, the density profile of the disk was assumed to be so that corotation absorption plays no role in determining the mode stability (Tsang & Lai 2008). An example wavefunction for an interface mode is shown in Figure 1, for typical disc parameters.

The numerical eigenvalues are shown in Figure 2 and Figure 3 for various disc and magnetosphere parameters, for .

### 3.2 Discussion of Numerical Results

Figure 2 shows the complex eigenvalues as a function of sound speed , for density contrasts corresponding to and , with magnetosphere rotation rate equal to the Kepler frequency at the interface []. For this case we see that there exists a cutoff in the disc sound speed below which no growing interface modes are found. This arises from the stabilizing effect of the background differential rotation, and can be understood as follows.

Setting and rewriting (22) in terms of the radial velocity perturbation , we have:

(25) |

where [see Eq. (4)] and is the vorticity, and where all quantities are evaluated at the interface . For the wave frequencies of interest, the waves are evanescent in the region of the disk just outside the interface. Let . Equation (25) can be solved in terms of , giving

(26) |

where . The terms inside the square root correspond to the mode growth due to Rayleigh-Taylor instability and the suppression due to vorticity, respectively. With we find the critical sound speed

(27) |

above which the perturbations will be unstable.

Figure 3 shows cases where the inner region is uniformly rotating at an angular frequency of one half the Kepler frequency at the interface []. When this leads to the development of the Kelvin-Helmholtz instability, and both this and the Rayleigh-Taylor instability play a role in the mode growth. In the Appendix we derive the expression for the plane-parallel Rayleigh-Taylor and Kelvin-Helmholtz instabilities for a compressible upper region (with density and horizontal velocity ), and incompressible lower region (with density and horizontal velocity ). For we have where is the horizontal wavenumber and

(28) |

Here , g is the acceleration due to gravity in the vertical direction, and is the vertical scale height in the upper region.The Kelvin-Helmholtz term is approximately

(29) |

For this reduces to the incompressible limit with . For we have , a factor of smaller than the incompressible result.

For the rotating system under consideration, the imaginary part of the mode frequency can be written schematically as [cf. equation (28)]

(30) |

We also have and so . Thus , and depends weakly on sound speed. On the other hand, from equation (26) we see that the vorticity suppresses mode growth through the term . For sufficiently small , equation (18) indicates , i.e. . Therefore the vorticity term scales with sound speed as , and can be dominated by the Kelvin-Helmholtz term for small enough sound speed. In the left panel of Figure 3, the mode growth () for small is mainly driven by the Kelvin-Helmholtz instability. For the sound speed ranges where and dominate over overlap, and hence the critical sound speed in equation (27) is not relevant. For larger values of these regions can overlap for all .

### 3.3 Effect of a Relativistic Potential

While in Section 3.2 and other sections of the paper we focus on Newtonian discs, it is of interest to consider how general relativity may modify our results. The effect of general relativity can be approximated by using the pseudo-Newtonian Paczynski & Wiita (1980) potential:

(31) |

with the Schwarzschild radius. This gives the Keplerian orbital frequency () and epicyclic frequency () as:

(32) |

with at .

For the interface modes are the same as for the Newtonian case. However, as the suppression effect of [see Eq. (26)] is reduced as the vorticity goes to zero at , so that for there is no cutoff sound speed (see Figures 2-3) for interface mode instability. This is illustrated in Figure 4. Thus if the magnetosphere boundary is at , the interface modes will always be present and highly unstable for any sound speed.

### 3.4 P-modes with Magnetosphere Boundary

The boundary condition given by equation (22) also provides an inner reflection boundary for disc p-modes, which were studied in detail in Lai & Tsang (2008). These modes have wavefunctions primarily “trapped” in the wave region between the disc boundary and the inner Lindblad resonance radius, , where . Figure 5 depicts an example of the p-mode wave function for the same disc model as in Figure 1. The growth rates of these p-modes are determined primarily by the outgoing flux at the outer boundary and the effect of the corotation resonance, as discussed in Lai & Tsang (2008). In Figure 6 the eigenfrequencies are shown for p-modes in a disc with the density profile , where so that wave absorption at the corotation resonance is inactive (since in this case the vortensity is constant). For the density profile , the corotation absorption tends to damp the mode, while for the corotation absorption enhances it.

## 4 Interface Modes at the Star-Disc Boundary

### 4.1 Star-Disc Boundary Condition

In the case of accretion on to a non-magnetic star, our model consists of a dense uniformly rotating compressible stellar atmosphere truncating the accretion disc. This model ignores the structure of the boundary layer. However the qualitative properties of the dynamics should be captured for modes with characteristic radial length scale much greater than the radial scale of the boundary layer.

Several studies of CVs (e.g. by examining the rotationally broadened line emissions from the stellar surface) have shown that the stellar rotation rates are significantly below the breakup rotation rate (see Warner, 2004), and we limit our examinations to systems with . As in the case of the disc, we consider only the effect of perturbations on a cylindrical equatorial surface of the stellar atmosphere (i.e., we are considering a ”cylindrical” star). In this region equation (18) also describes the enthalpy perturbations within the stellar atmosphere. Inside the atmosphere, we assume a small constant density scale height, . Equation (18) then becomes

(33) |

For and , this has the solution

(34) |

The Lagrangian pressure perturbation at the stellar surface is then

(35) | |||||

Once again matching the Lagrangian displacement and pressure perturbation at the interface gives the boundary condition for the interface modes for the star-disc boundary case:

(36) |

### 4.2 Numerical Results

We repeat the numerical procedure of Section 3 using the radiative outer boundary condition [equation (24)] and using equation (36) as the inner disc boundary condition. A sample wavefunction for the star-disc interface mode is shown in Figure 7. For typical mode frequencies the region just outside the boundary is an evanescent zone; wave propagation becomes possible only beyond the outer Lindblad resonance (). Figure 8 shows the eigenfrequencies for the lowest order modes with as a function of disc sound speed, for representative parameters , and . Figure 9 shows the dependence of the mode eigenfrequencies on the density (), rotation rate (), and scale height () of the star. We see that both the real mode frequency and the growth rate do not depend strongly on these parameters.

### 4.3 Discussion of Numerical Results

When the disc is truncated by the star’s surface, the effective gravity acts to stabilize the perturbations (since ), as does the vorticity. Thus compared to the interface mode in the magnetosphere-disc case (Sections 3.1 – 3.2), the mode growth rates here are much smaller and are primarily driven by wave propagation through the corotation, beyond the outer Lindblad resonance. In the left panels of Figure 8 the eigenfrequencies are shown for typical parameters [, , ], and disc density index , so that the corotation absorption plays no role. For other density indices, wave absorption at the corotation can act to either damp or grow the interface modes (Tsang & Lai 2008; Lai & Tsang 2008). For example, in the Shakura-Sunyaev -disc model the disc solution for the outer disc solution (with free-free opacity and gas pressure dominating) has the surface density , hence the modes would be stabilized by absorption at the corotation resonance. However, for models where the disc has density index at corotation, the corotational absorption acts to enhance mode growth, as shown in the right panels of Figure 8.

The mode eigenfrequencies have very little dependence on the properties of the stellar atmosphere (), as shown in Figure 9. The mode frequencies instead primarily depend on the disc sound speed, which in turn depends on the accretion rate. Observations of CVs indicate that DNOs are usually only detected in high states, with the oscillation period decreasing with increasing luminosity (Warner 2004). The downward trend of the as a function of in Figure 8 would appear to contradict the observed period-luminosity anti-correlation. But note that in our model, the interface mode frequency depends on the sound speed at the inner-most disc region and boundary layer, and it will be necessary to model the thermodynamic and radiative properties of the boundary layer in order to compare with observation directly. Also, the oscillations of the type considered here would yield periods shorter than the surface Keplerian period, except for the mode. Though higher- modes would be more difficult to observe due to the averaging out of the luminosity variation, most observed DNOs, even those with 1:2:3 harmonic structure (Warner & Woudt, 2005) occur with period roughly at or greater than the corresponding surface Keplerian period. These long-period oscillations cannot be explained by the model considered here.

## 5 Conclusions

We have studied the non-radial oscillation modes at the interface between an accretion disc and a magnetosphere or stellar surface. Although the models explored in this paper are perhaps too simplified compared to realistic situations, they offer some insight into the behavior of the interface modes in various astrophysical contexts.

Our study of the interface modes at the magnetosphere-disc boundary extended the work by Li & Narayan (2004), who considered incompressible disc flow (and therefore could not treat real discs). The model can have very strongly unstable modes due to Rayleigh-Taylor and Kelvin-Helmholtz instabilities. In systems where the magnetosphere has developed from advection of frozen magnetic flux, the magnetosphere is expected to be roughly rotating with at the Keplerian rate. Since there is no shear at the interface, only the Rayleigh-Taylor instability may occur. However the disc vorticty (due to differential rotation) acts to suppress the instability, leading to a cutoff below a critical sound speed. Thus a sufficiently hot disc is required to generate unstable low- modes. For magnetospheres rotating with the central star, shearing is expected between the magnetosphere and disc, and the Kelvin-Helmholtz instability becomes active. This can help to drive the instability for low- modes to overcome the vorticity in low sound-speed discs. In discs that terminate near the ISCO in a general relativistic potential, the vorticity approaches zero at the inner disc radius, and unstable interface modes can be found for any sound speed.

We can expect a strong dependence of the interface mode growth rates on the sound speed, and thus accretion rate, while the real frequencies of these oscillations remain close to , and do not depend strongly on the sound speed. It is worth noting that the same boundary condition that gives rise to the interface modes also gives rise to intertial-acoustic modes (or p-modes) in the disc.

Although higher- interface modes are more unstable in our model, these are less likely to be observed due to the averaging out of the luminosity variation over the observable emitting region. In addition, if the effect of viscous damping is considered (e.g., Wang & Robertson 1985), small wavelength or high-, perturbations are suppressed.

These results are for perturbations with no vertical structure, and are applicable mainly to the midplane of accretion discs interacting with a magnetosphere. Global 3D numerical studies of Rayleigh-Taylor instability induced accretion onto magnetized stars have been performed by Romanova, Kulkarni & Lovelace (2008) and Kulkarni & Romanova (2008), and show such small instabilities in the disk midplane. The low- oscillations at the magnetosphere-disc interface may be relevent to the high-frequency QPOs observed in some NS and BH X-ray binary systems (Li & Narayan, 2004; see Section 1 of Lai & Tsang for a critical review of various theoretical models), although to obtain the correct QPO frequencies for the BH systems, the disc inner radius must lie outside the inner-most stable circular orbit.

For the star-disc boundary model considered in in this paper, the interface mode growth rates are much smaller than for the magnetospheric case, since the effective gravity now acts to stabilize the system, and the Rayleigh-Taylor instability is inactive. The modes discussed here are unstable due primarily to propagation through the corotation. With sufficiently steep disc density profile ( with ), corotation absorption can also help to drive these modes, as studed previously by Tsang & Lai (2008) and Lai & Tsang (2008). Such modes may be responsible for the high-frequency (of order the Keplerian frequency at the stellar surface) dwarf nova oscillations observed in CVs, although oscillations with longer periods would require a different explanation.

## Acknowledgements

This work has been supported in part by NASA Grant NNX07AG81G, NSF Grant AST 0707628, and by Chandra Grant TM6-7004X (Smithsonian Astrophysical Observatory).

## Appendix: Plane Parallel Flow with a Compressible Upper Layer

Consider a system consisting of two fluids in the gravitational field . The upper fluid has density (with , where is the sound speed) and horizontal velocity along the x-axis; the lower fluid is incompressible with density and horizontal velocity .

The linear perturbation equations for the upper fluid are

(37) | |||||

(38) |

where . For perturbations of the form , these become

(39) | |||||

(40) | |||||

(41) |

where . Assuming the perturbation is isothermal, so that , we obtain

(42) |

The two independent solutions of equation (42)

(43) |

where . Obviously, the physically relevant solution is

(44) |

with . This gives, by equation (40) the Eulerian pressure perturbation for the upper region:

(45) |

For the lower region the fluid is an incompressible potential flow with with satisifying . For z, the appropriate solution is

(46) |

This gives . The Eulerian pressure perturbation in the lower region () is then:

(47) |

Matching the Lagrangian displacement and Lagrangian pressure perturbation across the boundary between the upper and lower regions we get

(48) |

which can be written as the quadratic:

(49) |

where has a non trivial dependence on . This has a solution

(50) |

where .

### Footnotes

- footnotemark:

### References

- Arons, J. & Lea, S. M., 1976, ApJ, 207, 914
- Bisnovatyi-Kogan, G. S. & Ruzmaikin, A. A., 1974, Ap&SS, 28, 45
- Bisnovatyi-Kogan, G. S. & Rusmaikin, A. A., 1976, Ap&SS, 42, 401
- Carroll, B. W., Cabot, W., McDermott, P. N., Savedoff, M. P. & Van Horn, H. M., 1985, ApJ, 296, 529
- Collins, T. J. B., Helfer, H. L. & Van Horn, H. M., 2000, ApJ, 534, 944
- Elsner, R. & Lamb, F. K., 1977, ApJ, 215, 897
- Fu, W., Lai, D. 2008, ApJ, 690, 1386
- Ghosh, P. & Lamb, F. K., 1978, ApJ, 223, L83
- Ghosh, P. & Lamb, F. K., 1979, ApJ, 232, 259
- Gilfanov, M., Revnivtsev, M. & Molkov, S., 2003, A&A, 410, 217
- Igumenshchev, I. V., Narayan, R. & Abramowicz, M. A., 2003, ApJ, 592, 1042
- Knigge, C., Drake, N., Long, K. S., Wade, R. A., Horne, K. & Baptista, R., 1998, ApJ, 499, 29
- Kulkarni, A.K. & Romanova, M.M., 2008, MNRAS, 386, 673
- Lai, D., & Tsang, D., 2008, MNRAS, in press (arXiv:0810.0203)
- Li, L., Narayan, R., 2004, ApJ, 601, 414
- Ortega-Rodriguez, M. & Wagoner, R. V., 2007, ApJ, 668, 1158
- Paczynski, B., Wiita, P.J. 1980, A&A, 88, 23
- Patterson, J., 1981 ApJS, 45, 517
- Piro, A. L. & Bildsten, L., 2004, ApJ, 616, L155
- Popham R., 1999, MNRAS, 308, 979
- Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 1998, Numerical Recipes (Cambridge Univ. Press)
- Remillard, R. A. & McClintock, J. E., 2006, ARAA, Vol. 44, pp. 49-92
- Romanova, M. M., Kulkarni, A.K. & Lovelace, R.V.E, 2008, ApJ, 673, L171
- Rothstein, D.M. & Lovelace, R.V.E., 2008, ApJ, 677, 1221
- Spruit, H. C. & Taam, R. E., 1990, A&A, 229, 475
- Spruit, H.C., Stehle, R. & Papaloizou, J.C.B., 1995, MNRAS, 275,1223
- Tsang, D. & Lai, D., 2008, MNRAS, 387, 446
- van der Klis, M. 2006, in Compact Stellar X-ray Sources, ed. W.H.G. Lewin and M. van der Klis (Cambridge Univ. Press) (astro-ph/0410551)
- Warner, B., 2004, PASP, 116, 115
- Warner, B. & Woudt, P. A., 2002, MNRAS, 335, 84
- Warner, B. & Woudt, P.A., 2005, The Astrophysics of Cataclysmic Variables and Related Objects, Proceedings of ASP Conference Vol. 330. Edited by J.-M. Hameury and J.-P. Lasota. San Francisco: Astronomical Society of the Pacific, 2005., p.227 (arXiv:astro-ph/0409287)