WAVE DYNAMICS OF ACCRETION DISKS AND OTHER ROTATING ASTROPHYSICAL FLOWS
A Dissertation Presented to the Faculty of the Graduate School
of Cornell University in Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy
by Wen Fu May 2012

c 2012 Wen Fu ALL RIGHTS RESERVED

WAVE DYNAMICS OF ACCRETION DISKS AND OTHER ROTATING ASTROPHYSICAL FLOWS Wen Fu, Ph.D. Cornell University 2012
A general consensus on the origin of quasi-periodic oscillations observed in accreting X-ray binaries is that they are related to the wave dynamics of the accretion disk. We conduct the first systematic study on the detailed dynamical properties of different discoseismic modes, in particular the effects of magnetic fields. We show through local analysis that even a weak magnetic field can “destroy” the self-trapping zone of inertial-gravity modes. The so-called corrugation modes are also strongly affected when the poloidal field approaches equal-partition. Whereas the basic wave properties of inertial-acoustic modes (p-modes) are not affected qualitatively by disk B fields. These modes become unstable large-scale oscillations when the disk vortensity (vorticity divided by surface density) profile has positive gradient at the corotation (where wave pattern speed matches background flow rotating speed) and they are not qualitatively affected by inner disk boundary as long as it has some reflectivity. With numerical simulation, we probe the nonlinear evolution of global overstable disk p-modes and demonstrate that they are quite robust even after non-linear saturation. We find, however, that disc B field can split the corotation resonance and significantly reduce the growth rates of these modes, thus challenging its viability as a model for observed high frequency quasi-periodic oscillations.
We employ similar technique to study the dynamics of other astrophysical flows, such as accreting tori and rotating protoneutron stars. Similar suppress-

ing effects from magnetic fields are also found in thin accreting tori and rotating protoneutron stars. We found that magnetic fields tend to suppress the Papaloizou-Pringle instability in relatively thin tori, while it can lead to a new instability in relatively thick tori where Papaloizou-Pringle instability cannot operate if the system is non-magnetized. Differentially rotating neutron stars have long been known to be subject to a global non-axisymmetric instability. We show that a toroidal stellar magnetic field can suppress this instability when the B field is a few × 1016 G or stronger.

BIOGRAPHICAL SKETCH
Wen Fu was born in 1986 and grew up in a small town along the Yangtze river. In his childhood, he liked to go fishing and swimming in the unbearably hot summer of central China, while in the frigid winter he enjoyed collecting leafs, hay, branches and starting a “bonfire” with his cousin in the middle of the dry paddy field. He started school early at the age of five because at that time the government policy on starting school age had not been implemented in his hometown yet, and his parents thought this could help save extra daycare expense. Under the influence of his father who was a history teacher and his big sister who majored in film study, Wen intended to study social science in high school, but ended up studying natural science after he was picked for a class of special talented – a class that was formed for the national math and physics Olympic competitions. He fell in love with astronomy in college and managed to get into Cornell for graduate study when he was twenty. The life in upstate New York proved to be a bittersweet experience for him. On one hand, Wen was fascinated by the gorgeous summer of Ithaca with all the beautiful trees and flowers on campus, clear and crisp summer sky, water tinkering down the gorges, birds chirping and deer wandering outside the living room; on the other hand, he was not so fond of the harsh, depressing Ithaca winter, with the bleak sky, snow storms, slippery and muddy roads, etc.. Wen met his love of life – Xu Li back in college, but they ended up in different cities for graduate school. During his five years at Cornell, Wen made a sizeable “contribution” to the US airline industry by frequently flying back and forth between Ithaca and Houston. Fortunately, Wen got a job offer from Rice University and will embark upon a new journey with his wife in Houston.
iii

To my parents. iv

ACKNOWLEDGEMENTS
First and foremost, I would like to thank my advisor Dong Lai. It has been a great pleasure to be one of his students. I appreciate all his contributions of time, ideas and funding to make my Ph.D. experience productive. I am grateful for his inspiration, motivation and encouragement, even during my tough days. I am also grateful for the wonderful example he has provided as a good astrophysicist and professor. Moreover, I want to thank him for his understanding of my two body issue and his efforts to accommodate my commuting needs. I do not think I could have made it through if it was not for his kind support.
I am greatly indebted to my committee members: David Chernoff, Richard Lovelace and Phil Nicholson for guiding me through my Ph.D. pursuit and for generously helping with all the paperworks. Phil, with Terry Herter had been very nice and really helped me out when I could not pass the graduate school English oral proficiency test for my first semester and had trouble keeping my TA job. I am also appreciative of all the assistance from our department staff members: Daniel O’Connor, Monica Armstrong, Mary Mulvanerton, Sherry Falletta. They are the most friendly and professional staff members that I have ever seen.
The other members of our group and my office mates throughout the years have been a wonderful source of friendship as well as good advice and collaboration. I want to thank them for the stimulating discussions and the fun they brought into the office life: Zach Medin, David Tsang, Chen Wang, James Fuller, Natalia Shabaltas, Francois Foucart, David Bernat, Nick Taylor, Abdul Mroue, Ruxandra Bondarescu, Richard Kipphorn, Nishant Agarwal. I also very much appreciated the valuable friendship from Rong Long, Yeyun Zhou, Mengqiao Wang, Bo Li, Heng Li, Xiang Ma, Shan Huang, Yiming Li, Yanling Wu, Lei Hao,
v

Min Long. A part of my thesis work was carried out at Los Alamos National Laboratory.
I would like to thank Hui Li for hosting me and helping me with everything at Los Alamos. I gratefully acknowledge the support from a Laboratory Directed Research and Development Program at LANL. During my Ph.D. study, I have also been supported by the these grants: NSF Grant AST-0707628, AST-1008245 and NASA Grant NNX07AG81G, NNX10AP19G.
Lastly, I want to thank my family for all their love and support. For my parents who raised me up and worked so hard to support me all the way to college. For my sister whose encouragements were always with me during the difficult times. And for my loving, supportive, and patient wife Xu Li who accompanied me throughout this memorable journey.
vi

TABLE OF CONTENTS

Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

iii iv v vii x xi

1 Introduction 1.1 Quasi-Periodic X-Ray Oscillations and discoseismic modes . . . . 1.2 Dynamics of a Magnetosphere-Disc System . . . . . . . . . . . . . 1.3 Dynamics of magnetized accretion tori . . . . . . . . . . . . . . . . 1.4 Dynamics of rotating protoneutron stars . . . . . . . . . . . . . . .

1 1 9 13 16

2 Effects of Magnetic Fields on the Discoseismic Modes of Accreting

Black Holes

20

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.2 Basic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.3 Hydrodynamic Limit: Discoseismic Modes . . . . . . . . . . . . . 23

2.4 Effect of Poloidal Fields . . . . . . . . . . . . . . . . . . . . . . . . 27

2.4.1 P-modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.4.2 G-modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.4.3 C-modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.5 Effect of Toroidal Fields . . . . . . . . . . . . . . . . . . . . . . . . 33

2.5.1 P-modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.5.2 G-modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3 Corotational Instability, Magnetic Resonances and Global Inertial-

Acoustic Oscillations in Magnetized Black-Hole Accretion Discs

38

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.2 Setup and Basic Equations . . . . . . . . . . . . . . . . . . . . . . . 39

3.3 Effective potential and Wave Propagation diagram . . . . . . . . . 46

3.3.1 Hydrodynamical discs . . . . . . . . . . . . . . . . . . . . . 47

3.3.2 Discs with pure toroidal magnetic fields . . . . . . . . . . . 49

3.3.3 Discs with pure vertical magnetic fields . . . . . . . . . . . 53

3.3.4 Discs with mixed magnetic fields . . . . . . . . . . . . . . . 55

3.4 Super-Reflection of the Corotation Barrier . . . . . . . . . . . . . . 55

3.4.1 Non-magnetic Discs . . . . . . . . . . . . . . . . . . . . . . 56

3.4.2 Numerical calculations of reflectivity . . . . . . . . . . . . 57

3.5 Global P-Modes of Black-Hole Accretion Discs . . . . . . . . . . . 59

3.5.1 Results for discs with pure toroidal magnetic fields . . . . 61

3.5.2 Results for discs with mixed magnetic fields . . . . . . . . 66

3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

vii

4 Corotational Instability in Black-Hole Accretion Discs: Nonlinear

Simulations

70

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.2 Numerical Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5 Dynamics of the Innermost Accretion Flows Around Compact Objects: Magnetosphere-Disc Interface, Global Oscillations and Instabilities 82 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.2 Equilibrium and Perturbation equations . . . . . . . . . . . . . . . 83 5.2.1 The magnetosphere . . . . . . . . . . . . . . . . . . . . . . . 85 5.2.2 The disc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.2.3 The interface . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.3 Interface modes: Newtonian potential . . . . . . . . . . . . . . . . 90 5.3.1 Numerical solutions . . . . . . . . . . . . . . . . . . . . . . 91 5.3.2 Discussion of results . . . . . . . . . . . . . . . . . . . . . . 94 5.4 Interface modes and Discoseismic modes: Pseudo-Newtonian potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.4.1 Non-magnetic discs . . . . . . . . . . . . . . . . . . . . . . 99 5.4.2 Magnetic discs . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.4.3 Effects of different inner disc radii . . . . . . . . . . . . . . 105 5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.5.1 Interface instabilities in a rotating, magnetized system . . 108 5.5.2 Inertial-Acoustic Modes of relativistic Disc . . . . . . . . . 109

6 Papaloizou-Pringle Instability of Magnetized Accretion Tori

111

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

6.2 Equilibrium Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

6.2.1 Model (a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

6.2.2 Model (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

6.3 MHD Equations for Perturbations . . . . . . . . . . . . . . . . . . 118

6.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 120

6.4.1 Model (a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

6.4.2 Model (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

6.5 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

6.5.1 Model (a): Pure Toroidal Field Configuration . . . . . . . . 129

6.5.2 Other Magnetic Field Configurations . . . . . . . . . . . . 132

6.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

7 Low-T/|W| instabilities in differentially rotating proto-neutron stars

with magnetic fields

136

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

7.2 Equilibrium Model of A Magnetized Rotating Cylinder . . . . . . 137

7.3 Linear Perturbation Analysis . . . . . . . . . . . . . . . . . . . . . 140

viii

7.3.1 Perturbation equations . . . . . . . . . . . . . . . . . . . . . 140 7.3.2 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . 143 7.3.3 Cowling approximation . . . . . . . . . . . . . . . . . . . . 145 7.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 7.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

8 Conclusion

156

A General Magnetic Field profiles in the Magnetosphere

162

B Plane-parallel flow with a compressible upper layer and a magnetized

lower layer

164

Bibliography

166

ix

LIST OF TABLES 4.1 Comparison of results from linear and nonlinear studies of un-
stable disk p-mode . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
x

LIST OF FIGURES

1.1 High-frequency quasi-periodic oscillations observed in blackhole binary and black-hole candidate systems. The traces in blue show power density spectrums (PDSs) for the range 1330 keV. Red traces indicate PDSs with a broader energy range, which may be either 2–30 keV or 6–30 keV (Image taken from Remillard & McClintock 2006). . . . . . . . . . . . . . . . . . . . . . . .

2

2.1 Wave propagation diagram showing various trapped modes in BH accretion disks: (a) axisymmetric p-mode; (b) axisymmetric g-mode; (c)non-axisymmetric p-mode; (d) non-axisymmetric gmod√e and c-mode. The cur√ves depict various critical frequencies (κ, nΩ⊥, Ω, Ω ± κ/m, Ω ± nΩ⊥/m), the vertical dotted lines denote the inner-most stable circular orbit (ISCO). The curvy horizontal lines specify wave propagation zones and the height of the line is ω (for the top panels) or ω/m (for the bottom panels) of the mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 The effect of poloidal magnetic field on the g-mode propagation zone for m = 0, η = 1. In each panel, the upper three curves are ω2 (eq. [2.32]) and the lower two curves are ω3 (eq. [2.33]). Axisymmetric g-modes of frequency ω can propagate in the region where ω3 < ω < ω2. The solid line refers to the case of b = 0, the short-dashed line b = 0.03 and the long-dashed line b = 0.1,
where b ≡ vAz/cs, with vAz = Bz/ 4πρ (Alfve´n speed) and cs the sound speed. The vertical dotted lines correspond to the inner disk radius at ISCO. The curvy horizontal lines specify the wave propagation zones, and the height of the line is ω of the mode. The lower and upper panels are for the case of a Schwarzschild BH (a = 0) and a Kerr BH (a = 0.8M), respectively. The angular frequencies are in units of M−1 = c3/(GM). . . . . . . . . . . . . . . 29

xi

2.3 The effect of poloidal magnetic field on the g-mode propagation zone for m 0, η = 1, a = 0. The three panels are for b = vAz/cs = 0, 0.005 and 0.015. The solid line, short-dashed lines and the long-dashed lines show Ω, Ω ± ω2 and Ω ± ω2/2, respectively. In the bottom panel, the dotted lines show Ω ± ω3 (In the upper and middle panels, Ω ± ω3 almost concide with Ω, since ω3 = 0 for b = 0 and ω3 ≪ Ω for b ≪ 1). Non-axisymmetric g-modes can propagate in the region where Ω − ω2/m < ω/m < Ω − ω3/m or Ω + ω3/m < ω/m < Ω + ω2/m. The vertical dotted lines correspond to the inner disk radius at ISCO. The curvy horizontal lines in top and middle panels specify wave propagation zones and the height of the line is ω/m of the mode. Note that the self-trapping zone (depicted in the upper and middle panels) disappears as b increases. The angular frequencies are in units of M−1 = c3/(GM). 31
2.4 The effect of poloidal magnetic field on c-modes with m = 1 and a = 0.2M. The upper panel shows the original c-modes since b is small and ω1 = Ω⊥; in the bottom panel, with a large b, the c-mode trapping zone is instead bounded by the inner reflection boudary and Ω − ω2. The vertical dotted line refers to the inner disk radius at ISCO. The curvy horizontal lines specify wave propagation zones and the height of the line is ω of the mode. . . 33
2.5 Wave propagation diagram for non-axisymmetric g-modes and c-modes, for b = vAz/cs = 0.7, m = 2, η = 1, and a = 0. The solid line, dot-short dashed lines, long-dashed lines, dotted lines and short-dashed lines show Ω, Ω ± ω4/2, Ω ± ω3/2, Ω ± ω1/2 and Ω ± ω2/2 (eqs. [2.31]-[2.34]), respectively. The vertical dotted line shows the inner disc radius. The curvy horizontal lines specify wave propagation zones and the height of the line is ω/2 of the mode. The singular points (where kr → ∞) are indicated by filled circles. The angular frequencies are in units of M−1 = c3/(GM). . . 34
xii

3.1 Wave propagation diagram in accretion discs when the background fluid vortensity gradient (evaluated at the corotation resonance radius) (dζ/dr)c = 0. In the upper panel, the solid line shows the effective potential V(r) as a function of r for B = 0 discs, and the dashed line shows V(r) for discs with finite toroidal magnetic fields. Waves can propagate only in the region where V(r) < 0. The labels ILR, OLR and CR denote the inner Lindblad resonance, outer Lindblad resonance and corotation resonance, respectively. The lower panel shows the blow-up of the region near the CR (located at rc = 1), and the long-dashed line shows the case with a stronger Bφ than the short-dashed line. The labels IMR and OMR denote the inner magnetic (slow) resonance and outer magnetic resonance, respectively. Note that the vertical scales in the upper and bottom panels differ by a large factor. A Newtonian potential is used in calculating V(r). (A color version of this figure is available in the online journal.) . . . 48
3.2 Same as Fig. 3.1, except for the case (dζ/dr)c < 0. . . . . . . . . . . 49 3.3 Same as Fig. 3.1, except for the case (dζ/dr)c > 0. . . . . . . . . . . 50 3.4 Numerical values of f2+ and f2− [see Eq. (3.48)] as a function
of the dimensionless ratio vAφ/(rΩ) (evaluated at r = rc), which specifies the toroidal magnetic field strength in the disc. Dimensionless units are adopted, where G = M = rc = 1, m = 2, and cs/(rΩ) ≪ 1. Different line types represent different disc density and magnetic field profiles. . . . . . . . . . . . . . . . . . . . . . . 51 3.5 Effective potential around the corotation/magnetic resonance region for a disc with mixed magnetic fields (solid line). The dotted line shows V(r) for the case with a pure toroidal field. The corotation resonance (CR), inner/outer magnetic (slow) resonances (IMR and OMR) and inner/outer Alfve´n resonances (IAR and OAR) are indicated. Note that the MRs represent secondorder singularities in V(r), while the ARs (which exist only for the mixed field case) are first-order singularities. . . . . . . . . . . 54 3.6 Reflection coefficient as a function of ν for a Newtonian Keplerian disk. Different lines represent different magnetic field strengths, measured by the dimensionless parameters bφ = (vAφ/rΩ)c and bz = (vAz/rΩ)c. The other parameters are m = 2, Bφ ∝ Bz ∝ rq with q = 0 and β = cs/(rΩ) = 0.1. The parameter ν is defined according to Eq. (3.60). . . . . . . . . . . . . . . . . . . . . 58
xiii

3.7 Example wavefunctions of m = 2 p-modes in BH accretion discs. The disc density profile is ρ ∝ r−1 (i.e., σ = 1) and the sound speed is cs = 0.1rΩ (i.e., β = 0.1). The left columns are for a B = 0 disc, with the eigenvalue ω = (0.9324 + i0.0027)Ωin (where Ωin is the disc rotation rate at r = rin = rISCO); the right panels are for a disc with bφ = (vAφ/rΩ)in = 0.01 [corresponding to (vAφ/rΩ)c = 0.0175] (with Bφ independent of r, i.e., q = 0), with ω = (0.9312 + i0.0018)Ωin. In each subfigure, the upper panels show different variables as a function of r in the full fluid zone. The solid and dotted lines depict the real and imaginary parts, respectively. The bottom panels zoom in the region near the corotation resonance, with the vertical lines indicating the locations of the CR (left) and magnetic slow resonances (right). The vertical scales of the wavefunctions are arbitrary, with δh(rout) = 1. 62
3.8 Angular momentum flux carried by the wave as a function of r. The disc parameters are σ = 1, q = 0 and β = 0.1, and the wave modes have m = 2. Panel (a) shows the result of a B = 0 disc [with the mode frequency ω = (0.9324 + i0.0027)Ωin] with the vertical line indicating the corotation resonance. Panel (b) is for a magnetized disc with bφ = 0.01 [corresponding to (vAφ/rΩ)c = 0.0175] [left column; with ω = (0.9312 + i0.0018)Ωin] or bφ = 0.015 [corresponding to (vAφ/rΩ)c = 0.026)] [right column; with ω = (0.93 + i0.00065)Ωin]. The solid and dotted curves show the total angular momentum flux and the flux carried by the fluid motion only (the Reynolds stress), respectively. The bottom panels are the blow-up of panel (b) near the magnetic resonances whose locations are indicated by the two vertical lines. . . . . . . . . . . 63
3.9 The net angular momentum flux jump across the two magnetic resonances (upper panel) and the mode growth rate (lower panel) as a function of the dimensionless toroidal magnetic field strength bφ = (vAφ/rΩ)in. The solid lines are for disc sound speed cs = 0.1rΩ, and the dashed line in the lower panel is for cs = 0.05rΩ. The other parameters are the same as in Figs. 7.5 and 7.6: m = 2, σ = 1, q = 0. . . . . . . . . . . . . . . . . . . . . . . 64
3.10 Frequencies of the m = 2 p-modes as a function of the dimensionless toroidal magnetic field strength, (vAφ/rΩ)c, with the upper and lower panels showing the real and imaginary parts, respectively. Different curves correspond to different values of q = d ln Bφ/d ln r (i.e., the slope of Bφ profile) and sound speed cs. The disc density profile is ρ ∝ r−1, i.e. σ = 1. . . . . . . . . . . . 65
xiv

3.11 Frequencies of the m = 2 p-modes as a function of the dimensionless vertical magnetic field strength, (vAz/rΩ)c, with the upper and lower panels showing the real and imaginary parts, respectively. Different curves correspond to different values of the disc toroidal fields. Note that the real mode frequency depends weakly on the magnetic field (for the range of field strength considered), thus bφ = 0.01 corresponds to (vAφ/rΩ)c ≃ 0.0175 and bφ = 0.015 corresponds to (vAφ/rΩ)c ≃ 0.026. The other disc parameters are: σ = 1, q = 0 and cs = 0.1rΩ. . . . . . . . . . . . . . .
3.12 The angular momentum flux carried by wave modes as a function of radius in a BH accretion disc with mixed magnetic fields. The bottom panel is the zoom-in of the upper panel near the resonances. The dotted curves are for a disc with a pure toroidal field, while the solid curves are for a disc with both toroidal and vertical fields. The short-dashed vertical lines indicate the inner and outer magnetic slow resonances (IMR and OMR). The dot-dashed vertical lines represent the corotation resonance (CR). The long-dashed vertical lines (very close to the shortdash lines) indicate the inner/outer Alfve´n resonances (IAR and OAR) which only exist in the mixed field case. Note that the mode frequency in the bz = 0 case is different from the bz = 0.01 case, thus there is a slight shift in location of the MRs and CR between the two cases. Note that bφ = 0.01 and bz = 0.01 correspond to (vAφ/rΩ)c = (vAz/rΩ)c = 0.0176. The other disc parameters are: σ = 1, q = 0 and cs = 0.1rΩ. . . . . . . . . . . . . . . . . .

66 67

4.1 Evolution of the radial velocity amplitude |ur| (evaluated at r = 1.1) for three runs with initial azimuthal mode number m = 2 (left panel), m = 3 (middle panel) and random m (right panel). The blue dashed lines are the fits for the exponential growth stage (between ∼ 10 and ∼ 30 orbits) of the mode amplitude. . . . . . . 72
4.2 Comparison of radial profiles of velocity perturbations from non-linear simulation and linear mode calculation. The top and bottom panels show the normalized radial and azimuthal velocity perturbations, respectively. And the left and right panels are for cases with azimuthal mode number m = 2 and m = 3, respectively. In each panel, the green line is taken from the real part of the complex wavefunction obtained in linear mode calculation, while the blue line is from points with fixed φ = 0.4π at T = 20 orbits, i.e., during the exponential growth stage of nonlinear simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

xv

4.3 Power density spectra of the radial velocity perturbations near the disk inner boundary. Each panel shows the normalized FFT magnitude as a function of frequency. The left and right columns are for runs with initial m = 3 and random m, respectively. In the top and bottom rows, the Fourier transforms are sampled for time periods of [10, 30] orbits and [30, 100] orbits, respectively. . 76
4.4 Power density spectra of the radial velocity perturbations near the disk inner boundary. Each panel shows the normalized FFT magnitude as a function of azimuthal mode number m. The left and right columns are for runs with initial m = 3 and random m, respectively. In top and bottom rows, the time point for FFT are 30, and 100 orbits, respectively. . . . . . . . . . . . . . . . . . . . . 78
4.5 Evolution of radial profiles of velocity perturbations from nonlinear simulations. The top and bottom panels show the radial and azimuthal components of velocity perturbation, respectively. And the left and right panels are for cases with azimuthal mode number m = 2 and m = 3, respectively. In each panel, the data are all taken from points with fixed φ = 0.4π. Different line colors represent different time points during the simulation. . . . 79
4.6 Evolution of radial velocity for runs with initial m = 2 (left), m = 3 (middle) and random m (right), respectively. From top row to bottom row, the time points are 10, 20, 30, 40 and 100 orbits, respectively. Note that the color scale varies from panel to panel. 81
5.1 Schematic description of the cylindrical magnetosphere-disc model considered in this chapter. The thick black line on the left represents the central compact object. Right next to it is a region of low density plasma threaded by both vertical and toroidal magnetic fields. The disc (on the right) consists of fluid of high density (compared to the magnetosphere) and zero or weak magnetic field. These two regions are separated by a thin interface/boundary layer. . . . . . . . . . . . . . . . . . . . . . . . 83
5.2 The wavefunctions for the unstable interface mode with m = 4, ρm/ρd = 1/99, Ωm/Ωd = 1 in a magnetosphere-disc system without (left) and with (right) magnetosphere toroidal magnetic field, respectively. The vertical long-dashed lines mark the position of the interface which separates magnetosphere region and disc region. The solid and dotted lines represent the real and imaginary parts of the wavefunctions, respectively. . . . . . . . . . . . . . . 91
xvi

5.3 Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at the interface) for unstable interface modes as a function of magnetosphere toroidal magnetic field strength. The top and bottom panels are for the cases without and with velocity shear at the interface, respectively, while the left and right panels depict models with the magnetosphere to disc density ratio being 1/99 and 1/9, respectively. Different line types are associated with different values of m, with the solid lines for m = 1, dotted lines for m = 2, short-dashed lines for m = 3, long-dashed lines for m = 4, and dot-dashed lines for m = 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.4 Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable inertial-acoustic disc modes (left panels) and unstable interface modes (right panels) as a function of the magnetosphere toroidal field strength. Different line types represent different azimuthal mode number m. The sound speed profile, density contrast and velocity shear across the interface are the same as in the bottom-left panel of Fig. 7.3. The only difference from Fig. 7.3 is that GR effect (using pseudo-Newtonian potential) is included in the disc rotation profile, which renders both types of modes unstable. On the right panles, the vertical lines (for m = 2 and 3) mark the point when corotation resonance moves into the flow (see Fig. 5.5). . . 100
xvii

5.5 Wave propagation diagram (upper panel) for both interface modes and disc inertial-acoustic modes accompanied by the disc vortensity profile (lower panel). In the upper panel, the three thick solid curves depict the disc rotation profile Ω and Ω ± κ/m, where m is the azimuthal mode number and κ is the radial epicyclic frequency. Note that since κ(rISCO) = 0 these three curves join each other at the Innermost Stable Circular Orbit (which is also the disc inner boundary rin in our setup of Sections 5.4.1 & 5.4.2). The dashed horizontal line marks the point when the corotation radius [where Re(ω)/m = Ω] happens to be exactly at the inner disc boundary. For modes with Re(ω)/m above this line, there is no corotation in the flow; for any modes with Re(ω)/m below this line, the corotation exists in the flow (rc > rin). Three types of wave modes (labeled as I, II and III) are shown in this diagram where the horizontal solid lines indicate the evanescent zones and the wavy curves denote the wave propagation zones. Note that for both Type I and II waves, there is also a tiny wave zone (hardly visible in the figure) in the disc region just outside the ISCO. Type I and II waves are both interface modes while type III waves are disc inertial-acoustic modes. Note that the magnetosphere region (r < rISCO) and the region between the inner and outer Lindblad resonances [where Re(ω)/m = Ω ± κ/m] are wave evanescent zones. . . . . . . . . . . . . . . . . . . . . . . 101
5.6 Wave functions for the three types of wave modes depicted in Fig. 5.5: the interface mode with no corotation resonance in the flow (Type I, left panels), interface mode with corotation resonance in the flow (Type II, middle panles) and disc inertialacoustic mode (Type III, right panels). The azimuthal mode number is m = 3 and all other parameters are the same in Fig. 5.4. The solid and dotted lines show the real and imaginary parts, respectively. In the middle and right panels, the vertical dotdashed lines represent the location of corotation resonance. . . . 102
5.7 Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable inertial-acoustic disc modes (left panels) and unstable interface modes (right panels) as a function of the disc toroidal field strength. The x-axis specifies the ratio of the Alfve´n velocity to disc rotation velocity at the inner boundary rin. The magnetosphere density is set to zero. The disc parameters are rin = rISCO, m = 2, cs = 0.1rΩ and ρ ∝ r−1. For the interface mode (right panels), the two dotted vertical lines from left to right mark the entering of the outer magnetic resonance and inner magnetic resonance into the flow, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
xviii

5.8 The real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable inertial-acoustic modes (left panels) and unstable interface modes (right panels) as a function of the inner disc radius rin (in units of rISCO). The magnetosphere inside rin has zero density and the disc toroidal magnetic field is set to zero. The other disc parameters are m = 2, cs = 0.1rΩ and ρ ∝ r−1. . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.9 Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable interface modes as a function of the disc toroidal field strength, with the inner disc boundary located at rin = 1.5rISCO. The x-axis specifies the ratio of the Alfve´n velocity to disc rotation velocity at rin. The magnetosphere density is set to zero. The disc parameters are rin = rISCO, m = 2, cs = 0.1rΩ and ρ ∝ r−1. The two dotted vertical lines from left to right mark the entering of outer magnetic resonance and inner magnetic resonance into the flow, respectively. . . . . . . . 107
6.1 The two magnetic field profiles adopted in the equilibrium torus model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
6.2 Some examples of Model (a) with a pure toroidal magnetic field (Bφ ∝ r in the fluid and Bz = 0) and constant angular momentum distribution (Ω ∝ r−2). The x-axis is the relative thickness of the torus with ∆r = r2 − r1 being the absolute width, and the y-axis shows the locations of the two boundaries. The different lines represent different values of vAφ0 = vAφ(r0)/(r0Ω0), as indicated. The horizontal line indicates the location of gas pressure maximum (r0 = 1). . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.3 The growth rate of Papaloizou-Pringle instability (in units of Ω0, the fluid rotation rate at the pressure maximum of the torus) as a function of the relative thickness of the torus for different values of m. The rotation profile is Ω ∝ r−2. This figure is similar to Fig. 1 in Blaes & Glatzel (1986) and to Fig. 1 in Abramowicz et al. (1987). 123
6.4 The instability growth rate as a function of the relative thickness of the torus for a range of toroidal magnetic field strengths. Different lines represent different vAφ0 with the solid line denoting the hydrodynamic case. The upper and bottom panels are for the m = 1 and m = 2 modes, respectively. . . . . . . . . . . . . . . 124
6.5 The instability region in the parameter space defined by the dimensionless torus thickness ∆r/r2 and the toroidal magnetic field strength. The other parameters are fixed to m = 2, Ω ∝ r−2 and Bφ ∝ r. The dotted area denotes the region where a growing mode can be found. . . . . . . . . . . . . . . . . . . . . . . . . . . 125
xix

6.6 The instability growth rate as a function of vAφ0 = vAφ(r0)/(r0Ω0) for tori with different thickness ∆r/r2. The upper and bottom panels depict the cases with m = 1 and m = 2, respectively. . . . . 126
6.7 Some special radii for the m = 1 overstable mode in a torus as a function of the dimensionless magnetic field strength vAφ0. The upper and bottom panels correspond to a thin ∆r/r2 = 0.2) and thick (∆r/r2 = 0.7) torus, respectively. The two solid lines show the inner and outer torus boundaries (r1 and r2). The dotted line represents the corotation radius, while the shot-dashed and long-dashed lines show the inner and outer magnetic resonances (IMR and OMR), respectively [see Eq. (6.59)]. The thin solid line is the mode growth rate with the scale shown on the right. The torus rotation and magnetic field profiles are Ω ∝ r−2 and Bφ ∝ r. 127
6.8 The instability growth rate as a function of vAz0 for a torus with relative thickness ∆r/r2 = 0.2. Different lines represent different values of vAφ0. The other parameters are fixed to m = 2, Ω ∝ r−2, Bφ ∝ r and Bz ∝ r in the fluid. . . . . . . . . . . . . . . . . . . . . . 128
6.9 The instability growth rate as a function of vAφ0 for accretion tori described by Model (b) with a pure toroidal magnetic field. The solid and dotted lines correspond to two different inner disc boundary radii. The bottom pane shows the corresponding torus relative thickness. . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
6.10 The instability growth rate as a function of vAz0 for a thin torus (with r1 = 0.9 and ∆r/r2 ≃ 0.2). The different lines represent different values of vAφ0. Unlike Fig. 6.8, here the magnetic field profile is described by Model (b). . . . . . . . . . . . . . . . . . . . 130
7.1 Density profiles for hydrodynamic and hydromagnetic equilibria. The polytropic index is N = 1. The solid, short-dashed and long-dashed lines represent different rotation profiles for nonmagnetic models, while the dotted line is for a magnetic model with WB/|W| = 0.03. . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
7.2 The m = 2 mode frequency as a function of T/|W| with and with out Cowling approximation. The upper and bottom panels show the real and imaginary parts of the frequency, respectively, with Ωc being the rotation frequency at the center. The star has no magnetic field and the polytropic index is N = 1. . . . . . . . . . . 146
7.3 The m = 2 mode frequency as a function of WB/|W| for stellar models with different rotation profiles (thus different T/|W|’s). The upper and bottom panels show the real and imaginary parts of the frequency, respectively, with Ωc being the rotation frequency at the center. The other parameters are the same as in Fig. 7.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
xx

7.4 The m = 2 mode growth rate as a function of T/|W| with and without toroidal magnetic field. The other parameters are A = 0.86 and N = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
7.5 Example wavefunctions of unstable low-T/|W| (≃ 0.095) mode with A = 0.86, C = 0.9, m = 2 and N = 1. The left column shows the radial displacement as a function of radius whereas the right column shows the radial velocity perturbation, with the solid and short-dashed lines representing the real and imaginary parts, respectively. The upper and lower panels are for nonmagnetic and magnetic stellar models, respectively. The dotted lines indicate the location of the corotation resonance (in the nonmagnetic case) or slow magnetosonic resonances (in the magnetic case). The vertical scales of the wavefunctions are arbitrary. . . . 149
7.6 Angular momentum carried by the wave as a function of r. The model parameters are the same as in Fig. 7.5. The upper and lower panels show the nonmagnetic and magnetic models, respectively. The locations of the corotation resonance and slow resonances are indicated by the vertical dotted lines. . . . . . . . 150
xxi

CHAPTER 1 INTRODUCTION
1.1 Quasi-Periodic X-Ray Oscillations and discoseismic modes
Quasi-periodic oscillations (QPOs) in black-hole (BH) X-ray binaries have been intensively studied since the launch of NASA’s Rossi X-ray Timing Explorer (Swank 1999). The observed QPOs can be divided into two classes. The lowfrequency QPOs (about 0.1–50 Hz) are common, observable when the systems are in the hard state and the steep power-law state (also called “very high state”; see Done et al. 2007). They typically have high amplitudes and high coherence (Q ∼> 10), and can vary in frequency on short timescales (minutes). Highfrequency QPOs (HFQPOs) (40-450 Hz) in BH X-ray binaries have also been intensively studied (Fig. 1.1; Remillard & McClintock 2006; Belloni et al. 2011; Altamirano & Belloni 2012), although the signals are much weaker and transient. They are only observed in the “transitional state” (or “steep power law state”) of the X-ray binaries, and have low amplitudes and low coherence. Their frequencies do not vary significantly in response to sizeable (factors of 3-4) luminosity changes. Several systems show pairs of QPOs with frequency ratios close to 2 : 3 (see Remillard & McClintock 2006 for a review). QPOs in X-ray emission have also been detected in the active galaxy RE J1034+396 (Gierlinski et al. 2008; Middleton et al. 2009) as well as in several Ultra-Luminous X-ray sources (ULXs) (e.g., Strohmayer & Mushotzky 2009; Feng, Rao & Kaaret 2010). These are likely the “supermassive” or “massive” analogs of the QPOs detected in Galactic BH X-ray binaries (see Middleton & Done 2010).
Despite much observational progress, the physical origin of HFQPOs re-
1

Figure 1.1: High-frequency quasi-periodic oscillations observed in blackhole binary and black-hole candidate systems. The traces in blue show power density spectrums (PDSs) for the range 1330 keV. Red traces indicate PDSs with a broader energy range, which may be either 2–30 keV or 6–30 keV (Image taken from Remillard & McClintock 2006).
mains elusive. Based on the general consensus that HFQPOs are associated with the dynamics of innermost accretion flows, a number of ideas/models with various degrees of sophistication have been proposed or studied. The general notions of hot blobs exhibiting (test-mass) orbital motion and relativistic precession (Stella et al. 1999), nonlinear resonances (Abramowicz & Kluzniak 2001), and in the case of NS systems, beats between the orbital frequency and stellar spin (Miller et al. 1998), have been the most popular among the observers, although it is not clear how these generic ideas are realized in real fluid dynamical models of accretion flows. Other ideas for NS kHz QPOs are discussed in Bachetti et al. (2010). In the case of HFQPOs in BH X-ray binaries, the only models that go beyond the “test-mass motion” idea are [see Sect. 1 of Lai & Tsang (2009) for a critical review] relativistic diskoseismic oscillation mod-
2

els, in which general relativity (GR) effects produce trapped oscillation modes in the inner region of the disc (e.g., Okazaki et al. 1987; Nowak & Wagoner 1991; see Wagoner 1999 and Kato 2001 for reviews). One can divide various diskoseismic modes according to their vertical structures (see Fig. 1 of Fu & Lai 2009): The so-called p-modes (also called inertial-acoustic modes) have zero node in their wavefunctions along the vertical direction, while the g-modes (also known as inertial modes or inertial-gravity modes) 1 have at least one vertical nodes2. Diskoseismic g-modes have received wide theoretical attentions because they can be trapped by the GR effect without relying on special inner disc boundary conditions (Okazaki et al. 1987), and they may be resonantly excited by global disc deformations (warp and eccentricity) through nonlinear effects (Kato 2003b, 2008; Ferreira & Ogilvie 2008). Kato (2003a) and Li, Goodman & Narayan (2003) (see also Zhang & Lai 2006 and Latter & Balbus 2009) showed that the non-axisymmetric g-mode that contains a corotation resonance (CR, where the wave patten frequency ω/m equals the background rotation rate Ω; here ω is the mode frequency and m is the azimuthal mode number) in the wave zone is heavily damped (see Tsang & Lai 2009a for a similar study for the cmodes). Thus the only non-axisymmetric g-modes of interest are those trapped around the maximum of the (Ω + κ/m) profile (where κ is the radial epicyclic frequency). However, the frequencies of such modes, ω ≃ mΩ(rISCO), 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 BH (Silbergleit & Wagoner 2008; Tassev & Bertschinger 2007).
In Chapter 2, we study analytically the effects of magnetic fields on the rel-
1Note that other than similar mathematical forms of dispersion relations, disk g-modes have no relation to stellar g-modes which are driven by buoyancy.
2The so-called c-modes (corrugation modes) also have at least one vertical node, but they have low frequencies and are not relevant to HFQPOs.
3

ativistic diskoseismic modes in accretion disks around BHs. We consider both poloidal and toroidal fields and use local analysis of the full MHD equations to examine how the magnetic field changes the radial wave propagation diagrams for various modes. We show that the trapping region of g-modes can be easily “destroyed” even when the disk field strength is such that the associated Alfve´n speed is much smaller than the sound speed. On the other hand, the propagation characteristics of p-modes (acoustic oscillations) and c-modes are largely unchanged. Note that since we assume laminar flows for our unperturbed disks, we do not directly address the effects of turbulence on disk modes. However, we believe that our work is relevant to this issue, since magnetic fields naturally arise in a turbulent disk. This elimination of the g-mode cavity is probably the reason for the disappearance of g-modes in recent MHD simulations (O’Neill, Reynolds & Miller 2009). Moreover, several numerical simulations suggest that g-modes may not survive vigorous disc turbulence driven by magneto-rotational instabilities (MRI) (e.g. Arras, Blaes & Turner 2006; Reynolds & Miller 2009).
In Chapter 3, our focus in on p-modes, because the basic wave properties (e.g., propagation diagram) are not affected much by either the disc magnetic fields or the corotation resonance as shown in Fu & Lai (2009) and Tsang & Lai (2009a). In an unmagnetized disc, non-axisymmetric p-modes are trapped between the inner disc edge (at the ISCO) and the inner Lindblad resonance (ILR) radius (where ω˜ = ω−mΩ = −κ). So the existence of global p-modes does require a partially reflective inner disc boundary (see LT09). LT09 showed that such non-axisymmetric p-modes can be overstable for two reasons: (1) Because the spiral density waves inside the ILR carry negative energies while those outside the outer Lindblad resonance (OLR) (where ω˜ = κ) carry positive energies, the
4

leakage of the p-mode energy when the wave tunnels through the corotation barrier (between ILR and OLR) leads to mode growth. (2) More importantly, p-modes can become overstable due to wave absorption at the corotation resonance (CR, where ω˜ = 0). Near the BH, the radial epicyclic frequency κ reaches a maximum before decreasing to zero at the ISCO. This causes a non-monotonic behavior in the fluid vortensity, ζ = κ2/(2ΣΩ) (where Σ is the surface density of the disc), such that dζ/dr > 0 inside the radius where ζ peaks. It can be shown that the sign of the corotational wave absorption depends on the sign of dζ/dr (Tsang & Lai 2008; LT09; see also Goldreich & Tremaine 1979) 3. Thus, p-modes with positive vortensity gradient at the corotation radius can grow in amplitudes due to corotational wave absorption. The resulting overstable non-axisymmetric p-modes were suggested to offer a possible explanation for HFQPOs (LT09; Tsang & Lai 2009b).
Realistic BH accretion disks are expected to contain appreciable magnetic fields as a result of the nonlinear development of the MRI (e.g., Balbus & Hawley 1998). Starting from Hawley et al. (1995), Brandenburg et al. (1995) and Stone et al. (1996) in the 1990s, the property of the MRI turbulence have been studied using numerical simulations under different conditions and setups (local shearing box with/without vertical stratification, with/without explicit viscosity and resistivity, different equations of state, global simulations; e.g., Miller & Stone 2000; Hirose et al. 2009; Guan et al. 2009; Simon et al. 2009; Davis et al. 2010; Fromang 2010; Longaretti & Lesur 2010; Sorathia et al. 2010). Although the precise mechanism of the nonlinear saturation of the MRI remains unclear (see Pessah &
3This applies to barotropic discs, for which the pressure is a function of the density ρ only. For a disc with radial stratification, the radial entropy gradient affects wave absorption at the CR (see Lovelace et al. 1999; Baruteau & Masset 2008), and one can define an effective vortensity, ζeff = ζ/S 2/Γ (where S = P/ΣΓ and Γ is the 2-dimensional adiabatic index), that plays a similar role as ζ (see Tsang & Lai 2009b).
5

Goodman 2009), the simulations generally show that turbulent discs are dominated by large-scale toroidal fields approaching (but less than) equipartition, with somewhat smaller poloidal fields. Therefore, it is important to understand how magnetic fields change the characters of the CR and the overstability of BH disc p-modes. These are the main issues we will address in Chapter 3.
We have recently shown (Fu & Lai 2009), based on local analysis of diskoseismic waves, that the propagation diagram of p-modes, (particularly the wave zone between inner disc edge (at rISCO) and the ILR are not significantly modified by disc magnetic fields. Thus, we expect that the p-mode frequencies in a magnetic disc are not much different from those in a zero-field disc. However, the effect of magnetic field on the corotational wave absorption is more complicated. As we will see, for a pure toroidal field, the CR is split into two magnetic slow resonances, where the Doppler-shifted mode frequency ω˜ matches the frequency of slow MHD waves propagating along the field line (see Sakurai et al. 1991; Goossens et al. 1992; Terquem 2003). With a mixed (poloidal and toroidal) magnetic field, two additional Alfve´n resonances appear, where ω˜ matches the Alfve´n frequency associated with the toroidal field (Keppens et al. 2002). Thus, in general, the original CR is split into four resonances (the inner/outer magnetic slow resonances and the inner/outer Alfve´n resonances), which exhibit a rather rich and complex dynamical behavior.
We note that Tagger and collaborators (Tagger & Pellat 1999; Varniere & Tagger 2002; Tagger & Varniere 2006) developed the theory of accretion-ejection instability for discs containing large-scale poloidal magnetic fields with strengths of order equipartition. Most of their analysis were restricted to vertical fields threading an infinitely thin disc with vacuum outside. They showed that the
6

long-range behavior of the perturbed vacuum field outside the disc greatly enhances the overstability of non-axisymmetric disc oscillation modes. It is currently unclear whether the magnetic field configurations adopted by Tagger et al. can be realized by MRI turbulence – such field configuration may be produced by magnetic field advection from large radii (e.g. Lovelace et al. 2009). It is also unclear if the accretion-ejection instability is robust when more general disc field configurations are assumed. By contrast, in Chaper 3 we focus on the effects of magnetic fields inside the disc. Since the p-modes we are interested in have no vertical structure (i.e., kz = 0), our disc model is essentially a cylinder, with the unperturbed disc properties (e.g., the density ρ, sound speed cs, rotation Ω and magnetic fields Bφ and Bz) depending only on r (the cylindrical radius). While the magnetic field configurations considered in Chapter 3 are prone to MRI, the p-modes are not directly affected by the MRI because kz = 0.
In Chapter 3, we have presented detailed study of the linear instability of non-axisymmetric inertial-acoustic modes (also called p-modes) trapped in the inner-most region of black-hole (BH) accretion discs. This global instability arises because of wave absorption at the corotation resonance (where the wave pattern rotation frequency matches the background disc rotation rate) and requires that the disc vortensity has a positive gradient at the corotation radius. The disc vortensity (vorticity divided by surface density) is given by ζ = κ2/(2ΩΣ), where Ω(r) is the disc rotation frequency, κ(r) is the radial epicyclic frequency, and Σ(r) is the surface density. General relativistic (GR) effect plays an important role in the instability: For a Newtonian disc, with Ω = κ ∝ r−3/2 and relatively flat Σ(r) profile, we have dζ/dr < 0, so the corotational wave absorption leads to mode damping. By contrast, κ is non-monotonic near a BH (e.g., for a Schwarzschild BH, κ reaches a maximum at r = 8GM/c2 and goes to
7

zero at rISCO = 6GM/c2), the vortensity is also non-monotonic. Thus, p-modes with frequencies such that dζ/dr > 0 at the corotation resonance are overstable. Our calculations based on several disk models showed that the lowest-order p-modes with m = 2, 3, 4, · · · have the largest growth rates, with the mode frequencies ω ≃ 0.6mΩISCO (thus giving commensurate frequency ratio 2 : 3 : 4 · · ·), consistent with the High-frequency Quasi-Periodic Oscillations (HFPQOs) observed in BH X-ray binaries (e.g., Remillard & McClintock 2006; Altamirano & Belloni 2012).
So far our works are all based on linear analysis. While these are useful for identifying the key physics and issues, the nonlinear evolution and saturation of the mode growth can only be studied by numerical simulations. It is known that fluid perturbations near the corotation resonane are particularly prone to become nonlinear (e.g., Balmforth & Korycansky 2001; Ogilvie & Lubow 2003). Moreover, real accretion disks are much more complex than any semi-analytic models considered in our previous works. Numerical MHD simulations (including GRMHD) are playing an increasingly important role in unraveling the nature of BH accretion flows (e.g., De Villiers & Hawley 2003; Machida & Matsumoto 2003; Fragile et al. 2007; Noble et al. 2009,2011; Reynolds & Miller 2009; Beckwith et al. 2008,2009; Moscibrodzka et al. 2009,2010; Penna et al. 2010; Kulkarni et al. 2011; Hawley et al. 2011; O’Neill et al. 2011). Despite much progress, global GRMHD simulations still lag far behind observations, and so far they have not revealed clear signs of HFQPOs compatible with the observations of BH x-ray binaries 4. If the corotation instability and its magnetic coun-
4Henisey et al. (2009) found evidence of excitation of wave modes in simulations of tilted BH accretion disks. Hydrodynamic simulations using α-viscosity (Chan 2009; O’Neill et al. 2009) showed wave generation in the inner disk region by viscous instability (Kato 2001). The MHD simulations by O’Neill et al. (2011) revealed possible LFQPOs due to disk dynamo cycles. Note that the observed HFQPOs are much weaker than LFQPOs, therefore much more difficult to obtain by brute-force simulations.
8

terparts studied in our recent papers play a role in HFQPOs, the length-scale involved would be small and a proper treatment of flow boundary conditions is important. It is necessary to carry out “controlled” numerical experiments to capture and evaluate these subtle effects.
In Chapter 4, we use two-dimensional hydrodynamic simulations to investigate the nonlinear evolution of corotational instability of p-modes. Our 2D model has obvious limitations, for example it does not include disc magnetic field and turbulence. However, we emphasize that since the p-modes we are studying are 2D density waves with no vertical structure, their basic radial “shapes” and real frequencies may be qualitatively unaffected by the turbulence (see Arras et al. 2006; Reynolds & Miller 2009; Fu & Lai 2009). Indeed, several local simulations have indicated that density waves can propagate in the presence of MRI turbulence (Gardiner & Stone 2005; Fromang et al. 2007; Heinemann & Papaloizou 2009). Our goal here is to investigate the saturation of overstable p-modes and the their nonlinear pattern frequencies.
1.2 Dynamics of a Magnetosphere-Disc System
In Chapter 5, we study a class of models involving global oscillations of inner accretion discs and magnetosphere boundary layers around NSs or BHs. Our model setup consists of an uniformly rotating magnetosphere with both toroidal and poloidal magnetic fields and low plasma density, surrounded by a geometrically thin (Keplerian) disc with high fluid density and zero (or low) magnetic field. To make semi-analytical linear perturbation analysis feasible, our model is two-dimensional with no vertical dependence (i.e., it is height av-
9

eraged). We consider both Newtonian discs and general relativistic (GR) discs (modelled by the pseudo-Newtonian potential) – we will see that GR can qualitatively affect the linear oscillation modes of the system by changing the disc vorticity profile. This kind of magnetosphere-disc systems can form in accreting NS systems when the NS magnetic field holds off the accretion disc at a certain distance from the stellar surface (e.g., Ghosh & Lamb 1978). They may also be applied to accreting BH systems when magnetic fields advect inwards in the accretion disc and accumulate near the inner edge of the disc (e.g. BisnovatyiKogan & Ruzmaikin 1974, 1976; Igumenshchev, Narayan & Abramowicz 2003; Rothstein & Lovelace 2008). It is usually thought that during the “transitional state” (when HFQPOs are observed) of BH X-ray binaries, the accretion flow consists of a thin, thermal disc, truncated by a hot, tenuous (and perhaps highly magnetized) inner corona (e.g., Done et al. 2007; Oda et al. 2010). While the true nature of this flow is unclear, our simple magnetosphere-disc setup may represent an idealization of such a flow.
Loosely speaking, two types of global oscillation modes exist in our 2D magnetosphere-disc systems: inertial-acoustic modes of the disc and interface modes. Both types of oscillations can be overstable, but driven by different physical mechanisms. Our work represents an extension of related previous studies, which have focused on disc oscillations or interface oscillations separately.
(i) Disc Oscillations: Although relativistic fluid discs can support different types of oscillation modes (Okazaki, Kato & Fukue 1987; Nowak & Wagoner 1991; see Wagoner 1999, Kato 2001 and Kato et al. 2008 for reviews), our focus in Chapter 5, as well in our previous chapters is on the inertial-acoustic modes (also called p-modes). These modes have no vertical structure (i.e., their wave
10

functions, such as the pressure perturbation, have no node in the vertical direction), and are adequately captured by our 2D model.
On the other hand, our magnetosphere-disc model does not capture other discoseismic modes, including g-modes (also called inertial-gravity modes), cmodes and vertical p-modes. All these modes have vertical structures.
In Chapter 3 we have shown that non-axisymmetric p-modes trapped in the inner-most region of a BH accretion disc are a promising candidate for explaining HFQPOs in BH x-ray binaries because their frequencies naturally match the observed values without fine-tuning of the disc properties, and because they become overstable due to wave absorption at corotation point when the general relativistic effect on the disc rotation profile is taken into account. Although a toroidal disc magnetic field tends to suppress the instability (Fu & Lai 2011a), a large-scale poloidal field may enhance the instability (Tagger & Pallet 1999; Tagger & Varniere 2006; Yu & Lai 2012). Most of these previous works adopt a reflective boundary condition at the inner edge of the disc – this is important for the trapping of disc p-modes (see Lai & Tsang 2009). How such a reflection can be achieved was not clear. We will show in Chapter 5 that the discmagnetosphere boundary serves as a “reflector” for the disc p-modes.
(ii) Interface Oscillations: Li & Narayan (2004) were the first to consider the interface modes between a magnetosphere (with a vertical magnetic field) and an incompressible disc in the cylindrical approximation (i.e., no z-dependence). They showed that the interface modes can be subject to Rayleigh-Taylor instability and/or Kelvin-Helmholtz instability, depending on the density contrast and velocity shear across the interface. The mode frequencies roughly scale as mΩin, where m = 1, 2, 3 . . . are the azimuthal mode number and Ωin is the angular
11

frequency of the disc flow at the interface, making the interface modes a viable candidate for explaining HFQPOs when the interface radius rin is suitably adjusted. Tsang & Lai (2009b) generalized the analysis of Li & Narayan (2004) by considering compressible discs. They showed that a relatively large disc sound speed is necessary to overcome the stabilizing effect of disc differential rotation and thereby maintain the mode growth. Besides linear calculations, numerical simulations of the magnetosphere-disc interface have also been presented in a number of papers (see Kulkarni & Romannova 2008; Romanova, Kulkarni & Lovelace 2008 and references therein).
One ingredient that has been missing from both Li & Narayan (2004) and Tsang & Lai (2009b) is toroidal magnetic fields in the magnetosphere, which could be a very important component (Ghosh & Lamb 1978; Ikhsanov & Pustil’nik 1996). In Chapter 5, we generalize previous calculations of global non-axisymmetric modes confined near the interface of the magnetosphere-disc system by taking into account toroidal magnetic fields in the magnetosphere. We aim to examine the possibility of large-scale Rayleigh-Taylor/KelvinHelmholtz instabilities of the interface in the presence of magnetic fields and disc differential rotation.
Overall, the goal of Chapter 5 is to present a complete analysis of large-scale non-axisymmetric modes (including their instabilities) in the magnetosphere boundary and in the surrounding disc, and to assess their viabilities for explaining HFQPOs in accreting compact binaries.
We note that in addition to Li & Narayan (2004) and Tsang & Lai (2009b), our study on the interface modes complements several other works on the instability of magnetized accretion discs. For example, Lubow & Spruit (1995) and
12

Spruit, Stehle & Papaloizou (1995) considered the magnetic interchange instability of a thin rotating disc when a poloidal magnetic field provides some radial support in the disc. Lovelace, Turner & Romanova (2009) studied the instability of a magnetopause for cases where the shear layer has appreciable radial width and found that the Rossby wave instability may arise in the shear layer. Lovelace, Romanova & Newman (2010) considered a setup similar to the one in Chapter 5, but they included only vertical field in the magnetosphere, and focused on small-scale (with radial wavelength ≪ r) Kelvin-Helmholtz modes in the boundary layer.
1.3 Dynamics of magnetized accretion tori
Differentially rotating flows can be found not only in BH & NS accretion discs, but also in other rotating astrophysical flows. Papaloizou & Pringle (1984) discovered that accretion tori can be subjected to a global non-axisymmetric instability that grows on a dynamical time-scale. Accretion tori are bagel-shaped discs with high internal temperatures and well-defined boundaries. They may be representative of certain stages or regions of the inner accretion flows around black holes, such as those found in active galactic nuclei and quasars (e.g., Begelman, Blandford & Rees 1984). They may also form in the gravitational collapse of the rotating core of massive stars (e.g., Woosley 1993) and after the merger of compact neutron star and black hole binaries (e.g., Duez et al. 2009; Etienne et al. 2009; Rezzolla et al. 2010; Montero et al. 2010), and thus are thought to be the central engine of gamma-ray bursts (e.g., Meszaros 2006). The PapaloizouPringle (PP) instability arises from the interaction between non-axisymmetric waves across the corotation radius (rc), where the wave pattern rotation fre-
13

quency Ωp equals the background fluid rotation rate Ω (e.g., Blaes & Glatzel 1986; Goldreich, Goodman & Narayan 1986; Glatzel 1987b). Waves outside the corotation radius (r > rc) have Ωp larger than Ω(r) and carry positive energy, while waves at r < rc have Ωp < Ω(r) and carry negative energy. Instability occurs when the negative-energy waves inside rc lose energy to the positiveenergy waves outside rc, leading the amplification of the wave amplitudes. To maintain the instability, the waves must be efficiently reflected at the inner and outer boundaries so that they are trapped in the torus. The growth rate of the PP instability is maximal for a constant-angular momentum torus. For a very thin torus (with the inner and outer radii close to each other), the instability dis-
√ appears when p = d ln Ω/d ln r > − 3; for a wider torus, the instability persists (with decreasing growth rate) as the Keplerian rotation profile (p = −3/2) is approached (e.g., Papaloizou & Pringle 1985, 1987; Goldreich et al. 1986; Zurek & Benz 1986; Sekiya & Miyama 1988). Other properties of the PP instability, such as its connection with the instability of vortices (e.g., Glatzel 1987a), its nonlinear evolution (e.g., Goodman, Narayan & Goldreich 1987; Hawley 1991) and the effect of accretion (Blaes 1987), have been studied.
Interest in the PP instability waned in the 1990s when Balbus & Hawley (1991) pointed out that the Magneto-Rotational Instability (MRI), originally studied for magnetized Taylor-Couette flows (Velikhov 1959; Chandrasekhar 1961), can be important for astrophysical accretion discs. Since the MRI is robust and requires only a weak magnetic field, the nonlinear development of MRI may lead to efficient angular momentum transport in accretion discs. Over the last two decades, numerous studies have been devoted to the MRI and related issues such as MHD turbulence in the disc (see, e.g., Balbus & Hawley 1998 and Balbus 2003 for reviews; a sample of recent numerical studies include
14

Hirose et al. 2009; Guan et al. 2009; Simon et al. 2009; Davis et al. 2010; Fromang 2010; Longaretti & Lesur 2010; Sorathia et al. 2010).
Nevertheless, the question remains as to what happens to the original PP instability in an accretion torus when a finite magnetic field is present. After all, the tori produced in various astrophysical situations (e.g., binary mergers; Rezzolla et al. 2010; Montero et al. 2010) are expected to be highly magnetized. One might dismiss this question as purely academic since such a magnetic torus is likely MRI unstable and therefore turbulent. We note, however, that the usual MRI operates on perturbations with vertical structure (i.e., with finite vertical wavenumber kz), while the PP instability operates on perturbations with kz = 0. That is, the PP instability pertains to the height-averaged behavior of the disc. Therefore one might expect that the PP instability will continue to operate even in the presence of MRI-induced turbulence. Furthermore, in connection with Galactic black-hole X-ray binaries, it has been suggested that accretion tori can support discrete, trapped oscillation modes, which might explain high-frequency quasi-periodic oscillations (e.g., Strohmayer 2001; Remillard & McClintock 2006) observed in a number of X-ray binary systems (e.g., Rezzolla et al. 2003; Lee et al. 2004; Schnittman & Rezzolla 2006; Blaes et al. 2007; Sramkova et al. 2007; Montero et al. 2007). Although it is currently not clear that pressure-supported tori provide a realistic model for the accretion flow around a black hole in any spectral state, structures resembling pressuresupported tori do appear to be present in some non-radiative global MHD simulations of MRI-driven turbulent accretion flows (e.g., Hawley & Balbus 2002; De Villiers et al. 2003; Machida et al. 2006).
There have been a number of previous studies on global MHD instabili-
15

ties in accretion flows. For example, Knobloch (1992), Kumar et al. (1994), Curry, Pudritz & Sutherland (1994) and Curry & Pudritz (1995) carried out global analysis for the axisymmetric modes with finite kz in differentially rotating flows threaded by vertical and/or azimuthal magnetic fields, thus establishing the robustness of MRI in these flows. Ogilvie & Pringle (1996) studied the non-axisymmetric instability of a cylindrical flow containing an azimuthal field, while Curry & Pudritz (1996) studied a similar flow containing a vertical field. Both studies focused on modes with finite vertical wavenumbers, which inevitably invite MRI. Although the effect of boundaries is emphasized, a somewhat arbitrary rigid boundary condition was adopted in these studies. As far we are aware, the behavior of the PP instability for finite tori with magnetic (both vertical and azimuthal) fields has not been clarified.
In Chapter 6, we carry out global stability analysis of magnetized accretion tori subjected to nonaxisymmetric perturbations as we do in Chapter 3 and Chapter 5. Since our main aim is to understand the effects of magnetic fields on the original PP instability, we focus on modes with no vertical structure (kz = 0) and we pay particular attention to the boundary conditions. As in many previous studies mentioned above, we model the torus by a cylindrical incompressible flow threaded by both vertical and toroidal magnetic fields.
1.4 Dynamics of rotating protoneutron stars
As another example of astrophysical rotating flows, rotating neutron stars (NSs) formed in the core collapse of a massive star or the accretion-induced collapse of a white dwarf maybe subject to nonaxisymmetric instabilities (e.g, Ander-
16

sson 2003; Stergioulas 2003; Ott 2009). The onset and development of these rotational instabilities are often parameterized by the ratio β ≡ T/|W|, where T is the rotational kinetic energy and W the gravitational potential energy of the star. In particular, the dynamical bar-mode (m = 2) instability sets in when β 0.27 and grows on the dynamical time-scale. This critical β, originally derived for incompressible Maclaurin spheroid in Newtonian gravity (Chandrasekhar 1969), is relatively insensitive to the stiffness of the equation of state as long as the degree of differential rotation is not too large (e.g., Toman et al. 1998; New, Centrella & Tohline 2000; Liu & Lindblom 2001), although simulations show that it tends to be reduced by general relativity effect as the compactness of the star M/R increases (Shibata, Baumgarte & Shapiro 2000; Saijo et al. 2001). There also exist secular instabilities, which are driven by some dissipative mechanisms, such as viscosity and gravitational radiation. In the latter case it is known as the Chandrasekhar-Friedman-Schutz instability (Chandrasekhar 1970; Friedman & Schutz 1978). Although the threshold of the secular bar-mode instability (β 0.14) is easier to be satisfied than the dynamical bar-mode instability, it grows on a much longer time-scale due to its dissipative nature (e.g., Lai & Shapiro 1995; Andersson 2003).
The nonlinear development of the dynamical bar-mode instability has been extensively studied in a large number of numerical simulations (Tohline, Durisen & McCollough 1985; Pickett, Durisen & Davis 1996; Cazes & Tohline 2000; Brown 2000; Liu 2002; Shibata & Sekiguchi 2005; Camarda et al. 2009). In the early 2000, it was found that for stars with sufficiently large differential rotation, dynamical instability can develop at significantly lower β than 0.27 (Centrella et al. 2001), even for β on the order of 0.01 (Shibata, Karino & Eriguchi 2002, 2003; Ott et al. 2005; Ou & Tohline 2006; Saijo & Yoshida 2006; Cedra-
17

Duran, Quilis & Font 2007; Corvino et al. 2010). These low-T/|W| appear to have quite different physical origin from the canonical bar-mode instability (see below). Most importantly, recent 3D simulations of a large sample of rotational core-collapse models carried out by Dimmelmeier et al. (2008), which include a state-of-the-art treatment of the microphysics during collapse and the initial rotational profiles obtained from models of precollapse evolution of massive stellar cores, have shown that in many of the models, the proto-NSs exhibit sufficient differential rotation to be subject to the low-T/|W| instability (see also Ott et al. 2007) 5. Such proto-NSs would generate strong gravitational waves (GWs), much stronger than a non-rotating core-collapse would produce, significantly increasing the possibility of detecting GWs from extra-galactic corecollapse supernovae by LIGO and other ground-based GW detectors (Ott 2009). We note that our current understanding of the angular momentum evolution of pre-supernova stars is uncertain, so one cannot predict the rotation profile of the collapsing core with great confidence (Heger, Woosley & Spruit 2005). Therefore the detection (or non-detection) of the rotational signature of proto-NSs by GW detectors (such as Advanced LIGO) may provide valuable information on massive star evolution and the mechanism of core-collapse supernova explosion.
Despite clear numerical evidence for their existence, the physical origin of the low -T/|W| instabilities remains unclear. It has been suggested (Watts, Andersson & Jones 2005; Saijo & Yoshida 2006) that the instabilities are associated with the existence of corotation resonance (where the wave pattern speed equals the background fluid rotation rate) inside the star and are thus likely to be a subclass of shear instabilities which require a certain degree of differential rotation
5By contrast, the threshold for the canonical bar-mode instability is never reached even when the precollapse core has a very large angular momentum, because in that case core bounce would occur at low densities.
18

(Watts, Anderson & Jones 2005; Corvino et al. 2010). Corotation resonance has long been known to be the key ingredient for some instabilities in other astrophysical fluid systems, such as the Papaloizou-Pringle instability for accretion torii (Papaloizou & Pringle 1984; Fu & Lai 2010b) and the corotational instability for thin accretion discs (Narayan, Goldreich & Goodman 1987; Tsang & Lai 2008, 2009; Lai & Tsang 2009; Fu & Lai 2010a). In addition, numerical calculations by Ou & Tohline (2006) suggested that the presence of a local minimum in the radial vortensity profile of the star is also needed to amplify the mode (see also Corvino et al. 2010).
An important issue concerning the low-T/|W| instability is the effects of magnetic fields. Proto-NSs are expected to contain appreciable magnetic fields. In particular, large toroidal fields can be generated from twisting relatively weak poloidal fields by differential rotation or from magneto-rotatioinal instabilities (e.g., Balbus & Hawley 1998; Akiyama et al. 2003; Obergaulinger et al. 2009). While magnetic fields have a negligible effect on the high T/|W| instability (Camarda et al. 2009), it is not clear whether low-T/|W| can survive in the presence of B fields. Indeed, our previous work on magnetized discs showed that even a weak magnetic field can change the structure of corotation resonance significantly (Fu & Lai 2010a). In Chapter 7, as a first step of clarifying this issue, we carry out eigenvalue calculation of the effects of purely toroidal B fields on low T/|W| instability by employing a cylindrical stellar model.
Several chapters of this thesis are based on our published results. Chapter 2 is from Fu & Lai (2009). Chapter 3 is from Fu & Lai (2011a). Chapter 6 is from Fu & Lai (2011b). Chapter 7 is from Fu & Lai (2011c).
19

CHAPTER 2 EFFECTS OF MAGNETIC FIELDS ON THE DISCOSEISMIC MODES OF
ACCRETING BLACK HOLES
2.1 Introduction
The origin of the rapid quasi-periodic variabilities observed in a number of accreting black hole X-ray binaries is not understood. It has been suggested that these variabilities are associated with discoseismic oscillation modes of the black hole accretion disk. In particular, in a disk with no magnetic field, the so-called g-modes (inertial oscillations) can be self-trapped at the inner region of the disk due to general relativistic effects. Real accretion disks, however, are expected to be turbulent and contain appreciable magnetic fields. We show in this chapter that even a weak magnetic field (with the magnetic energy much less than the thermal energy) can modify or “destroy” the self-trapping zone of disk g-modes, rendering their existence questionable in realistic black hole accretion disks. The so-called corrugation modes (c-modes) are also strongly affected when the poloidal field approaches equal-partition. On the other hand, acoustic oscillations (p-modes), which do not have vertical structure, are not affected qualitatively by the magnetic field, and therefore may survive in a turbulent, magnetic disk.
20

2.2 Basic Equations

We consider a non-self-gravitating accretion disk, satisfying the usual ideal

MHD equations:

∂ρ ∂t + ∇ · (ρv) = 0,

∂v ∂t

+

(v

·

∇)v

=

−

1 ρ

∇Π

−

∇Φ

+

1T, ρ

∂B ∂t = ∇ × (v × B).

(2.1) (2.2) (2.3)

Here ρ, P, v are the fluid density, pressure and velocity, Φ is the gravitational

potential, and

B2 1 Π ≡ P + , T ≡ (B · ∇)B
8π 4π

(2.4)

are the total pressure and the magnetic tension, respectively. The magnetic field

B also satisfies the equation ∇ · B = 0. We assume that the fluid obeys the

barotropic equation of state P = P(ρ).

We adopt the cylindrical coordinates (r, φ, z) which are centered on the cen-

tral BH and have the z-axis in the direction perpendicular to the disk plane. The

unperturbed background flow is assumed to be axisymmetric with a velocity field v = rΩ(r)φˆ, and magnetic field B = Bφ(r)φˆ + Bzzˆ, i.e., Bz is constant while Bφ has a radial dependance. Force balance in the unperturbed flow implies

G

≡

1 ∇Π −

1 T

=

Ω2rrˆ

−

∇Φ.

ρρ

(2.5)

Consider perturbations of the form eimφ−iωt. The linearized fluid equations

are

1 ∂ imρ ∂ −iω˜ δρ + r ∂r (ρrδvr) + r δvφ + ∂z(ρδvz) = 0,

δρ 1 ∂

1

−iω˜ δvr − 2Ωδvφ = Gr ρ − ρ ∂r δΠ + ρ(δT)r,

(2.6) (2.7)

21

where

κ2 im 1

−iω˜ δvφ

+

2Ω δvr

=

−

δΠ ρr

+

ρ (δT )φ,

δρ 1 ∂

1

−iω˜ δvz = Gz

ρ

−

ρ

δΠ ∂z

+

ρ (δT)z,

−iω˜ δBr =

imBφ r

+

∂ Bz ∂z

δvr,

−iω˜ δBφ

=

∂ − ∂r (Bφδvr)

+

Bz

∂ ∂z δvφ

−

Bφ

∂ ∂z δvz

+

r

dΩ dr

δBr,

−iω˜ δBz

=

− Bz r

∂ ∂r (rδvr)

−

imBz r

δvφ

+

imBφ r

δvz

,

(2.8) (2.9) (2.10) (2.11) (2.12)

ω˜ = ω − mΩ

(2.13)

is the comoving wave frequency, and

κ2 ≡ 2Ω d(r2Ω) r dr

(2.14)

is the radial epicyclic frequency. δρ, δΠ, δv, δB are Eulerian perturbations, and

ρ, B refer to the unperturbed flow variables. In addition, for barotropic fluid,

we have

11

1

δρ

=

c2s δP

=

c2s (δΠ

−

B 4π

·

δB),

where cs is the sound speed.

(2.15)

To perform local (WKB) analysis, we consider perturbations with spatial de-

pendence eikrr+ikzz. In the leading-order approximation, we keep only the radial

gradient of Ω(r) and Bφ(r) while assuming that the variation scales of all the

other background quantities are much larger than the wavelength of the pertur-

bation, i.e., kr, kz ≫ 1/r. The linearized MHD equations then reduce to

iω˜ − ρc2s

δΠ

+

ikrδvr

+

ikφδvφ

+

ikzδvz

+

iω˜ Bφ 4πρc2s

δBφ

+

iω˜ Bz 4πρc2s

δBz

=

0,

− ikr ρ

δΠ

+

iω˜ δvr

+

2Ωδvφ

+

ikzBz + ikφBφ 4πρ 4πρ

δBr

−

Bφ 2πρr

δBφ

=

0,

(2.16) (2.17)

− ikφ δΠ ρ

−

κ2 2Ω δvr

+

iω˜ δvφ

+

(q

+ 1)Bφ 4πρr

δBr

+

ikzBz + ikφBφ 4πρ 4πρ

δBφ = 0, (2.18)

22

− ikz ρ

δΠ

+

iω˜ δvz

+

ikzBz + ikφBφ 4πρ 4πρ

δBz = 0,

(ikφBφ + ikzBz)δvr + iω˜ δBr = 0,

ikr Bφδvr − ikzBzδvφ + ikzBφδvz − pΩδBr − iω˜ δBφ = 0,

ikr Bzδvr + ikφBzδvφ − ikφBφδvz − iω˜ δBz = 0,

(2.19) (2.20) (2.21) (2.22)

where kφ ≡ m/r. We have assumed Bφ ∼ rq and p ≡ d ln Ω/d ln r. Note that in deriving eqs. (2.16)-(2.22), we have dropped the terms proportional to Gr and Gz in eqs. (2.7) and (2.9): since Gr = (Ω2 − Ω2K)r ∼ Ω2Kr(H/r)2 and Gz ∼ Ω2Kz (where Ωk is the Keplerian frequency, i.e., the angular frequency in the absence of pressure force), Gr is much smaller than the other terms in eq. (2.7) provided that krr ≫ 1 + v2Aφ/c2s, and Gz is also negligible if we focus on the mid-plane of the disk.

2.3 Hydrodynamic Limit: Discoseismic Modes

In the absence of magnetic fields, for kz, kr ≫ kφ, the perturbed MHD equations (16)-(22) lead to the dispersion relation:

(ω˜ 2 − κ2)(ω˜ 2 − kz2c2s ) = kr2c2s ω˜ 2.

(2.23)

For kz = 0 (or ω˜ 2 ≫ kz2c2s), this becomes ω˜ 2 = kr2c2s + κ2, the usual dispersion relation for spiral density wave; for ω˜ 2 ≪ kz2c2s, this becomes ω˜ = ±κkz/(kr2 + kz2)1/2, describing inertial oscillations (e.g., Goodman 1993).

For an accretion disk, with scale height H ≪ r, the vertical dependence of the perturbation is not well described by the plane wave eikzz unless kzH ≫ 1. Okazaki et al. (1987) showed that for a thin disk, the perturbation equations can be approximately separated in r and z (see also Nowak & Wagoner

23

1991, 1992; Ipser 1994). For example, for vertically isothermal disks with constant scale height H, one finds δP(r, z), δvr(r, z), δvφ(r, z) ∝ Hn(z/H), while δvz(r, z) ∝ Hn′ (z/H), where Hn(Z) (with n = 0, 1, · · ·) is the Hermite polynomials and Hn′ (Z) = dHn(Z)/dZ. With this separation of variables, Okazaki et al. (1987) obtained the dispersion relation for a given n:

(ω˜ 2 − κ2)(ω˜ 2 − nΩ2⊥) = kr2c2s ω˜ 2,

(2.24)

where Ω⊥ is the vertical epicyclic frequency and is related to H by H = cs/Ω⊥. An important property of relativistic disks around BHs is that κ is non-monotonic. Three types of trapped modes can be identified (see Fig. 2.1; see also Wagoner 1999, Kato 2001 and Ortega-Rodriguez et al. 2006 for reviews):

(i) P-modes. For n = 0, waves can propagate in the region where ω˜ 2 > κ2. These are acoustic waves (modified by disk rotation), and have also been termed inertial-acoustic modes. If waves can be reflected at the disk inner radius (rISCO, the inner-most stable circular orbit), discrete p-modes can be self-trapped at the inner-most region of the disk (see Fig. 2.1a, 2.1c).

(ii) G-modes. For n ≥ 1, waves can propagate in the region where ω˜ 2 < κ2 <
nΩ2⊥ or ω˜ 2 > Ω2⊥ > κ2 (note that κ < Ω⊥ in GR). The former specifies the g-mode
propagation zone: self-trapped g-modes can be maintained in the region where
κ peaks (for m = 0: see Fig. 2.1b) or in the region where Ω − κ/m < ω/m <
Ω + κ/m (for m 0: see Fig. 2.1d).1 Because these discrete, self-trapped modes
do not require special boundary conditions (e.g., wave reflection at r = rISCO),
they have been the focus of most studies of relativistic diskoseismology. Note
that although we call these g-modes (following the terminology of Kato 2001
1Note that non-axisymmetric g-modes with ω/m < Ω(rISCO) contain corotation resonance in the wave zone, leading to strong damping of the mode (Kato 2003; Li et al. 2003; Zhang & Lai 2006). On the other hand, modes with Ω(rISCO) < ω/m < max(Ω + κ/m) do not suffer corotational damping, and are therefore of great interest.

24

Figure 2.1: Wave propagation diagram showing various trapped modes in BH accretion disks: (a) axisymmetric p-mode; (b) axisymmetric g-mode; (c)non-axisymmetric p-mode; (d) non-axisymmetric g-mode and √c-mode. The curves d√epict various critical frequencies (κ, nΩ⊥, Ω, Ω ± κ/m, Ω ± nΩ⊥/m), the vertical dotted lines denote the inner-most stable circular orbit (ISCO). The curvy horizontal lines specify wave propagation zones and the height of the line is ω (for the top panels) or ω/m (for the bottom panels) of the mode.
and Wagoner 1999), they have no relation to gravity waves, which are driven by buoyancy. Instead, these modes describe inertial oscillations, and have also been termed inertial modes (or inertial-gravity modes).
(iii) C-modes. For n ≥ 1 and m ≥ 1, the wave propagation condition ω˜ 2 >
25

nΩ2⊥

>

κ2

leads

to

an

additional

wave

trapping

region,

where

ω/m

<

Ω−

√ nΩ⊥/m

(see Fig. 2.1d). Note that for spinning BHs, Ω⊥ < Ω. Clearly, these modes exist √
only when Ω − nΩ⊥/m > 0 and wave reflection occurs at r = rISCO. Following

the previous works (e.g., Kato 1990 and Silbergleit et al. 2001, who focused

on the “fundamental” n = m = 1 mode, corresponding to the Lense-Thirring

precession of the inner disk), we call these (“corrugation”) c-modes.

Comparing eqs. (2.23) and (2.24), we see that we can obtain the radial dis-
persion relation of different modes by adopting the vertical “quantization” con√
dition kz = n/H in eq. (2.23), with kz = 0 specifying p-modes. In a generic
disk (e.g., when the disk is not isothermal vertically), the same “quantization”
condition would not hold, but we still expect kz ∼ 1/H ∼ Ω⊥/cs for the (vertically) lowest-order g-mode or c-mode. In the next sections, we will adopt kz = √ηΩ⊥/cs, with η of order unity, when we study how magnetic fields mod-
ify low-order g-modes and c-modes.

Our approach in this paper is based on Newtonian theory. GR effect can be incorporated into our analysis by using the Paczynski-Witta pseudo-Newtonian potential, Φ = −M/(r − 2M). Alternatively, we could simply replace the Newtonian Ω, Ω⊥, κ by their exact general relativistic counterparts (e.g., Okazaki et al. 1987):

M/r3 Ω= ,
1 + a M/r3

4aM1/2 3a2 1/2

Ω⊥ = Ω 1 − r3/2 + r2

,

M(r2 − 6Mr + 8aM1/2r1/2 − 3a2) 1/2 κ = r2(r3/2 + aM1/2)2

(2.25) (2.26) (2.27)

(in geometric units such that G = c = 1), where a = Js/M is the spin parameter of the black hole. In general, Ω ≥ Ω⊥ > κ. In the case of a Schwarzschild BH,

26

Ω = Ω⊥ > κ, with κ peaks at r = 8M and becomes zero at rISCO = 6M. This non-monotonic behavior of the radial epicyclic frequency is preserved for Kerr BHs, and, as discussed above, is the key ingredient for the existence of trapped diskoseismic modes.

2.4 Effect of Poloidal Fields

We first consider the case of a pure poloidal field, with Bφ = 0. For kz, kr ≫ kφ, equations (2.16)-(2.22) then lead the dispersion relation:

ω˜ 6 − [(kz2 + kr2)(c2s + v2Az) + kz2v2Az + κ2]ω˜ 4

+

kz2v2Az

(kz2

+

kr2)(2c2s

+

v2Az)

+

dΩ2 d ln r

+ κ2kz2c2s

−kz4v2Azc2s

(kz2

+

kr2)v2Az

+

dΩ2 d ln r

= 0,

ω˜ 2

(2.28)

where vAz ≡ Bz/ 4πρ. In the incompressible limit, this reduces to the dispersion relation found in, e.g., Balbus & Hawley (1991). For a given k = (kr, kφ, kz), equation (2.28) admits three branches, corresponding to fast, slow magnetosonic waves and Alfve´n wave, all modified by differential rotation. For kz ≫ kr, the Alfve´n branch can become unstable when kz2v2Az < −dΩ2/d ln r. This is the wellknown MRI (e.g., Balbus & Hawley 1998).

2.4.1 P-modes
If kz = 0, equation (2.28) reduces to ω˜ 2 = κ2 + kr2(c2s + v2Az). 27

(2.29)

This is almost the same expression as in pure hydrodynamic case (ω˜ 2 = κ2 +kr2c2s), the only difference being that the sound speed is replaced by fast magnetosonic wave speed, c2s + v2Az. Thus the basic property of p-modes is not affected by poloidal magnetic fields.

2.4.2 G-modes

For a fixed kz = √η/H = √ηΩ⊥/cs (see §3), we can rewrite eq. (2.28) as an expres-

sion for kr2:

(c2s

+

v2Az)kr2

=

(ω˜ 2

− ω21)(ω˜ 2 − ω22)(ω˜ 2 − (ω˜ 2 − ω23)(ω˜ 2 − ω24)

ω25) .

The five critical frequencies are given by

(2.30)

where b ≡ vAz/cs.

ω21 = η(Ω⊥)2,

ω22

=

1 2

κ2 + 2η(Ω⊥)2b2 +

ω23 = η(Ω⊥)2b2,

ω24

=

η(Ω⊥)2 1

b2 + b2

,

ω25

=

1 2

κ2 + 2η(Ω⊥)2b2 −

κ4 + 16η(Ω⊥Ω)2b2 κ4 + 16η(Ω⊥Ω)2b2

(2.31) (2.32) (2.33) (2.34) (2.35)

Equation (2.30) allows us to identify various wave propagation regions (kr2 > 0). We first consider subthermal fields, with b < 1. When b 0.4 (and with η = 1 for the lowest order g-modes), the five critical frequencies satisfy ω21 > ω22 > ω23 > ω24 > 0 > ω25 in the inner region of the disk. Thus there are three wave propagation regions:

Region I : ω˜ 2 > ω21, 28

(2.36)

Figure 2.2: The effect of poloidal magnetic field on the g-mode propagation zone for m = 0, η = 1. In each panel, the upper three curves are ω2 (eq. [2.32]) and the lower two curves are ω3 (eq. [2.33]). Axisymmetric g-modes of frequency ω can propagate in the region where ω3 < ω < ω2. The solid line refers to the case of b = 0, the short-dashed line b = 0.03 and the long-dashed line
b = 0.1, where b ≡ vAz/cs, with vAz = Bz/ 4πρ (Alfve´n speed) and cs the sound speed. The vertical dotted lines correspond to the inner disk radius at ISCO. The curvy horizontal lines specify the wave propagation zones, and the height of the line is ω of the mode. The lower and upper panels are for the case of a Schwarzschild BH (a = 0) and a Kerr BH (a = 0.8M), respectively. The angular frequencies are in units of M−1 = c3/(GM).

Region II : ω23 < ω˜ 2 < ω22, Region III : ω˜ 2 < ω24.

(2.37) (2.38)

Region II corresponds to the original g-mode cavity modified by the mag-

29

netic field; in the zero field limit, ω3 = 0, ω2 = κ and eq. (2.37) reduces to ω˜ 2 < κ2. Fig. 2.2 depicts the critical frequencies ω2 and ω3 for several values of b. This also serves as the propagation diagram for m = 0 g-modes (wave can propagate in region where ω3 < ω < ω2). We see that as the magnetic field increases, the g-mode self-trapping zone gradually shrinks and disappears even when the magnetic field is still very subthermal (for a Schwarzschild BH, this occurs for b 0.08). More precisely, the g-mode cavity can still exist for large b, but it now requires a reflection boundary at rISCO. This behavior can be easily understood by inspecting eq. (32): While κ peaks at some radius rmax, Ω⊥ and Ω both increase monotonically with decreasing r. Since Ω⊥ and Ω are much larger than κ in the inner region of the disk, the 2(Ω⊥b)2 term or the 4Ω⊥Ωb term can dominate over κ2 even when b is still small, therefore making the self-trapping zone disappear. Roughly, this occurs at b bcrit ∼ (κ2/2Ω⊥Ω)rmax.
For non-axisymmetric perturbations (m 0), the wave propagation region II is determined by (i) Ω − ω2/m < ω/m < Ω + ω2/m, and (ii) ω/m > Ω + ω3/m or ω/m < Ω − ω3/m. Fig. 2.3 shows the propagation diagram. As mentioned before (see Footnote 1), for b = 0, only the modes with ω > mΩ(rISCO) are of interest, since otherwise there is a corotation resonance in the wave zone, leading to strong mode damping (Kato 2003; Li, Narayan & Goodman 2003; Zhang & Lai 2006). Thus, self-trapped g-modes reside around the radius where Ω + ω2/m is the maximum (and this maximum arises because κ depends nonmonotonically on r). We see from Fig. 2.3 that this g-mode self-trapping region disappears as b increases. The larger m is, the more fragile is the cavity. For example, the m = 1 cavity disappears for b 0.015, while for m = 2, this occurs for b 0.005.
30

Figure 2.3: The effect of poloidal magnetic field on the g-mode propagation zone for m 0, η = 1, a = 0. The three panels are for b = vAz/cs = 0, 0.005 and 0.015. The solid line, shortdashed lines and the long-dashed lines show Ω, Ω ± ω2 and Ω ± ω2/2, respectively. In the bottom panel, the dotted lines show Ω ± ω3 (In the upper and middle panels, Ω ± ω3 almost concide with Ω, since ω3 = 0 for b = 0 and ω3 ≪ Ω for b ≪ 1). Non-axisymmetric g-modes can propagate in the region where Ω − ω2/m < ω/m < Ω − ω3/m or Ω + ω3/m < ω/m < Ω + ω2/m. The vertical dotted lines correspond to the inner disk radius at ISCO. The curvy horizontal lines in top and middle panels specify wave propagation zones and the height of the line is ω/m of the mode. Note that the self-trapping zone (depicted in the upper and middle panels) disappears as b increases. The angular frequencies are in units of M−1 = c3/(GM).
31

2.4.3 C-modes
For m 0 and η ∼ 1, equation (2.30) also describes trapped c-modes. When b 0.4, ω21 is the largest among all the critical frequencies and the c-mode propagation zone corresponds to Region I (see eq. [2.36]). Note that since ω1 is not affected by the magnetic field, the trapping region is determined by ω/m < Ω − ω1/m = Ω − Ω⊥/m (for η = 1, see the upper panel of Fig. 2.4; cf. Fig. 2.1d). When b 0.4, the ordering between ω1 and ω2 switches and c-modes propagate in the region where ω˜ 2 > ω22, with the trapping zone determined by ω/m < Ω − ω2/m (see the bottom panel of Fig. 2.4). Thus, in the presence of a reflection boundary at rISCO, c-modes are not affected by the poloidal magnetic field when b 0.4, but can be appreciably modified when b 0.4.
From eq. (2.30) we can identify other wave propagation zones (see eq. [2.38]). Fig. 2.5 gives an example, for m = 2, b = 0.7. Note that, except for the c-mode trapping zone discussed above, all the propagation zones are bounded by at least one “singular point” (where kr → ∞). Unlike the turning point (kr → 0) associated with wave reflection, wave absorption is expected to occur at these singular points (see Zhang & Lai 2006; Tsang & Lai 2008a and references therein). Thus, the new wave trapping regions given by eq. (2.38) will not lead to interesting global oscillation modes (Note that in the case of b = 0.7, the ordering of five critical frequencies is different from the one described in §4.2. However, our conclusion still holds true, i.e., there is no chance to form a wave trapping zone bounded by two reflection points other than the c-mode oscillation region, which is bounded by a reflection point and the ISCO).
32

Figure 2.4: The effect of poloidal magnetic field on c-modes with m = 1 and a = 0.2M. The upper panel shows the original c-modes since b is small and ω1 = Ω⊥; in the bottom panel, with a large b, the cmode trapping zone is instead bounded by the inner reflection boudary and Ω − ω2. The vertical dotted line refers to the inner disk radius at ISCO. The curvy horizontal lines specify wave propagation zones and the height of the line is ω of the mode.
2.5 Effect of Toroidal Fields
In this section, we consider the effect of a pure toroidal field, with Bz = 0. Various instabilities may exist for such field geometry, depending upon the rotation profile Ω(r) and the magnetic field profile Bφ(r) (e.g., Acheson & Gibbons 1978; Terquem & Papaloizou 1996). Here we focus on how Bφ affects the diskoseismic modes.
33

Figure 2.5: Wave propagation diagram for non-axisymmetric g-modes and c-modes, for b = vAz/cs = 0.7, m = 2, η = 1, and a = 0. The solid line, dot-short dashed lines, long-dashed lines, dotted lines and short-dashed lines show Ω, Ω ± ω4/2, Ω ± ω3/2, Ω ± ω1/2 and Ω ± ω2/2 (eqs. [2.31]-[2.34]), respectively. The vertical dotted line shows the inner disc radius. The curvy horizontal lines specify wave propagation zones and the height of the line is ω/2 of the mode. The singular points (where kr → ∞) are indicated by filled circles. The angular frequencies are in units of M−1 = c3/(GM).

2.5.1 P-modes

With kz = 0, equations (2.16)-(2.22) reduces to

ω˜ 2 = κ2 + kr2(c2s + v2Aφ),

(2.39)

where vAφ ≡ Bφ/ 4πρ. Thus, the toroidal field affects p-modes in the same way as the poloidal field does (see §4.1).

34

2.5.2 G-modes

Since the general dispersion relation for m 0 is quite complicated, here we focus on axisymmetric perturbations.2 With m = 0, Equations (2.16)-(2.22) lead to

ω4 − [κ2 + (kz2 + kr2)(c2s + v2Aφ)]ω2 + κ2kz2(c2s + v2Aφ) +2(1 − q)v2Aφc2skz2/r2 = 0,

(2.40)

where q = d ln Bφ/d ln r. Solving for kr2, we have

kr2

=

(ω2

− ω2+)(ω2 − (c2s + v2Aφ)ω2

ω2−)

,

(2.41)

with the two critical frequencies given by

ω2±

=

κ2

+

ηΩ2⊥(1 2

+

b2φ)

±

1 2

[κ2 − ηΩ2⊥(1 + b2φ)]2 − 8(1 − q)ηv2AφΩ2⊥/r2,

(2.42)

where bφ ≡ vAφ/cs and we have used kz = √η/H = √ηΩ⊥/cs as in §4. Clearly, for

bφ = 0, eq. (2.42) reduces to eq. (2.24).

When q = 1 (i.e., Bφ ∝ r), eq. (2.42) gives ω2+ = η(Ω⊥)2(1 + b2φ), and ω2− = κ2. Since ω2− is independent of Bφ, the g-mode propagation zone is unaffected no matter how strong the field is. When q 1, as long as vAφ ≪ Ω⊥r (which is valid in most disk situations), the 8(1 − q)ηv2AφΩ2⊥/r2 term in eq. (2.42) represents only a small correction, i.e., ω2− is still very close to κ2. Thus for general toroidal field satisfying vAφ ≪ Ω⊥r, the axisymmetric g-mode propagation zone is not affected by the magnetic field.
2Since c-modes necessarily require m > 0, our analysis here cannot be applied to c-modes.

35

2.6 Summary
In this chapter we have studied the effects of both poloidal and toroidal magnetic fields on the diskoseismic modes in BH accretion disks. Previous works by Kato, Wagoner and others have been based on hydrodynamic disks with no magnetic field. The key finding of our study is that the g-mode self-trapping zone (which arises from GR effect) disappears when the disk contains even a small poloidal magnetic field, corresponding to vAz/cs = 0.01 − 0.1 (see Fig. 2.22.3; vAz is the Alfve´n speed and cs is the sound speed). It is well-known that the combination of a weak poloridal field and differential rotation gives rise to MRI, making real astrophysical disks turbulent. Earlier numerical simulations indicated that the magnetic field grows as MRI develops, until it saturates at vAz/cs ∼ 0.1 − 1, with the toroidal field stronger than the poloidal field by a factor of a few (see, e.g., Hawley et al. 1996; Balbus & Hawley 1998). Recent simulations showed that the turbulent state depends strongly on the net magnetic flux through the disk (e.g., Fromang & Papaloizou 2007; Simon, Hawley & Beckwith 2008). In any case, it is likely that the magnetic field in a turbulent disk is large enough to “destroy” the g-mode self-trapping zone.
Thus, the g-mode properties (including the frequencies and excitations) derived from hydrodynamical models are unlikely to be applicable to real BH accretion disks. The disappearance of the g-mode trapping zone might also explain why Arras et al. (2006) and Reynolds & Miller (2008) did not see any global g-modes in their MHD simulations.
As mentioned in Section 2.1, g-mode oscillations have been considered a promising candidate to explain QPOs in BH X-ray binaries. Theoretically, these
36

modes are appealing because in hydrodynamic disks their existence depends on general relativistic effect and does not require special disk boundary conditions. Our analytical results presented in this chapter, together with recent numerical simulations (Arras et al. 2006; Reynolds & Miller 2008), suggest that magnetic fields and turbulence associated with real accretion disks can change this picture significantly.
While g-modes can be easily modified or “destroyed” by magnetic fields, our analysis showed that p-modes are not affected qualitatively. The magnetic field simply changes the sound speed to the fast magnetosonic wave speed and leaves the p-mode propagation diagram unchanged. We also showed that a weak poloidal field (vAz/cs ≪ 1) does not affect the c-mode propagation zone, although a stronger field modifies it. Our results therefore suggests that global p-mode oscillation is robust and may exist in real BH accretion disks, provided that partial wave reflection at the disk inner edge can be achieved.3 Of particular interest is the non-axisymmetric p-modes, since they may be excited by instabilites associated with corotation resonance (Tsang & Lai 2008, Lai & Tsang 2009).
3Kato (2001) has discussed why such reflection may be possible.
37

CHAPTER 3 COROTATIONAL INSTABILITY, MAGNETIC RESONANCES AND GLOBAL INERTIAL-ACOUSTIC OSCILLATIONS IN MAGNETIZED
BLACK-HOLE ACCRETION DISCS
3.1 Introduction
Low-order, non-axisymmetric p-modes (also referred as inertial-acoustic modes) in hydrodynamic accretion discs around black holes are plausible candidates for high-frequency quasi-periodic oscillations (QPOs) observed in a number of accreting black-hole systems. These modes are trapped in the innermost region of the accretion disc, and are subject to global instabilities due to wave absorption at the corotation resonance (where the wave pattern frequency ω/m equals the disc rotation rate Ω), when the fluid vortensity, ζ = κ2/(2ΩΣ) (where κ and Σ are the radial epicyclic frequency and disc surface density, respectively), has a positive gradient. In this chapter, we investigate the effects of disc magnetic fields on the wave absorption at corotation and the related wave super-reflection of the corotation barrier, and on the overstability of disc p-modes. In general, in the presence of magnetic fields, the p-modes have the character of inertial-fast magnetosonic waves in their propagation zone. For discs with a pure toroidal field, the corotation resonance is split into two magnetic resonances, where the wave frequency in the corotating frame of the fluid, ω˜ = ω − mΩ, matches the slow magnetosonic wave frequency. Significant wave energy/angular momentum absorption occurs at both magnetic resonances, but with opposite signs, such that one of them enhances the super-reflection while the other diminishes it. The combined effect of the two magnetic resonances is
38

to reduce the super-reflection and the growth rate of the overstable p-modes. Our calculations show that even a subthermal toroidal field (with the magnetic pressure less than the gas pressure) may suppress the overstability of hydrodynamic (B = 0) p-modes. For accretion discs with mixed (toroidal and vertical) magnetic fields, two additional Alfve´n resonances appear, where ω˜ matches the local Alfve´n wave frequency. The effect of these additional resonances is to further reduce or diminish the growth rate of p-modes. Our results suggest that in order for the non-axisymmetric p-modes to be a viable candidate for the observed high-frequency QPOs, the disc magnetic field must be appreciably subthermal, or other mode excitation mechanisms are at work.

3.2 Setup and Basic Equations
We consider a non-self-gravitating accretion disc, satisfying the usual ideal MHD equations

∂ρ ∂t + ∇ · (ρv) = 0,

∂v ∂t

+

(v

·

∇)v

=

−

1 ρ

∇Π

−

∇Φ

+

1 (B
4πρ

·

∇)B,

∂B ∂t = ∇ × (v × B),

(3.1) (3.2) (3.3)

where ρ, P, v are the fluid density, pressure and velocity, Φ is the gravitational

potential, and

B2 Π≡P+
8π

(3.4)

is the total pressure. The magnetic field B also satisfies the equation ∇ · B = 0.

We assume the flow obeys the barotropic equation of state P = P(ρ).

39

We adopt the cylindrical coordinates (r, φ, z) which are centered on the central BH and have the z-axis in the direction perpendicular to the disc plane. The unperturbed background disc has a velocity field v = rΩ(r)φˆ, and the magnetic field consists of both toroidal and vertical components B = Bφ(r)φˆ + Bz(r)zˆ. The gravitational acceleration in radial direction is defined as

g

=

dΦ dr

(3.5)

so that −∇Φ = −grˆ and g = rΩ2K > 0, where ΩK is the angular frequency for a test mass (the Keplerian frequency)1 Thus the radial force balance equation reads

ρg = ρrΩ2 − dΠ − B2φ . dr 4πr

(3.6)

To probe the dynamical properties of the magnetized flow, we apply linear

perturbations to the ideal MHD Eqs. (3.1)-(3.3) by assuming the perturbation

of any physical variable f to have the form δ f ∝ eimφ−iωt with m being the az-

imuthal mode number and ω the wave frequency. The background flow and

magnetic field have no z-dependance and we will assume that the perturbation

also has no z-dependance (kz = 0). The resulting linearized perturbation equa-

tions clearly contain variables δv, δρ, δP, δΠ and δB. To simply the algebra, we

define a new variable

δh = δΠ = δP + B · δB . ρ ρ 4πρ

(3.7)

Using ∆v = δv + ξ · ∇v = dξ/dt = −iωξ + (v · ∇)ξ, we find that the Eulerian

perturbation δv is related to the Lagrangian displacement vector ξ by δv = −iω˜ ξ− rΩ′ξrφˆ (prime denotes radial derivative). Also, we have (for barotropic flow)

δρ = δP/c2s

(3.8)

1Here “Keplerian” does not necessarily mean the gravitational potential is Newtonian. See the end of Section 3.2.

40

with cs = dP/dρ being the sound speed. Therefore, we can express all the perturbation quantities in terms of ξr and δh, then further combine the perturbation equations into two equations for these two variables:

dξr dr

=

A11ξr

+

A12δh,

dδh dr = A21ξr + A22δh,

(3.9) (3.10)

where

A11

=

rω˜ 2 (ω2Aφ − Ω2)ω˜ 2 + ω2Aφω2 (c2s + v2A)(ω˜ 2 − m2ω2Aφ)(ω˜ 2 − ω2s)

+ (c2s

+

gω˜ 2 v2A)(ω˜ 2

−

ω2s )

−

ω˜ 2

+ 2mω˜ Ω + m2ω2Aφ r(ω˜ 2 − m2ω2Aφ)

,

(3.11)

ω˜ 4 m2

A12

=

− (c2s

+

v2A)(ω˜ 2

−

m2ω2Aφ)(ω˜ 2

−

ω2s )

+

r2(ω˜ 2

−

, m2ω2Aφ)

(3.12)

A21

=

ω˜ 2

−

m2ω2Aφ

−

4(mω2Aφ + ω˜ Ω)2 ω˜ 2 − m2ω2Aφ

+

r

d dr

(ω2Aφ

−

Ω2)

+

(ω2Aφ

−

Ω2) r ρ

dρ dr

+

g ρ

dρ dr

1

+ (c2s

+

v2A)(ω˜ 2

−

m2ω2Aφ)(ω˜ 2

−

ω2s )

r

(ω2Aφ − Ω2)ω˜ 2 + ω2Aφω2

+ g(ω˜ 2 − m2ω2Aφ) 2 ,

(3.13)

A22

=

−

rω˜ 2 (c2s +

(ω2Aφ − Ω2)ω˜ 2 + ω2Aφω2 v2A)(ω˜ 2 − m2ω2Aφ)(ω˜ 2 − ω2s

)

−

(c2s

+

gω˜ 2 v2A)(ω˜ 2

−

ω2s )

+

2m(mω2Aφ + ω˜ Ω) r(ω˜ 2 − m2ω2Aφ)

−

1 ρ

dρ .
dr

(3.14)

In the above expressions,

ω˜ = ω − mΩ

(3.15)

is the wave frequency in the co-rotating frame,

is the Alfve´n velocity,

|B| vA =
4πρ

ωAφ =

vAφ r

=

r

Bφ 4πρ

(3.16) (3.17)

41

is the toroidal Alfve´n frequency, and

ωs =

cs mωAφ = kφ c2s + v2A

csvAφ , c2s + v2A

is the slow magnetosonic wave frequency for k = (m/r)φˆ.

(3.18)

Equations (7.18)-(7.19) describe the linear perturbations (without vertical wavenumber) in a compressible 2D flow threaded by a mixture of toroidal and vertical magnetic fields. In the hydrodynamic disc limit (B → 0), they reduce to the related equations in Goldreich & Tremaine (1979) with no external forcing potential [see also Tsang & Lai 2008 (TL08 hereafter) and LT09]. They also agree with system (14) in Blokland et al. (2005) 2 in the kz → 0 and barotropic limit, although their equations are cast on slightly different variables (δΠ instead of δh). The Hameiri-Bondeson-Iacono-Bhattacharjee (HBIB) type of equations presented in Blokland et al. (2005) have been derived before by Hameiri (1981) and Bondeson et al. (1987) [see also Mikhailovskii et al. (2009); Goossens et al. (1992) (hereafter GHS92) and reference therein] for a magnetized rotating flow in the absence of external gravity. The equations including gravity were obtained by Keppens et al. (2002) (hereafter KCP02) using the same Frieman-Rotenberg technique (Frieman & Rotenberg 1960). System (5) in KCP02 is by far the most general formulation (with finite Bφ, Bz, m, kz, vφ, vz and adiabatic equation of state) that in principle governs all the MHD waves and instabilities for a compressible MHD fluid with an equilibrium flow in cylindrical geometry. In Chapter 3, we focus on the effects of magnetic fields on the corotational instability of the inertial-acoustic modes (p-modes) in accretion discs, thus Eqs. (7.18)-(7.19) are our main working equations. Recently, Yu & Li (2009) considered a similar system (with finite kz and pure toroidal B field) to study the Rossby wave insta-
2There is a typo in their Eq. (20), where in the last line, the term (B2θω˜ + Fvθ) should be (B2θω˜ + Fvθ)2.
42

bility and they obtained equations with similar structure as Eqs. (7.18)-(7.19). Our equations in the pure toroidal B field limit (Bz → 0) are the same as theirs in the kz → 0 limit. Also, the disc system that we study here are similar to the one considered by Terquem (2003), where instead of solving for global modes the author looked into the interaction between the disc and a planet and the consequences on planet migrations.

Evidently, the coefficients in Eqs. (7.22)-(7.27) are singular when

ω˜ 2 = m2ω2Aφ,

(3.19)

and

ω˜ 2 = ω2s.

(3.20)

In more general situations (with kz 0), mωAφ in Eq. (3.19) generalizes to (k · B)/ 4πρ, which is the Alfve´n frequency, and the slow magnetosonic frequency ωs generalizes to c2s/(c2s + v2A)(k · B)/ 4πρ. Following Sakurai et al. (1991) and KCP02, we shall call (3.19) the Alfve´n Resonance (AR) and (3.20) the Magnetic slow Resonance (MR). Note that these singularities/resonances have also been studied in the specific context of disc-planet interactions. Terquem (2003) considered a 2D disk with a pure toroidal magnetic field and found that a subthermal toroidal B-field is capable of stopping the inward planetary migration as long as the Bφ(r) gradient is negative enough. This analytic finding was later confirmed by 2D MHD simulation results (Fromang et al. 2005). On the other hand, Muto et al. (2008), by employing shearing sheet approximation, performed 3D (i.e., kz 0) calculations for a disk with a pure poloidal magnetic field and showed that type I planetary migration could be outward if the B-filed is strong enough (superthermal).

In the special case of kz = Bz = 0, Eq. (3.19) is no longer a real singularity

43

of the differential equations. In this case, one can show, through some subtle

mathematical manipulations, that the terms ω˜ 2 − m2ω2Aφ in A11, A12, A21 and A22 all get canceled. These four coefficients then reduce to

A11

=

v2Aφ c2s + v2Aφ

ω˜ 2

1 −

ω2s

 

ω˜ 2 r

−

m2c2s r3

−

c2s v2Aφ

ω˜ 2

+

2mω˜ Ω r

−

ω˜ 2 v2Aφ

(rΩ2

−

 g) ,

(3.21)

A12

=

1 − v2Aφ +

c2s

ω˜ 2 − ω˜ 2

m2c2s /r2 − ω2s

,

A21

=

ω˜ 2

−

m2ω2Aφ

−

κ2

−

d

ln(ρr/B2φ) d ln r

c2s

v2Aφ + v2Aφ

(Ω2

−

g/r

−

2c2s /r2)

1 − ω˜ 2 − ω2s

4Ω2ω2s

+

 

mv2Aφ c2s + v2Aφ

(Ω2

2 − g/r + 2c2s/r2)

+4mω˜ Ω

c2s

v2Aφ + v2Aφ

(Ω2

−

g/r

+

2c2s /r2)

,

A22

=

m r

c2s

c2s + v2Aφ

ω˜ 2

1 −

ω2s

 2ω˜ Ω

+

mv2Aφ c2s + v2Aφ

(Ω2

−

g/r

+

 2c2s /r2 )

−d d

ln ρ ln r

+

c2s

1 +

v2Aφ

(rΩ2

−

g

−

2v2Aφ/r),

where the radial epicyclic frequency κ is given by

(3.22)
(3.23) (3.24)

κ = 2Ω d (r2Ω) 1/2 . r dr

(3.25)

This peculiar disappearance of Alfve´n singularity can also be verified from a different perspective. GHS92 derived the jump conditions for ξr and δΠ (called P′ in their paper) across the Alfve´n resonance

[ξr] ∝ gB = mBz/r − kzBφ,

[δΠ] ∝ Bz
(see their Eqs. [42]-[43]). We see that when kz = Bz = 0, both ξr and δΠ are continuous across the Alfve´n resonance, thus the Alfve´n singularity disappears

44

3. In a non-magnetic disc, Eqs. (3.19)-(3.20) reduce to the corotation singularity (TL08).

In Sec. 3.3 and Sec. 3.4, we will simply employ Newtonian gravitational potential, which has free-particle orbital frequency ΩK ∝ r−3/2 and epicyclic frequency κ = ΩK. This way a dimension-free analysis can be easily made without missing the essential physics. However, in Sec. 3.5, in order to make a direct comparison with results in LT09, we will use the Paczynski-Wiita pseudoNewtonian potential (Paczynski & Wiita 1980)

GM Φ=− ,
r − rs

(3.26)

where M is the central BH mass and rs = 2GM/c2 is the Schwarzschild radius. The corresponding free-particle orbital frequency and radial epicyclic frequency

are

1 dΦ 1/2

GM 1

ΩK = r dr

=

, r r − rs

(3.27)

κ = ΩK

r − 3rs . r − rs

(3.28)

√

Note that κ2 peaks at r = (2 + 3)rs and decreases to zero at r = 3rs. Through-

out Chapter 3, we will consider thin discs with the sound speed cs ≪ rΩK and

magnetic field satisfying B2/(4πρ) ≪ (rΩK)2, so that Ω ≃ ΩK.

3Although the analysis in GHS92 does not consider external gravity, one can show that the
effect of gravity can be easily included through replacing the term ρv2φ/r in their Eq. (5) by ρv2φ/r+ ρg. This modification only changes the relations between the background variables (Bφ, vφ, cs, etc.) The jump conditions for ξr and P′ would remain the same.

45

3.3 Effective potential and Wave Propagation diagram

We can combine Eqs. (7.18)-(7.19) into a single second-order differential equa-

tion on δh

d2 d dr2 δh + C1(r) dr δh + C0(r)δh = 0,

where

C1(r)

=

−A11

−

A22

−

1 A21

dA21 , dr

C0(r)

=

A11 A22

−

A12 A21

−

dA22 dr

+

A22 A21

dA21 . dr

We introduce a new variable

(3.29)
(3.30) (3.31)

r1

η = exp

2C1(r)dr δh,

(3.32)

so that Eq. (3.29) can be rewritten as a wave equation for η:

d2η dr2 − V(r)η = 0,

(3.33)

where

V (r)

=

1 4

C1

(r)2

+

1 2

dC1(r) dr

−

C0(r)

(3.34)

is the effective potential for wave propagation. Considering a local plane wave

solution

r
η ∝ exp i kr(s)ds ,

(3.35)

Eq. (3.33) then yields

kr2 = −V(r).

(3.36)

Formally, wave can propagate only in the region with V(r) < 0. Our derivation above is quite general and can be applied to complicated systems as the details have all been captured in the four coefficients Aij. Figures 3.1-3.3 show the effective potential for various situations discussed below.

46

3.3.1 Hydrodynamical discs

For a non-magnetic disc, the four coefficients Aij simplify to

d ln(rρ) 2mΩ

A11 = − dr

−

, rω˜

1 m2 A12 = − c2s + r2ω˜ 2 ,

A21 = −D,

2mΩ

A22 =

, rω˜

where D = κ2 − ω˜ 2. Thus C0(r) and C1(r) now read

D m2 2mΩ d Ωρ

C0(r) = − c2s − r2 − rω˜

ln , dr D

dD

C1(r)

=

− dr

ln

rρ

.

The effective potential then takes the following form

(3.37) (3.38) (3.39) (3.40)
(3.41) (3.42)

D m2 2mΩ V(r) = c2s + r2 + rω˜

d ln Ωρ dr D

+

S

1/2

d2 dr2

S

−1/2,

(3.43)

where S = D/(rρ). The effective potential V(r) is singular at both the corota-

tion resonance (CR), where ω˜ = 0, and the Lindblad resonances (LRs), where

D = 0 (see Figs. 3.1-3.3). However, the singularities at the LRs are spurious

singularities of the wave equation (i.e. they can be removed by rewriting the

equation in alternative variables [e.g., Narayan et al. 1987; TL08]) and only the

CR is a non-removable singularity. The behavior of V(r) around the CR depends

on dζ/dr (evaluated at the corotation radius rc), the gradient of the background

fluid vortensity, ζ, defined by

κ2 ζ= .
2Ωρ

(3.44)

For dζ/dr < 0 a narrow “Rossby wave zone” lies just inside rc, while for

dζ/dr > 0 the Rossby zone lies just outside rc; in the special case of dζ/dr = 0,

the corotation singularity disappears (see TL08).

47

Figure 3.1: Wave propagation diagram in accretion discs when the background fluid vortensity gradient (evaluated at the corotation resonance radius) (dζ/dr)c = 0. In the upper panel, the solid line shows the effective potential V(r) as a function of r for B = 0 discs, and the dashed line shows V(r) for discs with finite toroidal magnetic fields. Waves can propagate only in the region where V(r) < 0. The labels ILR, OLR and CR denote the inner Lindblad resonance, outer Lindblad resonance and corotation resonance, respectively. The lower panel shows the blow-up of the region near the CR (located at rc = 1), and the long-dashed line shows the case with a stronger Bφ than the short-dashed line. The labels IMR and OMR denote the inner magnetic (slow) resonance and outer magnetic resonance, respectively. Note that the vertical scales in the upper and bottom panels differ by a large factor. A Newtonian potential is used in calculating V(r). (A color version of this figure is available in the online journal.)
48

Figure 3.2: Same as Fig. 3.1, except for the case (dζ/dr)c < 0.
3.3.2 Discs with pure toroidal magnetic fields
As mentioned in Sec. 3.2, for a disc with a pure toroidal magnetic field, the CR splits into two magnetic slow resonances (MRs). This can be seen clearly from the wave propagation diagram shown in Figs. 3.1-3.3. For B = 0, the effective potentials reduce to those shown in TL08. In the upper panels of Figs. 3.1-3.3, we see that for small magnetic fields (with ωAφ ≪ Ω), the effective potential profile is not affected much by the B-field except near the corotation. The blow-ups in the bottom panels show that the toroidal magnetic field splits the CR into two MRs, the inner magnetic resonance (IMR) and outer magnetic resonance (OMR), one on each side of the CR. A larger toroidal field (see the long-dashed lines) gives a
49

Figure 3.3: Same as Fig. 3.1, except for the case (dζ/dr)c > 0.

wider radial separation between the MRs and a shallower potential well around rc. Note that independent of the sign of (dζ/dr)c, there always exists an effective wave zone centered around corotation point and bounded by the two MRs.

Far away from the resonances (the regions on the left and right sides of the upper panels of Figs. 3.1-3.3), the WKB dispersion relation (3.36) reduces to (see Terquem 2003; Fu & Lai 2009)

ω˜ 2 = κ2 + kr2(c2s + v2Aφ).

(3.45)

This describes fast magnetosonic waves modified by the disc rotation. Thus the spiral density waves (inertial-acoustic waves) inside the ILR and outside the OLR are inertial-fast magnetosonic waves.

50

Figure 3.4: Numerical values of f2+ and f2− [see Eq. (3.47)] as a function of the dimensionless ratio vAφ/(rΩ) (evaluated at r = rc), which specifies the toroidal magnetic field strength in the disc. Dimensionless units are adopted, where G = M = rc = 1, m = 2, and cs/(rΩ) ≪ 1. Different line types represent different disc density and magnetic field profiles.

The behavior of the effective potential around the MRs is more difficult to

analyze. The IMR and OMR are separated from the CR by the distance

rmr − rc ≃ ±

2ωs 3mΩ

rc,

(3.46)

(the upper sign is for the OMR and the lower sign for the IMR). Because of the complexity of the equations involved, no simple analytical expression for V(r) can be derived, even in the asymptotic limit (e.g., ω˜ → ±ωs and vAφ/rΩ ≪ 1). Nevertheless, a careful examination of Eqs. (3.21)-(3.24), (3.30)-(3.31), (3.34) and numerical calculation of V(r) show that in the neighborhood of the MRs, the

51

effective potential can be written in the following form:

V (r)

=

f0

+

f1 ω˜ 2 − ω2s

+

f2+ (ω˜ − ωs)2

+

(ω˜

f2− + ωs)2

,

(3.47)

where f0, f1, f2± vary slowly around the CR/MRs region. Thus, the effective

potential is dominated by second-order singularities at the MRs, i.e., V(r) ∝

(r − rmr)−2 (see Figs. 3.1-3.3). The coefficients f2± can be calculated numerically,

and some examples of how f2+ and f2− depend the magnetic field strength are

shown in Fig. 3.4. We see that in dimensionless units (with G = M = rc = 1,

m = 2 and cs ≪ rΩ, vAφ ≪ rΩ), both f2+ and f2− are negative and approximately equal to −2. Thus, near one of the MRs (e.g., the OMR), we have kr2 = −V(r) ∼ 2/(ω˜ − ωs)2 ∼ 0.2/(r − rmr)2, or kr ∼ 0.4/|r − rmr|. Since kr2 is smaller than |dkr/dr|, the WKB analysis is not really valid near the MRs. On the other hand, close to

the CR radius (r = rc), we have V(r) ∼ −4/ω2s, which implies kr ∼ 2ω−s 1. Since the separation between the IMR and OMR is 2ωs/3Ω (for m = 2), we find that the

WKB phase variation across the CR region (but avoiding the MRs) is of order

unity. Therefore, despite the deep effective potential at the MRs, physically there

is not much a wave zone around the MRs/CR region. This feature is borne out

in our numerical calculations of global wave modes presented in Sec. 3.5.

It is worth noting with the second-order singularities associated with the MRs, the B → 0 limit is approached in a non-trivial way. From Fig. 3.4, we see that as the magnetic field decreases, f2+ and f2− reach the same values (≃ −2.25). Thus the singularities associated with f2+/(ω˜ − ωs)2 and f2−/(ω˜ + ωs)2 simply become 1/ω˜ 2 as B → 0. On the other hand, f1 can be written as f1 = f10ω˜ 0 + f11ω˜ + f12ω˜ 2 + · · ·, with f1i being near constant around the MRs. As B → 0, we have f10 + f2+ + f2− = 0, and the second-order singularity term (∝ 1/ω˜ 2) in V(r) disappears. What is left is the term f11/ω˜ , a first-order singularity associated with the CR in a non-magnetic disc. Evidently, as B → 0, f11 is proportional to

52

the vortensity gradient, dζ/dr.

3.3.3 Discs with pure vertical magnetic fields

Before examining the case of mixed magnetic fields, it is useful to consider discs

with pure vertical fields. Here we assume that Bz is constant (independent of r)

to simplify the analysis. The coefficients Aij of the wave equations (7.18)-(7.19)

are

A11

=

− c2s

c2s + v2A

d ln ρ dr

−

2mΩ rω˜

−

1 ,
r

1 m2

A12

=

− c2s

+

v2A

+

, r2ω˜ 2

A21

=

ω˜ 2

−

κ2

−

c2s v2A c2s + v2A

d ln ρ 2 ,
dr

A22

=

− c2s

v2A + v2A

d ln ρ dr

+

2mΩ . rω˜

We then have

d

C1

=

− dr

ln

A21 rρ

,

+

m2 r2ω˜ 2

c2s v2A c2s + v2A

D m2 2mΩ d Ωρλ

C0

=

− c2s

+

v2A

−

r2

−

rω˜

ln dr A21

d ln ρ dr

2

−

c2s

v2A +

v2A

d ln ρ dr

d ln(A21/r)

d +

dr dr

v2A d ln ρ c2s + v2A dr

,

where

(3.48) (3.49) (3.50) (3.51) (3.52)
(3.53)

λ = (c2s − v2A)/(c2s + v2A)

(3.54)

and D = κ2 − ω˜ 2. The resulting effective potential reads

D m2 2mΩ V(r) = c2s + v2A + r2 + rω˜

d Ωρλ ln
dr A21

+

S

1/2

d2 dr2

S

−1/2

−

m2 r2ω˜ 2

c2s v2A c2s + v2A

d ln ρ dr

2

+

c2s

v2A + v2A

d ln ρ dr

d ln(A21/r) dr

d −
dr

v2A d ln ρ c2s + v2A dr

,

(3.55)

53

Figure 3.5: Effective potential around the corotation/magnetic resonance region for a disc with mixed magnetic fields (solid line). The dotted line shows V(r) for the case with a pure toroidal field. The corotation resonance (CR), inner/outer magnetic (slow) resonances (IMR and OMR) and inner/outer Alfve´n resonances (IAR and OAR) are indicated. Note that the MRs represent second-order singularities in V(r), while the ARs (which exist only for the mixed field case) are first-order singularities.
where S = A21/ (rρ). Compared with the V(r) in the B = 0 case, we see that a constant vertical field does not split the CR. However, in contrast to the firstorder corotation singularity (1/ω˜ ) in the B = 0 case, the CR in discs with Bz 0 manifests mainly as a second-order singularity (∼ 1/ω˜ 2).
54

3.3.4 Discs with mixed magnetic fields

For a disc with mixed magnetic fields, the expression for the effective potential V(r) is more complicated than the pure toroidal field case. In Fig. 3.5, we show a sketch of V(r) when both Bz and Bφ are non-zero. We see that in addition to the two MRs (IMR and OMR) (already present in the pure toroidal field case), two Alfve´n resonances (ARs) appear when ω˜ = ±mωAφ. Since mωAφ > ωs, these ARs lie farther away from the corotation radius than the MRs, at the distance

rar − rc = ±

2ωAφ 3Ω

rc.

(3.56)

Unlike the MRs, the ARs are first-order singularities of the wave equation, i.e., V(r) ∝ 1/(ω˜ − mωAφ) ∝ 1/(r − rar) near the ARs.

3.4 Super-Reflection of the Corotation Barrier

As shown in Sec. 3.3, away from the CR/MRs/ARs region, the WKB dispersion

relation is simply

ω˜ 2 = κ2 + kr2(c2s + v2A).

(3.57)

Thus waves can propagate either inside the ILR (where ω˜ = −κ) or outside the OLR (where ω˜ = κ), and the region between the ILR and OLR represents the “corotation barrier”, where waves are evanescent except for various possible singularities associated with the CR/MRs/ARs. We will be interested in the wave modes trapped between the disc inner edge and ILR. As in the hydrodynamic (B = 0) case, the growth rates of such p-modes depend directly on the wave reflectivity of the corotation barrier (TL08, LT09).

55

Consider launching a density wave δh ∝ exp(−i r krdr) from a radius inside the ILR (r < rILR) (assuming kr > 0; the minus sign arises from the fact that group velocity and phase velocity have different signs where r < rILR, see TL08). The wave carries negative energy since its pattern speed ω/m is smaller than background rotation rate Ω. When the wave impinges on the corotation barrier, it is partly bounced back as a reflection wave δhR ∝ R exp(i r krdr), which carries negative energy, and partly transmitted through the evanescent zone and travels beyond the OLR as a transmission wave δhT ∝ T exp(i r krdr), which carries positive energy. R and T are reflection and transmission coefficients, respectively (see TL08). Let the net absorption of wave energy at the CR, MRs and ARs be Dabs. Energy conservation then implies

|R|2 = 1 + |T |2 + Dabs.

(3.58)

Obviously, wave transmission through the corotation barrier (the |T |2 term) always tends to increase |R|2 beyond unity. However, the resonant wave absorption term Dabs can be either positive or negative.

3.4.1 Non-magnetic Discs

TL08 showed that, for a B = 0 barotropic disc, the wave super-reflection (|R|2 −

1) is dominated by the corotational wave absorption, which depends on the

dimensionless parameter

ν=

cs d ln ζ

23 = β σ− ,

pκ dr c 3

2

(3.59)

where ζ is the vortensity [see Eq. (3.44)], β, p and σ are defined by β = cs/(rΩ),

Ω ∝ r−p and ρ ∝ r−σ (around the CR), respectively, and the subscript “c” im-

plies that the quantity is evaluated at r = rc. The second equality in Eq. (3.59)

56

assumes (Newtonian) Keplerian discs (p = 3/2). Importantly, the corotational

wave absorption is proportional to ν or dζ/dr. This can be understood from the

fact that for ν > 0 (ν < 0), the Rossby zone lies outside (inside) rc (see Figs. 3.1-

3.3), implying a positive (negative) wave energy absorption. Thus, if neglecting

|T |2, we have

dζ

|R| − 1 ∝ ν ∝

.

dr c

(3.60)

Tsang & Lai (2009b) further generalizes the analysis to non-barotropic flows

with radial stratification, in which case an effective vortensity (which depends

on the radial entropy profile of the disc) plays a similar role as ζ.

3.4.2 Numerical calculations of reflectivity

For discs with magnetic fields, we are not able to derive analytical expressions for the reflectivity. Thus we calculate R numerically.

For definiteness, we consider a (Newtonian) Keplerian disc with

Ω = κ ∝ r−3/2,

ρ ∝ r−σ,

Bφ ∝ Bz ∝ rq,

cs = β, rΩ

(3.61)

where σ, q and β are constants. The strength of the field is characterized by

the ratio bφ = vAφ/(rΩ) and bz = vAz/(rΩ), evaluated at rc. We choose the wave

frequency ω = ωr + iωi, with 0 < ωi ≪ ωr. The transmitted wave outside the

OLR takes the form

r
δh ∝ A exp i kr dr ,

(3.62)

where kr > 0 is given by Eq. (3.36) or approximately by Eq. (3.45), and A(r) can be obtained from Eq. (3.32) with η ∝ kr−1/2. Thus the boundary condition at

rout > rOLR is

dδh dr |rout = δh

ikr + A′/A

|rout ,

(3.63)

57

Figure 3.6: Reflection coefficient as a function of ν for a Newtonian Keplerian disk. Different lines represent different magnetic field strengths, measured by the dimensionless parameters bφ = (vAφ/rΩ)c and bz = (vAz/rΩ)c. The other parameters are m = 2, Bφ ∝ Bz ∝ rq with q = 0 and β = cs/(rΩ) = 0.1. The parameter ν is defined according to Eq. (3.59).

where A′ = dA/dr. Starting from r = rout, we integrate Eqs. (7.18)-(7.19) inwards

to a radius inside the ILR. At r = rin < rILR, the wavefunction has the form

rr
δh ∝ A exp −i kr dr + R exp i kr dr .

(3.64)

The reflectivity R can then be extracted from the equation

(−ik + A′/A)δh − δh′

|R| =

(ik + A′/A)δh − δh′

.
rin

(3.65)

In practice, at rin (rout) sufficiently inside rILR (outside rOLR), the A′/A term can be

dropped since kr ≫ |A′/A|.

Figure 3.6 shows the numerical results of the reflection coefficient as a function of ν, defined by Eq. (3.59). For B = 0, our result recovers that obtained

58

by TL08. We see even a relatively weak (sub-thermal, with vAφ ∼< cs) toroidal field can significantly affect the reflectivity. As the magnetic field increases, the reflectivity decreases regardless of the sign of ν. A finite vertical field has a similar effect, although it is sub-dominant compared to the toroidal field effect. We will show in the following section that this result is consistent with the stability property of global disc p-modes.

3.5 Global P-Modes of Black-Hole Accretion Discs

In this section, we present our calculations of the global non-axisymmetric pmodes in BH accretion discs. Recall that for a nonmagnetic disc, these are the inertial-acoustic modes partially trapped between the inner disc edge (at rin = rISCO) and the ILR. For a disc with finite magnetic fields, these become inertialfast magnetosonic waves.

We employ the standard shooting method (Press et al. 1992) to solve

Eqs. (7.18)-(7.19) for the complex eigenfrequency ω = ωr + iωi. The boundary

condition at rout > rOLR is the outgoing-radiative condition discussed in Sec. 3.4.2

[see Eq. (3.63)]. At the inner boundary rin = rISCO, we apply a free surface bound-

ary condition,

∆h = ∆Π/ρ = 0,

(3.66)

i.e., the Lagrangian perturbation of the total pressure vanishes. This boundary condition corresponds to a idealized situation in which no wave energy is lost at the inner disc radius. As discussed in LT09 (and references therein), the inner edge of a real BH disc is more complicated and may involve wave energy loss due to radial infall of the accreting gas. Since our focus here is to understand

59

how disc magnetic fields affect the corotational instability, we will adopt the simple boundary condition (3.66) for all the numerical calculations presented in this section.

As mentioned in Sec. 3.2, we use the Paczynski-Wiita pseudo-Newtonian potential [Eq. (3.26)] to mimic the GR effect. The disk density profile and magnetic field profile take the power-law forms

ρ ∝ r−σ, Bφ ∝ Bz ∝ rq.

(3.67)

The thermal effect and the magnetic field effect are characterized by the dimen-

sionless parameters

β = cs , rΩ

bφ

=

vAφ rΩ

|rin

,

bz

=

vAz rΩ

|rin

,

(3.68)

with β being constant.

The angular momentum flux carried by the wave across the magnetized disc

is given by (e.g. Pessah et al. 2006)

F(r)

=

πr2ρRe(δvrδv∗φ)

−

1 4

r2Re(δBrδB∗φ

),

(3.69)

where ∗ denotes complex conjugate. The first and second terms in (3.69) are

related to the Reynolds stress and Maxwell stress, respectively. The solution of

Eqs. (7.18)-(7.19) gives ξr(r) and δh(r). The expressions for δvr, δvφ, δBr and δBφ can be derived from the original linearized perturbation equations (3.1)-(3.3):

δvr = −iω˜ ξr,

δvφ

=

 −

κ2 2Ω

−

mω2Aφ ω˜

+

mv2Aφ rω˜

 A11

ξr

+





m rω˜

+

mv2Aφ rω˜

A12

 

δh,

δBr

=

imBφ r

ξr

,

δBφ =

−q Bφ r

−

A11 Bφ

ξr − A12Bφδh.

(3.70) (3.71) (3.72) (3.73)

60

3.5.1 Results for discs with pure toroidal magnetic fields

Figure 7.5 gives an example of the eigenfunctions for overstable p-modes in magnetized BH accretion discs. For comparison, the corresponding results for a B = 0 disc are also shown. Clearly, with a small B field (bφ ≪ 1), the wavefunctions away from the corotation/magnetic resonances (CR/MRs) are only slightly modified. The wave action is still trapped in the inner disc region. The most noticeable effects occur near the MRs, where both δh and δvr experience significant variations. In this example, the B = 0 disc has dζ/dr > 0 at the CR, so that corotational wave absorption leads to a mode growth rate ωi ≃ 0.0027Ωin. The inclusion of a small toroidal B field, bφ = 0.01 [corresponding to (vAφ/rΩ)c = 0.0175], however, reduces the mode growth rate to ωi ≃ 0.0018Ωin. Note that the disc sound speed in the example is cs = 0.1rΩ, so the bφ = 0.01 corresponds to (vAφ/cs)c = 0.175.

The origin of the diminished mode growth rate due to magnetic fields can be understood by examining the angular momentum flux carried by the wave modes. Figure 7.6 gives some examples. We see that in the B = 0 case, there is a jump of the angular momentum flux across the CR,

∆F = F(rc−) − F(rc+), (B = 0).

(3.74)

The fact that ∆F > 0 implies a positive wave energy (angular momentum) absorbed at the CR [i.e., Dabs > 0 in Eq. (3.58)]; this leads to super-reflection and is the main driver for the overstability of hydrodynamic p-modes. For discs with a finite Bφ, the CR is split into two MRs (the IMR and OMR), and we see from Fig. 7.6 that significant flux jumps occur at both the IMR and OMR. The flux jump at the IMR, ∆FIMR = F(rIMR−) − F(rIMR+), is negative, implying a negative angular momentum absorption. The flux jump at the OMR,

61

Figure 3.7: Example wavefunctions of m = 2 p-modes in BH accretion discs. The disc density profile is ρ ∝ r−1 (i.e., σ = 1) and the sound speed is cs = 0.1rΩ (i.e., β = 0.1). The left columns are for a B = 0 disc, with the eigenvalue ω = (0.9324 + i0.0027)Ωin (where Ωin is the disc rotation rate at r = rin = rISCO); the right panels are for a disc with bφ = (vAφ/rΩ)in = 0.01 [corresponding to (vAφ/rΩ)c = 0.0175] (with Bφ independent of r, i.e., q = 0), with ω = (0.9312 + i0.0018)Ωin. In each subfigure, the upper panels show different variables as a function of r in the full fluid zone. The solid and dotted lines depict the real and imaginary parts, respectively. The bottom panels zoom in the region near the corotation resonance, with the vertical lines indicating the locations of the CR (left) and magnetic slow resonances (right). The vertical scales of the wavefunctions are arbitrary, with δh(rout) = 1.
62

Figure 3.8: Angular momentum flux carried by the wave as a function of r. The disc parameters are σ = 1, q = 0 and β = 0.1, and the wave modes have m = 2. Panel (a) shows the result of a B = 0 disc [with the mode frequency ω = (0.9324 + i0.0027)Ωin] with the vertical line indicating the corotation resonance. Panel (b) is for a magnetized disc with bφ = 0.01 [corresponding to (vAφ/rΩ)c = 0.0175] [left column; with ω = (0.9312 + i0.0018)Ωin] or bφ = 0.015 [corresponding to (vAφ/rΩ)c = 0.026)] [right column; with ω = (0.93+i0.00065)Ωin]. The solid and dotted curves show the total angular momentum flux and the flux carried by the fluid motion only (the Reynolds stress), respectively. The bottom panels are the blow-up of panel (b) near the magnetic resonances whose locations are indicated by the two vertical lines.

∆FOMR = F(rOMR−) − F(rOMR+), however, is positive, implying a positive angu-

lar momentum absorption. The net angular momentum jump across the MRs

region is then

∆F = F(rIMR−) − F(rOMR+) = ∆FIMR + ∆FOMR.

(3.75)

Note that although both |∆FIMR| and |∆FOMR| are much larger (by orders of magnitude) than ∆F in the B = 0 case, the net jump ∆F is smaller for discs with finite Bφ. As Bφ increases, |∆FIMR| and |∆FOMR| both decrease, and ∆F becomes smaller and may even change sign (See Fig. 3.9). This decrease of wave absorption at

63

Figure 3.9: The net angular momentum flux jump across the two magnetic resonances (upper panel) and the mode growth rate (lower panel) as a function of the dimensionless toroidal magnetic field strength bφ = (vAφ/rΩ)in. The solid lines are for disc sound speed cs = 0.1rΩ, and the dashed line in the lower panel is for cs = 0.05rΩ. The other parameters are the same as in Figs. 7.5 and 7.6: m = 2, σ = 1, q = 0.
the MRs due to finite magnetic fields results in a reduced or even diminished super-reflection (see Sec. 3.4), leading to a reduction of the growth rate and even stabilization of disc p-modes.
Figure 3.9 shows the p-mode growth rate as a function of bφ and the corresponding net angular momentum jump across the CR/MRs region. Obviously, the decrease of ∆F with increasing bφ directly correlates with the decrease of the mode growth rate. Note that when ∆F is zero or slightly negative, the p-mode still has a finite growth rate. This is because wave transmission across the corotation barrier always tends to increase the reflectivity [see Eq. (3.58)], thereby
64

Figure 3.10: Frequencies of the m = 2 p-modes as a function of the dimensionless toroidal magnetic field strength, (vAφ/rΩ)c, with the upper and lower panels showing the real and imaginary parts, respectively. Different curves correspond to different values of q = d ln Bφ/d ln r (i.e., the slope of Bφ profile) and sound speed cs. The disc density profile is ρ ∝ r−1, i.e. σ = 1.
promoting the the mode growth. With a further increase of bφ, the mode growth becomes completely suppressed. For example, for discs with cs = 0.1rΩ, the angular momentum absorption at the CR/MRs change sign at bφ = 0.0132, while the mode growth rate becomes zero at bφ ≃ 0.017 , or (vAφ/cs)in ≃ 0.17 [corresponding to (vAφ/cs)c = 0.3]. For discs with cs = 0.05rΩ, the mode growth rate vanishes at bφ ≃ 0.01, or (vAφ/cs)in ≃ 0.2 [corresponding to (vAφ/cs)c = 0.34].
Figure 3.10 shows the real and imaginary parts of the m = 2 mode frequency as a function of the dimensionless field strength vAφ/(rΩ) (evaluated at rc) for different disc sound speeds and different values of q (measuring the slope of Bφ). Clearly, for such a small field strength, the real mode frequency is approximately
65

Figure 3.11: Frequencies of the m = 2 p-modes as a function of the dimensionless vertical magnetic field strength, (vAz/rΩ)c, with the upper and lower panels showing the real and imaginary parts, respectively. Different curves correspond to different values of the disc toroidal fields. Note that the real mode frequency depends weakly on the magnetic field (for the range of field strength considered), thus bφ = 0.01 corresponds to (vAφ/rΩ)c ≃ 0.0175 and bφ = 0.015 corresponds to (vAφ/rΩ)c ≃ 0.026. The other disc parameters are: σ = 1, q = 0 and cs = 0.1rΩ.
independent of Bφ and q, since the propagation zone of the wave mode is hardly modified by the magnetic field. The mode growth rate, however, is significantly affected because of the modification to the CR. Not surprisingly, Im(ω) depends on Bφ mainly through the ratio (vAφ/rΩ)c.
3.5.2 Results for discs with mixed magnetic fields
The effects of vertical magnetic fields and mixed fields on the disc p-mode growth rate are illustrated in Fig. 3.11 and Fig. 3.12. We see that the real mode
66

Figure 3.12: The angular momentum flux carried by wave modes as a function of radius in a BH accretion disc with mixed magnetic fields. The bottom panel is the zoom-in of the upper panel near the resonances. The dotted curves are for a disc with a pure toroidal field, while the solid curves are for a disc with both toroidal and vertical fields. The short-dashed vertical lines indicate the inner and outer magnetic slow resonances (IMR and OMR). The dot-dashed vertical lines represent the corotation resonance (CR). The long-dashed vertical lines (very close to the short-dash lines) indicate the inner/outer Alfve´n resonances (IAR and OAR) which only exist in the mixed field case. Note that the mode frequency in the bz = 0 case is different from the bz = 0.01 case, thus there is a slight shift in location of the MRs and CR between the two cases. Note that bφ = 0.01 and bz = 0.01 correspond to (vAφ/rΩ)c = (vAz/rΩ)c = 0.0176. The other disc parameters are: σ = 1, q = 0 and cs = 0.1rΩ.
67

frequency depends very weakly on the field strength unless vA becomes comparable to or larger than cs. The solid line in Fig. 3.11b shows that for the pure vertical field case, the p-mode growth rate decreases with increasing Bz. This is consistent with our finding in Sec. ?? that a finite Bz tends to reduce the reflectivity (see Fig. 3.6). When the toroidal field is also included, we see the similar trend (see the dotted and dashed lines). Thus both vertical and toroidal magnetic fields tend to suppress the corotational instability. The effect of the toroidal magnetic field is somewhat larger that the vertical field, given that the required toroidal field strength Bφ to completely suppress the mode growth is smaller than the required Bz. This is understandable since the p-modes have no vertical structure (with kz = 0), and fluid motion in the disc does not directly bend the vertical field lines.
In Fig. 3.12, we compare the angular momentum flux profiles for the pure toroidal field case and the mixed field case. From the upper panel of this figure we see that for the mixed field case, the net flux jump across the resonance zone is smaller than the pure toroidal field case. This explains the somewhat smaller mode growth rate (ωi/Ωin = 0.0017) for the mixed field case as compared to the pure toroidal field case (ωi/Ωin = 0.0018). Note that the real part of the eigenfrequency is also modified slightly by the vertical field. This gives rise to the shift of “plateau” in the bottom panel of Fig. 3.12.
3.6 Summary
The main finding of this chapter is that overstable non-axisymmetric p-modes tend to be stabilized by disc magnetic fields. This implies that, in order for p-
68

modes to explain HFQPOs observed in black-hole X-ray binaries, the disc magnetic field must be sufficiently weak (with the magnetic pressure appreciably less than the thermal pressure) so that the corotational instability is not completely suppressed. Alternatively, some other mode excitation mechanisms are needed. One possibility could be mode excitations by turbulent viscosity (see Kato 2001 for a general discussion of viscous driving of disc oscillation modes) or by vorticity perturbations associated with the turbulence (see Heinemann & Papaloizou 2009). Another possibility could be the accretion-ejection instability studied by Tagger et al. (see Tagger & Pellat 1999; Varniere & Tagger 2002; Tagger & Varniere 2006), which involves large-scale equipartition vertical magnetic fields threading a thin disc embedded in a vacuum or a tenuous corona. Whether these other possible excitation mechanisms can compete with the stabilizing effects of the internal disc magnetic fields remains to be studied. Nevertheless, our finding is limited by the fact that the disc model we considered is essentially an infinite cylinder without vertical structure, whereas in reality disc background (density, magnetic field, etc.) varies in vertical direction. More realistic disc models certainly need to be investigated in future studies.
Finally, we note that although in this chapter we have studied the role of magnetic fields in the context of black-hole diskoseismic oscillations and QPOs, the dynamical effects of magnetic fields on the corotation resonance examined here are quite general. Our finding in this chapter suggests that dynamical instabilities in other rotating astrophysical flows where the corotation resonance plays an important role may be significantly affected by magnetic fields. Examples include the Papaloizou-Pringle instability in accretion tori and the global rotational instability in differentially rotating stars. We will study some of these issues in separate chapters (Chapters 6 and 7).
69

CHAPTER 4 COROTATIONAL INSTABILITY IN BLACK-HOLE ACCRETION DISCS:
NONLINEAR SIMULATIONS
4.1 Introduction
In Chapter 3, we have presented detailed linear theory of a global, non-axisymmetric hydrodynamic/hydromagnetic instability in thin, twodimensional disks. In the hydrodynamic case, the disk becomes unstable when the vortensity gradient at the corotational radius is positive. In this chapter, we use two-dimensional hydrodynamic simulations to investigate the nonlinear evolution of corotational instability.

4.2 Numerical Setup

The disk is assumed to be inviscid and geometrically thin so that the hydrdynamical equations can be reduced to two-dimension with vertically integrated quantities. We adopt an isothermal EOS P = c2sΣ throughout this study, where P is the vertically integrated pressure, Σ is the surface density and cs is the constant sound speed. Self-gravity and magnetic fields are omitted.

We use the Paczynski-Witta Pseudo-Newtonian potential (Paczynski & Witta

1980) to mimic the GR effect:

GM Φ=− ,
r − rS

(4.1)

where rS = 2GM/c2 is the Schwarzschild radius. The corresponding Keplerian

70

rotation frequency and epicyclic frequency are

ΩK = κ=Ω

GM 1 ,
r r − rS r − 3rS . r − rS

(4.2) (4.3)

In our compuation, we will adopt the units such that the inner radius (at the

Inner-most Stable Circular Orbit or ISCO) of the disc is at r = 1.0 and the Kepe-

rian frequency at the ISCO is ΩK = 1. In these dimensionless units, rs = 1/3,

and

21 1

ΩK

=

3 r − rs

√. r

(4.4)

Our computation domain extends from r = 1.0 to r = 4.0 in th radial direction

and from φ = 0 to φ = 2π in the azimuthal direction. We also use the Keplerian

orbital period (T = 2π/ΩK = 2π) at r = 1 as the unit for time. The equilibrium

state of the disk is axisymmetric. The surface density profile has a simple power-

law form

Σ0 = r−1,

(4.5)

and the equilibrium rotation frequency of the disk is given by the force balance

in radial direction

Ω0 =

r(r

4/9 − 1/3)2

−

c2s r2

.

(4.6)

Throughout our simulation, we will consider cs = 0.1 so that Ω0 ≃ ΩK.

We solve the Euler equations with the PLUTO code 1 (Mignone et al. 2007), a Godunov-type code with multiphysics and multialgorithm modules. For this study, we choose a Runge-Kutta scheme (for temporal integration) and piecewise linear reconstruction (for space integration) to achieve second order accuracy, and Roe solver as the solution of Riemann problems. The gird resolution
1publicly available at http://plutocode.ph.unito.it/

71

log|ur |max(r =1.1) log|ur |max(r =1.1) log|ur |max(r =1.1)

0 m=2 −1 −2 −3 −4 −5 −6 straight line with slope = 0.17
0 20 40 60 80 100
Time (orbits)

0 m=3 straight line with slope = 0.2
−1
−2
−3
−4
−5
−6 0 20 40 60 80 100
Time (orbits)

0 random m straight line with slope = 0.2
−1
−2
−3
−4
−5
−6 0 20 40 60 80 100
Time (orbits)

Figure 4.1: Evolution of the radial velocity amplitude |ur| (evaluated at r = 1.1) for three runs with initial azimuthal mode number m = 2 (left panel), m = 3 (middle panel) and random m (right panel). The blue dashed lines are the fits for the exponential growth stage (between ∼ 10 and ∼ 30 orbits) of the mode amplitude.

we adopt is (Nr × Nφ) = (1024 × 2048) so that each grid cell is almost square. Each one of our runs lasts 100 orbits (Keplerian orbits at inner disk boundary) and takes about 10 hours on 64 CPUs. To compare with the linear calculations (Lai & Tsang 2009), the inner disk boundary is reflective with zero radial velecity. At the outer disc boundary, we adopt the outgoing wave radiative boundary condition. This is implemented using the wave damping method (de Val-Borro et al. 2006) to reduce wave reflection at the outer boundary.

4.3 Results
In general, we carry out simulations with two different types of initial conditions for the surface density perturbation. In the first types of runs, we choose δΣ(r, φ) = Random(r, φ) to be purely random noise; in in the second types of runs, we impose δΣ(r, φ) = Random(r) cos(mφ) so that it has an initial azimuthal mode number m. In all cases, the initial surface density perturbation has a small amplitude ( |δΣ/Σ0| ≤ 0.0001 at r = 1).
72

Figure 4.2: Comparison of radial profiles of velocity perturbations from non-linear simulation and linear mode calculation. The top and bottom panels show the normalized radial and azimuthal velocity perturbations, respectively. And the left and right panels are for cases with azimuthal mode number m = 2 and m = 3, respectively. In each panel, the green line is taken from the real part of the complex wavefunction obtained in linear mode calculation, while the blue line is from points with fixed φ = 0.4π at T = 20 orbits, i.e., during the exponential growth stage of non-linear simulation.
73

Table 4.1: Comparison of results from linear and nonlinear studies of unstable disk p-mode

ma ωrb

ωic ωr/mΩd ωr1e |ωr − ωr1|/ωrf ωr2g |ωr2 − ωr1|/ωr2h

2 1.4066 0.0632 0.7033 1.3998

0.5%

1.4296

2.1%

3 2.1942 0.0733 0.7314 2.1997

0.3%

2.2445

2.0%

4 3.0051 0.0763 0.7512 2.9996

0.2%

3.0594

2.0%

5 3.8294 0.0751 0.7659 3.7995

0.8%

3.8886

2.3%

6 4.6621 0.0714 0.7770 4.6494

0.3%

4.7749

2.6%

7 5.5007 0.0664 0.7858 5.4992 0.03% 5.6756

3.1%

8 6.3436 0.0607 0.7930 6.3492 0.09% 6.3189

0.5%

a Azimuthal mode number b Mode frequency from linear calculation (in units of Keplerian orbital fre-
quency at inner disk boundary; same for ωi and ν) c Mode growth rate from linear calculation d Ratio of wave pattern speed to the Keplerian orbital frequency at inner
disk boundary e Mode frequency during the exponential growth stage of nonlinear simu-
lation (peak frequency of the power density between ∼ 10 orbits and ∼ 30 orbits) f Difference between ωr (linear result) and ωr1 (nonlinear result) g Mode frequency during the saturation stage of nonlinear simulation (peak frequency of the power density between ∼ 30 orbits and ∼ 100 orbits) h Difference between ωr1 and ωr2

74

Figure 4.1 shows the evolution of the radial velocity amplitude at near the inner disc radius for runs with initial azimuthal mode number m = 2, m = 3 and random m. This velocity amplitude is obtained by searching for the maximum ur at a given r by varying φ. We see that in all cases there are three different stages for the amplitude evolution. The first stage occupies roughly the first 10 orbits, during which the initial perturbation starts to affect the flow and presumably excites many modes/oscillations in the disc. In the second stage which runs from ∼ 10 orbits to ∼ 30 orbits, the fastest growing mode becomes dominant and undergoes an exponential growth with its amplitude increasing by about 5 orders of magnitude. In the last stage (beyond ∼ 30 orbits), the perturbation growth saturates and its amplitude remains at approximately the same level. A fit to the exponential growth stage gives the growth rate of the fastest growing perturbation. For the m = 2 run, we find 0.17/ log10(e)/2π ≃ 0.0637 (in units of the orbital frequency at rin) as the growth rate, which is quite consistent with the result from our linear eigenmode calculation (Lai & Tsang 2009; Fu & Lai 2011a), ωi ≃ 0.0632 (the imaginary part of the eigenfrequency). For the m = 3 run, our simulation gives 0.074 as the mode growth rate, clsoe to the value 0.0733 from the linear calculation. In Fig. 4.5, we plot the radial profile of velocity perturbation at a fixed φ coordinate from the exponential growth stage of the simulation. And we compare it with the wavefunctions obtained from linear mode calculation. Note that in each panel of Fig. 4.5 we have normalized different sets of data to make sure they have the same scale. We can see that wavefunctions obtained from two studies also agree quite well. This agreement plus the growth rate in Fig. 4.1 affirms that our non-linear simulations capture accurately the the same unstable modes predicted from the linear perturbation analysis.
To explore the time variability of the perturbation, we carry out Fourier
75

Figure 4.3: Power density spectra of the radial velocity perturbations near the disk inner boundary. Each panel shows the normalized FFT magnitude as a function of frequency. The left and right columns are for runs with initial m = 3 and random m, respectively. In the top and bottom rows, the Fourier transforms are sampled for time periods of [10, 30] orbits and [30, 100] orbits, respectively.
transform of the radial velocity ur(r, φ, t) at fixed r and φ for different time segments. In Fig. 4.3 we show some example of the resulting power density spectra (normalized to the maximum value of unity). Different rows correspond to different total sampling times and different columns represent runs with different initial surface density perturbation. In the left column, the disk undergoes an initial perturbation with a fixed azimuthal mode number m = 3. A mixture of different modes/oscillations are excited in the flow. After another 10 orbits of evolution, one of them (the fastest growing mode) has its oscillation ampli-
76

tude grew by a large amount such that it dominates over other modes. This corresponds to the primary spike in the top-left panel. The rest spikes in the same panel are simply harmonics of this primary spike. Table 4.1 shows that the frequency2 of this fastest growing mode (ωr1) is different from the frequency obtained in linear mode calculation by only 0.3%, which again demonstrates the consistency of these two studies. After the perturbation saturates (bottom-left panel), we see that the basic structure of the power density spectrum does not vary much except spikes are not as clean as before probably due to the interaction of different modes. The location of the primary spike (ωr2 in Table 4.1) increased slightly by 2.0%. In Table 4.1 we also include the results from both linear mode calculation and numerical simulation for modes with other mode number m. The comparison illustrates two points: Firstly, the frequencies of the fastest growing modes during the exponentially growth stage of numerical simulations are different from the linear calculation results by less than 1%; Secondly, the frequencies of the fastest growing modes during the saturation stage are only slightly higher (except for m = 8 mode) than ones during the exponential growing stage. These imply that the mode frequencies obtained in linear mode calculation are quite robust and can be reliably applied to explain the observed HFQPOs.
In the right column of Fig. 4.3, the simulation starts with a random initial perturbation (i.e., mixture of modes with different m). During the exponential growth stage, there are four prominent modes (four spikes) in the flow, whose mode number are m = 3, 4, 5, and 6, respectively. They happen to be the modes with highest growth rates based on linear mode calculations (Table 4.1). Their frequencies do not change much in the saturation regime, and more modes
2All the frequencies in this chapter are angular frequencies unless otherwise noted.
77

Figure 4.4: Power density spectra of the radial velocity perturbations near the disk inner boundary. Each panel shows the normalized FFT magnitude as a function of azimuthal mode number m. The left and right columns are for runs with initial m = 3 and random m, respectively. In top and bottom rows, the time point for FFT are 30, and 100 orbits, respectively.

show up in the disk.

To explore the azimuthal variability of the perturbation, we do Fourier trans-

form of radial velocity perturbation near the disk inner boundary at different

time points of the simulation. Fig. 4.4 shows the normalized FT magnitude

|F(m)|2 as a function of azimuthal mode number m where F(m) is computed as

follows

2π

F(m) =

e−imφur(1.1, φ, t)dφ.

0

(4.7)

In the left column, we see the fastest growing mode dominates in the expo-

78

Radial vecocity perturbation

Azimuthal vecocity perturbation

m=2

0.04

T=20 orbits (enlarged by 62 times) T=30 orbits

0.02

0.00

−0.02

−0.04 1.0 1.5 2.0 2.5 3.0 3.5 4.0
Radius [rISCO]
m=2 T=20 orbits (enlarged by 45 times) T=30 orbits
0.04
0.03
0.02
0.01
0.00
−0.011.0 1.5 2.0 2.5 3.0 3.5 4.0
Radius [rISCO]

Azimuthal vecocity perturbation

Radial vecocity perturbation

m=3 0.06 T=20 orbits (enlarged by 14 times)
T=30 orbits
0.04
0.02
0.00
−0.02
−0.04
1.0 1.5 2.0 2.5 3.0 3.5 4.0
Radius [rISCO]
m=3 0.04 T=20 orbits (enlarged by 26 times)
T=30 orbits
0.02 0.00 −0.02 −0.04 −0.06 −0.08 −0.10 −0.121.0 1.5 2.0 2.5 3.0 3.5 4.0
Radius [rISCO]

Figure 4.5: Evolution of radial profiles of velocity perturbations from nonlinear simulations. The top and bottom panels show the radial and azimuthal components of velocity perturbation, respectively. And the left and right panels are for cases with azimuthal mode number m = 2 and m = 3, respectively. In each panel, the data are all taken from points with fixed φ = 0.4π. Different line colors represent different time points during the simulation.

nential growth stage. Modes with frequencies that are simply harmonics of the primary frequency also exist in the disk. After saturation, interactions between these modes bring out other modes with all kinds of m. In the right column, the simulation starts with random initial perturbation. This means a mixture of modes with different m start to grow with approximately the same initial amplitude. Yet they have different growth rates. At the end of the exponential growth
79

stage, the modes with higher growth rates stand out. Fig. 4.5 shows the comparison of wavefunctions during the mode growing
stage (T = 20 orbits) and at the end of the growing stage (T = 30 orbits) for simulations with different initial perturbations. At T = 20 orbits, the oscillation mainly comes from the single fastest growing mode and takes a smooth radial profile. At T = 30 orbits, perturbation starts to saturate, oscillation now consists of many different modes, and its radial profile exhibits sharp variations at several places which may be the sites of shock formations.
To see this evolution from a different point of view, in Fig. 4.6 we show the color contour of redial velocity for runs with different initial perturbations (different columns) at different time points (different rows). Very visibly we see the spiral feature gradually appears and develops, as the system evolves it becomes less and less smooth, and more sharper variations emerges (compare the second row with the fourth row). We suspect this might indicate the formation of shock waves and may have something to do with the perturbation saturation. At the end (the last row) the spiral arms becomes fairly irregular which we think are related to the emergence and interaction of many modes with multiple values of m (see fig. 4.4).
80

Figure 4.6: Evolution of radial velocity for runs with initial m = 2 (left), m = 3 (middle) and random m (right), respectively. From top row to bottom row, the time points are 10, 20, 30, 40 and 100 orbits, respectively. Note that the color scale varies from panel to panel.
81

CHAPTER 5 DYNAMICS OF THE INNERMOST ACCRETION FLOWS AROUND COMPACT OBJECTS: MAGNETOSPHERE-DISC INTERFACE, GLOBAL
OSCILLATIONS AND INSTABILITIES
5.1 Introduction
In this chapter, we study global non-axisymmetric oscillation modes and instabilities in magnetosphere-disc systems, as expected in the innermost accretion flows in neutron star X-ray binaries and possibly also in accreting black bole systems. Our two-dimensional magnetosphere-disc model consists of a Keplerian disc in contact with an uniformly rotating magnetosphere with low plasma density. Two types of global modes exist in such systems, the interface modes and the disc inertial-acoustic modes, both of which can become unstable. We examine various physical effects and parameters that influence the properties of these oscillation modes, particularly their growth rates, including the magnetosphere field configuration, the velocity and density contrasts across the magnetosphere-disc interface, the rotation profile (Newtonian vs General Relativistic), sound speed and magnetic field of the disc. The interface modes are driven unstable by Rayleigh-Taylor and Kelvin-Helmholtz instabilities, but can be stabilized by the toroidal field (through magnetic tension) and disc differential rotation (through finite vorticity). General relativity increases the their growth rates by modifying the disc vorticity outside the magnetosphere boundary. The interface modes may also be affected by wave absorption associated with corotation resonance in the disc. In the presence of a magnetosphere, the inertial-acoustic modes of a relativistic disc are effectively trapped at the inner-
82

Figure 5.1: Schematic description of the cylindrical magnetosphere-disc model considered in this chapter. The thick black line on the left represents the central compact object. Right next to it is a region of low density plasma threaded by both vertical and toroidal magnetic fields. The disc (on the right) consists of fluid of high density (compared to the magnetosphere) and zero or weak magnetic field. These two regions are separated by a thin interface/boundary layer.
most region of the disc just outside the interface. They are driven unstable by wave absorption at the corotation resonance, but can be stabilized by modest disc magnetic fields. The overstable oscillation modes studied in this chapter have characteristic properties that make them possible candidates for the highfrequency quasi-periodic oscillations observed in X-ray binaries.
5.2 Equilibrium and Perturbation equations
We consider a magnetosphere-disc model similar to the one in Tsang & Lai (2009b). It consists of a magnetosphere region where magnetic pressure dominates over gas pressure and a disc region which has high density compared to the magnetosphere (see Fig. 5.1). These two regions are separated by an inter-
83

face (boundary layer). Unlike Tsang & Lai (2009b), who considered only vertical

magnetic field for the magnetosphere, we take into account both vertical and

toroidal fields. In the disc region, gas pressure dominates over magnetic pres-

sure. Since any initial poloidal field is likely to generate a dominating toroidal

field due to the disc differential rotation, for simplicity we take the disc to be

threaded by toroidal B field only. We assume that flows in both regions are

non-self-gravitating, satisfying the usual ideal MHD equations

∂ρ ∂t + ∇ · (ρv) = 0,

(5.1)

∂v 1 1 ∂t + (v · ∇)v = −ρ ∇Π − ∇Φ + 4πρ (B · ∇)B,
∂B = ∇ × (v × B),
∂t

(5.2) (5.3)

where ρ, P, v are the fluid density, pressure and velocity, Φ is the gravitational

potential due to the central compact object, and B2
Π≡P+ 8π

(5.4)

is the total pressure. The magnetic field B also satisfies the equation ∇ · B = 0.

We adopt the cylindrical coordinates (r, φ, z) which are centered on the

central object and have the z-axis in the direction perpendicular to the disc

plane. The unperturbed background disc has a velocity field v = rΩ(r)φˆ,

and the magnetic field may consist of both toroidal and vertical components

B = Bφ(r)φˆ + Bz(r)zˆ. The gravitational acceleration in radial direction is defined

as

g = dΦ dr

(5.5)

so that −∇Φ = −grˆ and g = rΩ2K > 0, where ΩK is the angular frequency for a test

mass (the Keplerian frequency). Thus the radial force balance equation reads

ρg = ρrΩ2 − dΠ − B2φ . dr 4πr

(5.6)

84

To investigate the dynamical properties of the flow, we perturb the MHD

Eqs. (7.13)-(7.15) by rewriting any physical variable f as f + δ f with |δ f | ≪ | f |.

Since the unperturbed state is axisymmetric and steady, we consider all per-

turbation variables having the form δ f ∝ eimφ−iωt, with m being the azimuthal

mode number and ω the wave frequency. Note that the background flow and

magnetic field have no z-dependance and we assume that the perturbations also

have no z-dependance. The resulting linearized perturbation equations contain

the variables δv, δρ, δP, δΠ and δB. For mathematical convenience, we define a

new variable

δh = δΠ = δP + B · δB . ρ ρ 4πρ

(5.7)

Moreover, using ∆v = δv + ξ · ∇v = dξ/dt = −iωξ + (v · ∇)ξ, we find that the

Eulerian perturbation δv is related to the Lagrangian displacement vector ξ by

δv = −iω˜ ξ − rΩ′ξrφˆ, where prime denotes radial derivative and

ω˜ = ω − mΩ,

(5.8)

is the Doppler-shifted wave frequency. In the next two subsections, we will combine the perturbations equations into two first-order differential equations (ODEs) in terms of ξr and δh for the magnetosphere region and the disc region, respectively.

5.2.1 The magnetosphere
In the magnetosphere region (r < rin), the flow is assumed to be incompressible and have uniform density (ρ is constant). For this particular magnetized fluid system, the detailed linearized perturbation equations have been given in Fu &

85

Lai (2011b). Here we just display the final two ODEs for ξr and δh:

dξr dr

=

A11ξr

+

A12δh,

(5.9)

dδh dr = A21ξr + A22δh,

(5.10)

where

A11

=

1 −
r

ω˜ 2

− 2mω˜ Ω + m2ω2Aφ ω˜ 2 − m2ω2Aφ

,

m2

A12 =

, r2

A21

=

ω˜ 2

−

m2ω2Aφ

−

2rΩ dΩ dr

+

d 2

ln

Bφ

−

1

d ln r

ω2Aφ

−

4

(ω˜ Ω + mω2Aφ)2 (ω˜ 2 − m2ω2Aφ)

,

A22

=

2m r

ω˜ Ω + mω2Aφ , ω˜ 2 − m2ω2Aφ

(5.11) (5.12) (5.13) (5.14)

and ωAφ ≡ vAφ/r = Bφ/(r 4πρ) is the toroidal Alfve´n frequency. Note that al-

though we assume low plasma density for the magnetosphere, we still require

that density is not too low in order to keep Alfve´n speed far less than the speed

of light. Otherwise, our original MHD equations break down (e.g., Lovelace,

Romanova & Newman 2010). Also note that Eqs. (5.9) and (5.10) are the same

as Eqs. (119) and (120) (derived for a purely toroidal magnetic field) of §83 in

Chandrasekhar (1961). The vertical magnetic field Bz does not appear in our

perturbation equations because we assumed kz = 0, i.e., vertical field lines are

not perturbed.

Defining σ2 = ω˜ 2 − m2ω2Aφ, we can further combine Eqs. (5.9) and (5.10) into one single equation

ξr′′

+

d dr

ln(r3σ2)ξr′

+

1

− m2 r2

ξr

=

0.

(5.15)

Two special cases are of interest:

86

(i) For systems with Bφ = 0, Eq. (5.15) reduces to

ξr′′

+

3 ( r

−

2mΩ′ ω˜

)ξr′

+

1

− m2 r2 ξr

=

0,

(5.16)

or equivalently

W ′′

+

W′ r

−

m2 r2

rd 1−
mω˜ dr

κ2 2Ω

W = 0,

(5.17)

where W = rδvr and

κ = 2Ω d (r2Ω) 1/2 r dr

(5.18)

is the radial epicyclic frequency. Note that Eq. (5.17) recovers Eq. (12) in Tsang

& Lai (2009b). For either uniform rotation profile [Ω = const, κ = 2Ω, thus

(κ2/2Ω)′ = 0] or uniform angular momentum profile (Ω ∝ r−2, thus κ = 0), it has

the same simple solution

W ∝ rm or δvr ∝ rm−1.

(5.19)

The corresponding solution for ξr is

ξr ∝ rm−1

if

Ω = const;

rm−1 ξr ∝ ω − mΩ(r)

if

Ω ∝ r−2.

(5.20)

(ii) In the special case of Ω = const and Bφ ∝ r so that ωAφ ∝ Bφ/r √ρ is constant (note that we have assumed constant density in the magnetosphere), [ln(r3σ2)]′ in Eq. (5.15) reduces to [ln(r3ω˜ 2)]′. Therefore magnetic fields in Eq. (5.15) (which enters through term σ2) completely disappear and the equation is exactly the same as the one for Ω = const and Bφ = 0. So the solutions for the magnetosphere region in this case are also ξr ∝ rm−1, δh ∝ rm and the relation between these two wavefunctions can be obtained by substituting them back into Eqs. (5.9)-(5.10). Our calculations will focus on this particular magnetosphere setup.

87

5.2.2 The disc

In the disc region (r > rin), we assume the flow is barotropic so that δP = c2sδρ with cs = dP/dρ being the sound speed, which will be parameterized by a constant cˆs ≡ cs/(rΩ) in our computation. For concreteness, we will also assume a power-law disc density profile ρ ∝ r−p with p being a constant. For simplicity, we take the disc toroidal magnetic field to be uniformly distributed, i.e., Bφ = const. From Eq. (5.6), we obtain the disc rotation profile as

Ω(r) = ΩK(r) . 1 + pcˆ2s

(5.21)

The perturbation equations in component form can also be presented in the

same form as Eqs. (5.9) and (5.10) with the detailed expressions of A11, A12, A21

and A22 given in Fu & Lai (2011a). For a non-magnetic disc, these equations

reduce to

dξr = dr

−

2mΩ rω˜

−

d ln(rρ) dr

ξr +

m2 1 r2ω˜ 2 − c2s

δh,

dδh dr

=

(ω˜ 2

−

κ2)ξr

+

2mΩ δh. rω˜

(5.22) (5.23)

5.2.3 The interface

In the equilibrium state, pressure balance at the interface reads

Pm + Pbz + Pbφ = Pd,

(5.24)

where Pm, Pbz and Pbφ are the magnetosphere gas pressure, magnetic pressure of Bz and magnetic pressure of Bφ just inside the interface, respectively, while Pd is the disc total pressure just outside the interface. This equality imposes an upper limit on the strength of the magnetosphere toroidal B field, Pbφ/Pd < 1. Defining

88

b = (ωAφ)|rin/Ωd and µ = (ρd − ρm)/(ρd + ρm), where Ωd is the disc rotation rate at the interface, ρm and ρd are the fluid densities of the magnetosphere and the disc

at the interface, respectively, then the above inequality becomes

Pbφ Pd

≃

(1 − µ)b2 (1 + µ)2cˆ2s

(5.25)

Note that in the magnetosphere, magnetic pressure dominates over gas pressure

(Pbz + Pbφ ≫ Pm). So Pbφ/Pd ≃ Pbφ/(Pbφ + Pbz) approximately characterizes the

relative strength of the toroidal field compared to the vertical field.

In the perturbed state, we demand that both the radial Lagrangian displacement ξr and the Lagrangian perturbation of total pressure ∆Π = δΠ + ξrΠ′ be continuous across the interface. The latter of these requirements yields

ρm δh + ξrr(Ω2 − Ω2K − ω2Aφ) rin− = ρd δh − ξrr pcˆ2s Ω2 rin+ .

(5.26)

In the magnetosphere, we have the analytical solutions for δh (∝ rm) and ξr (∝

rm−1) (see Sec. 5.2.1), which are related via

r δh =
m

ω˜ 2 + 2ω˜ Ω + m(2 − m)ω2Aφ

ξr.

(5.27)

We substitute this relation into the left hand side of Eq. (7.59) and consolidate

it with another requirement (continuity of ξr) to obtain the matching condition

across the interface

(ω

−

mΩm)[ω − m

(m

−

2)Ωm]

+

Ω2m

=

1+µ (1 − µ)rin

δh ξr

+

1 − 2µpcˆ2s + (m − 1)b2 1−µ

Ω2d,

(5.28)

where δh and ξr are disc solutions at the interface, Ωm and Ωd are the rotation

rates of the magnetosphere and the disc at the interface, respectively. In the

case of Bφ = 0 in the magnetosphere, this matching condition reduces to Eq. (22)

in Tsang & Lai (2009b). Note that Eq. (5.28) has no explicit dependence on the

disc magnetic field. This has to do with the fact that we chose an uniformly

distributed disc B field.

89

5.3 Interface modes: Newtonian potential
Li & Narayan (2004) and Tsang & Lai (2009b), who considered incompressible and compressible discs respectively, have shown that interface modes of magnetosphere-disc boundary can be subject to Rayleigh-Taylor and/or KevinHelmholtz instabilities, depending on the density contrast and velocity shear across the interface. In order to investigate the effects of magnetosphere toroidal magnetic fields on these instabilities, in this section we assume the disc region to be magnetic-field free and employ the standard shooting method (Press et al. 1992) to solve Eqs. (5.22)-(5.23) using the matching condition Eq. (5.28) as the inner boundary condition at rin and outgoing wave condition (Tsang & Lai 2009b) at the outer boundary of the disc (rout). The complex mode frequency ω can then be determined as an eigenvalue of the system. Note that in general we also need to apply shooting method in the magnetosphere region (see Appendix A), but for the specific magnetic field profile (and uniform rotation) of the magnetosphere we chose in this study, the perturbation equations in the magnetosphere have simple analytic solutions (see §2.2) and the effect of magnetosphere on the interface modes is embodied by the interface matching condition (i.e, inner boundary condition for disc perturbation equations).
In our computation, we use p = 1.5, cˆs = 0.15, and put the outer boundary at rout = 2.85rin. We vary the strength of magnetosphere toroidal magnetic field b (therefore Pbφ/Pd; Eq. [5.25]) for different sets of (Ωm/Ωd, ρm/ρd) to see how the eigenvalue ω will be modified. Throughout this section, we use Newtonian potential Φ = −GM/r, so that the Keplerian frequency is ΩK = (GM/r3)1/2, and the epicyclic frequency κ equals Ω.
90

Figure 5.2: The wavefunctions for the unstable interface mode with m = 4, ρm/ρd = 1/99, Ωm/Ωd = 1 in a magnetosphere-disc system without (left) and with (right) magnetosphere toroidal magnetic field, respectively. The vertical long-dashed lines mark the position of the interface which separates magnetosphere region and disc region. The solid and dotted lines represent the real and imaginary parts of the wavefunctions, respectively.
5.3.1 Numerical solutions
Fig. 7.5 shows two example wavefunctions for the unstable interface mode. The left panel depicts the case with no toroidal magnetic field in the magnetosphere (but keep in mind that there is always finite vertical B field) while in the right panel the magnetosphere has a non-zero Bφ. Since there is no velocity shear across the boundary in both cases (Ωm/Ωd = 1), the instability is of RayleighTaylor type. We see that the addition of toroidal field has a small effect on the eigenfunction. However, as the numbers in the figure indicate, the growth rate of the unstable mode is reduced by more than 50% (from 0.256Ωd to 0.11Ωd) even though the the toroidal magnetic field pressure is only 10% of the disc gas pressure (both evaluated at the interface), i.e. even though the toroidal field is
91

Figure 5.3: Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at the interface) for unstable interface modes as a function of magnetosphere toroidal magnetic field strength. The top and bottom panels are for the cases without and with velocity shear at the interface, respectively, while the left and right panels depict models with the magnetosphere to disc density ratio being 1/99 and 1/9, respectively. Different line types are associated with different values of m, with the solid lines for m = 1, dotted lines for m = 2, short-dashed lines for m = 3, long-dashed lines for m = 4, and dot-dashed lines for m = 5.
92

much weaker than vertical field in the magnetosphere.
The suppressing effect of toroidal magnetic field on the instability can be more easily seen in Fig. 7.3, where we plot the real and imaginary parts of the eigen-frequencies for the interface mode as a function of Bφ. We consider four different sets of (Ωm/Ωd, ρm/ρd) such that the panels in the same row have the same velocity shear but different density contrast while the panels in the same column have the same density contrast but different velocity shear. The modes in the top two panels are subject to Rayleigh-Taylor type instability as there are no velocity shear at the interface. The bottom two panels, however, have non-zero velocity shear, thus are subject to instabilities of both RayleighTaylor and Kelvin-Helmholtz types. In all cases, we observe that the inclusion of toroidal magnetic fields can significantly diminish the growth rate of the unstable modes, even completely kill the instability. The critical Bφ for absolute shutdown of unstable modes depends on the detailed interface parameters (sound speed in the disc, velocity shear, density contrast, etc.). In the case with the most unstable modes (bottom-right panel) that we have calculated, the required Pbφ/Pd to fully turn off the mode growth is close to one. This means that the toroidal magnetic field needs to dominate over the vertical field in the magnetosphere (remember that Pbφ/Pd ≃ Pbφ/(Pbφ + Pbz)). The mode frequencies (real part of ω), on the other hand, are barely affected by Bφ (at least when the growth rate is not very small). By comparing panels in the same column, we see that increasing velocity shear at the interface leads to larger growth rates, which is not surprising as more velocity shear means more free energy that the system can tap on to drive the instability. Comparing panels in the same row reveals an interesting features. We found that larger density contrasts leads to smaller growth rates, which apparently contradicts the standard Rayleigh-
93

Taylor instability result. In particular, the m = 1 mode is totally stable when ρm/ρd = 1/99 while becomes unstable when ρm/ρd = 1/9 and Bφ does not affect the m = 1 mode at all. Comparing different lines in each individual panel shows that modes with larger m generally have larger growth rates, but they are also more easily stabilized by the toroidal field.
5.3.2 Discussion of results
The zero effect of Bφ on the m = 1 mode (whether it is stable or unstable) can be easily understood by looking at the matching condition Eq. (5.28). This equation (i.e., inner boundary condition for disc equations) is the only place that Bφ affects the eigenvalue problem, and the magnetic effect comes in via the last term on the right hand side of the equation [(m−1)b2Ω2d]. When m = 1, this term vanishes so that the equation is the same as in the case with zero Bφ.
To explain the other features observed in Fig. 7.3, we carry out a local analysis of the effects of magnetic fields on Rayleigh-Taylor/Kelvin-Helmholtz instability in a plane-parallel flow (see Appendix B). From the last term in the final expression of wave frequency (Eq. [B.1]), we see that magnetic tension (k · B) provides a suppressing force against Rayleigh-Taylor/Kelvin-Helmholtz instability of a two-layered fluid system. This is because extra work needs to be done in order to increase the boundary layer deformation (i.e., growing perturbation). Our cylindrical magnetosphere-disc model resembles the simple parallel two-layered flow in that our wavenumber in azimuthal direction m/r acts like k in the x direction and our Bφ also lies along the direction of the wave vector. Thus, the same stabilizing mechanism also applies in our system, i.e., the
94

tension force of the toroidal magnetic field suppresses the magnetosphere-disc interface instability.

Qualitatively, the growth rate of interface modes in our system is deter-

mined by combined effects of four factors: density contrast (ρm/ρd), velocity shear across the interface (Ωm/Ωd), degree of differential rotation in the disc, and

toroidal field of the magnetosphere. The first two tend to enhance instability,

while the latter two tend to suppress instability. In different limiting cases, ap-

proximate analytic expressions for the mode growth rate can be derived: Li &

Narayan (2004) considered a case where the disc is incompressible with con-

stant density and came up with their Eq. (45), while Tsang & Lai (2009b) stud-

ied a compressible disc whose density is much larger than the magnetosphere

density, and summarized the results in their Eqs. (26) and (28). Although for

the generic conditions studied in this chapter, an analytical expression cannot

be rigorously derived, we can write the mode growth rate schematically as fol-

lows:

ωi ∼ ω2RT + ω2KH − ω2vort − ω2b

(5.29)

with

ω2RT ≃ 2(1 + µ)mΩ2eff,d − 2(1 − µ)mΩ2eff,m,

(5.30)

ω2KH ≃ (1 − µ2)m2(Ωd − Ωm)2 + 2(1 − µ2)m(Ωd − Ωm)(ζd − ζm),

(5.31)

ω2vort ≃ (ζd − ζm)2 + 2µ(ζd2 − ζm2 ) + µ2(ζd + ζm)2,

ω2b

≃

m2 B2φ 4πr2(ρd +

ρm)

≃

(1

−

µ)m2ω2Aφ/2,

(5.32) (5.33)

where rΩ2eff = dΦ/dr − rΩ2 is the effective gravitational acceleration, ζ = κ2/(2Ω)

is the fluid vorticity and as before the subscripts d and m denote disc and mag-

netosphere, respectively. Several features can be noted in Eq. (5.29):

95

(i) The first two terms under the square root are the Rayleigh-Taylor and Kelvin-Helmholtz types of destabilizing factors while the last two terms characterize the stabilizing effects due to flow vorticity and toroidal magnetic field in the magnetosphere, respectively. The suppressing effect of finite fluid vorticity has been discussed in both Li & Narayan (2004) and Tsang & Lai (2009b). We see that the various terms have different dependences on the density contrast µ. As µ increases (µ = 0 when ρm/ρd = 1, and µ = 1 when ρm/ρd = 0), ω2RT increases, making the system more Rayleigh-Taylor unstable; at the same time, ω2KH decreases, making the system less Kelvin-Helmholtz unstable. In addition, ω2vort also becomes larger, leading to stronger suppressing effect. This explains why in Fig. 7.3, the mode with ρm/ρd = 1/99 is less unstable than the one with ρm/ρd = 1/9.
(ii) Since Ω2eff,d = (1/r)(dΦ/dr) − Ω2 = −(c2s/ρ)(dρ/dr), when cs is too small (i.e., small effective gravity in the disc), the destabilizing terms would not be able to compete with the stabilizing terms, resulting in a stable system. Thus, a sufficiently high disc sound speed is needed to attain the interface instability (see Tsang & Lai 2009b).
(iii) When there is no Bφ, we see that ω2RT and ω2KH both depend on m while ω2vort does not. Thus, modes with higher m tend to be more unstable, as shown in Fig. 7.3 (the Pbφ/Pd = 0 case; see also Li & Narayan 2004 and Tsang & Lai 2009b). The m = 1 mode would be even less unstable if the density contrast is too big [see (i) above]. This explains why in Fig. 7.3 the m = 1 mode with ρm/ρd = 1/99 is stable while the one with ρm/ρd = 1/9 is unstable.
(iv) When Bφ 0, besides ω2b, the magnetic field strength also appears in ω2RT where Ω2eff,m contains a term that is proportional to −B2φ/4πρmr2 (note the mi-
96

nus sign). This shows that the toroidal field in the magnetosphere plays two different roles in determining the stability of system. On the one hand, the magnetic tension resists perturbation growth, thus suppressing any instability; on the other hand, the magnetic force increases the effective gravity (pointing towards the center) of the background flow, thus promoting Rayleigh-Taylor type instability. The former effect is proportional to m while the latter to m2. Hence, in general the suppressing effect is more important. This is consistent with the fact that in Fig. 7.3 the growth rates of modes with higher m decrease more rapidly as Bφ increases. As noted above, Bφ does not affect the m = 1 mode because the last term in Eq. (5.28) vanishes for m = 1. Now we have a better understanding of what this means physically: For the m = 1 mode, the aforementioned two opposing effects associated with Bφ happen to cancel each other. This exact cancellation, however, cannot be captured by the approximate expression Eq. (5.29).
Overall, we see that Eq. (5.29), though schematic, is quite useful in explaining most of the numerical results of this section. Note that various flow parameters, such as the density contrast µ, azimuthal mode number m, toroidal magnetic field Bφ and disc vorticity ζd, appear in more than one of the four terms in this equation. Thus, they are associated with both the stabilizing and destabilizing effects. The only exception is perhaps the velocity shear at the interface, which always facilitates the interface instability.
Finally, we note that the interface instability associated with magnetospheredisc boundary studied in this chapter is qualitatively different from those of Lubow & Spruit (1995) and Spruit, Stehle & Papaloizou (1995), who carried out local analysis of a thin rotating disc threaded by a large-scale poloidal field.
97

Lovelace, Romanova & Newman (2010) considered similar interface instability as in our chapter, but they focused on small-scale modes and included only vertical field in the magnetosphere.
5.4 Interface modes and Discoseismic modes: Pseudo-
Newtonian potential
In this section, we consider how the global oscillation modes in our discmagnetosphere system (Fig. 1) are modified by general relativistic (GR) effects. In particular, GR changes the disc rotation profile ΩK(r) and makes the radial epicyclic frequency κ(r) a non-monotonic function of r: As r decreases, κ first increases, attains a maximum value and then falls to zero at the Innermost Stable Circular Orbit (ISCO), rISCO = 6GM/c2 (for a non-spinning compact object). If the inner disc boundary is close to rISCO, then this non-monotonic κ profile can have two consequences: (1) It significantly reduces disc vorticity (κ2/2Ω) near the inner disc boundary, therefore helps the unstable interface modes grow as there are less disc vorticity suppressing effects (see the discussion in Section 5.3.2; see also Tsang & Lai 2009b); (2) It could lead to a vortensity (κ2/2Ωρ a.k.a potential vorticity) profile that has positive gradient in the inner disc region. As found in our previous studies (Tsang & Lai 2008; Lai & Tsang 2009; Fu & Lai 2011a), positive vortensity gradient renders the disc inertial-acoustic modes (p-modes) unstable. In Section 5.3, we chose a density profile ρ ∝ r−3/2 so that the vortensity profile (with the Newtonian potential) is completely flat, and thus we only found unstable (overstable) interface modes. In this section, with the GR effect included, both the interface modes and disc p-modes can be unstable and they
98

co-exist in our disc-magnetoshere system. The main goal of this section is then to study these two modes and the effects of the magnetosphere toroidal B field, disc toroidal B field and inner disc boundary location.

We employ the pseudo-Newtonian potential (Paczynski & Wiita 1980) to

mimic the GR effect:

GM Φ=− ,
r − rS

(5.34)

where rS = 2GM/c2 is the Schwarzschild radius. The corresponding Keplerian

rotation frequency and epicyclic frequency are

ΩK = κ=Ω

GM 1 ,
r r − rS r − 3rS . r − rS

(5.35) (5.36)

Except for Section 5.4.3, we will take the magnetosphere-disc interface to be

located at the ISCO, rin = rISCO = 3rS.

5.4.1 Non-magnetic discs
We first consider the case where the disc outside the magnetosphere is nonmagnetic. We solve the same equations and use the same interface boundary condition as in Section 3 except that we replace the Newtonian potential with Pseudo-Newtonian potential.
Figure 5.4 shows the complex eigenfrequencies of disc inertial-acoustic modes and interface modes as a function of the dimensionless magnetosphere toroidal field strength Pbφ/Pd. The other parameters are indicated in the figure. As noted above, because of the GR effect, the disc vortensity profile has a positive slope near the ISCO (see Fig. 5.5). The disc inertial-acoustic mode
99

Figure 5.4: Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable inertial-acoustic disc modes (left panels) and unstable interface modes (right panels) as a function of the magnetosphere toroidal field strength. Different line types represent different azimuthal mode number m. The sound speed profile, density contrast and velocity shear across the interface are the same as in the bottom-left panel of Fig. 7.3. The only difference from Fig. 7.3 is that GR effect (using pseudo-Newtonian potential) is included in the disc rotation profile, which renders both types of modes unstable. On the right panles, the vertical lines (for m = 2 and 3) mark the point when corotation resonance moves into the flow (see Fig. 5.5).
(p-mode) of lowest radial order (for a given m) has Re(ω) < mΩd, and the coro-
tation resonance lies in the disc (see the right panels of Fig. 5.6 for an example).
Corotational wave absorption then makes these modes unstable. We see from
Fig. 5.4 that as long as there is a magnetosphere to serve as an inner boundary
condition for the disc mode, the mode frquency and growth rate depend rather
weakly on Pbφ/Pd.
The interface mode is also strongly influenced by the the GR effect. Com-
paring the right panels of Fig. 5.4 with the bottom-left panel of Fig. 7.3, we see
100

Figure 5.5: Wave propagation diagram (upper panel) for both interface modes and disc inertial-acoustic modes accompanied by the disc vortensity profile (lower panel). In the upper panel, the three thick solid curves depict the disc rotation profile Ω and Ω ± κ/m, where m is the azimuthal mode number and κ is the radial epicyclic frequency. Note that since κ(rISCO) = 0 these three curves join each other at the Innermost Stable Circular Orbit (which is also the disc inner boundary rin in our setup of Sections 5.4.1 & 5.4.2). The dashed horizontal line marks the point when the corotation radius [where Re(ω)/m = Ω] happens to be exactly at the inner disc boundary. For modes with Re(ω)/m above this line, there is no corotation in the flow; for any modes with Re(ω)/m below this line, the corotation exists in the flow (rc > rin). Three types of wave modes (labeled as I, II and III) are shown in this diagram where the horizontal solid lines indicate the evanescent zones and the wavy curves denote the wave propagation zones. Note that for both Type I and II waves, there is also a tiny wave zone (hardly visible in the figure) in the disc region just outside the ISCO. Type I and II waves are both interface modes while type III waves are disc inertial-acoustic modes. Note that the magnetosphere region (r < rISCO) and the region between the inner and outer Lindblad resonances [where Re(ω)/m = Ω ± κ/m] are wave evanescent zones.
101

Figure 5.6: Wave functions for the three types of wave modes depicted in Fig. 5.5: the interface mode with no corotation resonance in the flow (Type I, left panels), interface mode with corotation resonance in the flow (Type II, middle panles) and disc inertial-acoustic mode (Type III, right panels). The azimuthal mode number is m = 3 and all other parameters are the same in Fig. 5.4. The solid and dotted lines show the real and imaginary parts, respectively. In the middle and right panels, the vertical dot-dashed lines represent the location of corotation resonance.
that for Pbφ = 0, the mode growth rate (for a given m) is larger when the GR effect is included. Again, this is because in GR, the disc vorticity is smaller near the ISCO, leading to less rotational suppression of the RT instability (see Section 3.2). As Pbφ/Pd increases, the interface mode growth rate first decreases (just as in the Newtonian case) due to magnetic tension, then starts to increase beyond certain critical value of Pbφ/Pd (for a given m). This behaviour arises from the effect of corotation resonance (see Fig. 5.5): For small Pbφ/Pd, the interface mode has Re(ω) > mΩd, thus no corotation resonance exists in the flow; as Pbφ/Pd becomes larger, the mode frequency drops below mΩd, and corotation resonance comes into play, which overwhelms the suppression effect of the magnetic tension (see the left and middle panels of Fig. 5.6 for two examples of mode wavefunctions). The growth rate of such interface mode (labeled as Type II in Fig. 5.5) is much larger than the coresponding disc p-mode because
102

the corotation resonance lies much closer to rin for the interface mode than for the p-mode – this gives rise to much stronger wave absorption at the corotation resonance.
5.4.2 Magnetic discs
Now we consider the case in which the disc outside the magnetosphere has a finite toroidal magnetic field. For simplicity, we take the limit that the density in the magnetosphere ρm goes to zero. In this case, the interface boundary condition Eq. (5.28) reduces to ∆Π = 0, i.e., the Lagrangian perturbation of total disc pressure equals zero at rin. With this boundary condition, the setup is the same as in Fu & Lai (2011a) and we only need to consider the dynamical equations for the disc.
Again, two types of oscillation modes exist in our system (see Fig. 5.7): the disc inertial-acoustic modes and interface modes. The effect of the disc magnetic field on the inertial-acoustic modes has already been studied in detail by Fu & Lai (2011a). We see from the left panels of Fig. 5.7 that the mode frequnecy is only slightly modified by the disc Bφ, but the growth rate is reduced so that the mode become stable even for modest (sub-thermal) disc toroidal fields.
The interface mode (the right panels of Fig. 5.7) was not considered by Fu & Lai (2011a). We see that as the disc magnetic field increases, the mode frequency varies modestly (about 10%) for a wide range of Bφ, but the growth rate changes more significantly. The dependence of the growth rare as a function of disc Bφ can be understood as follows: (i) For Bφ = 0, RT instability drives the mode growth; (ii) As Bφ increases, the magnetic tension tends to suppress the growth;
103

Figure 5.7: Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable inertial-acoustic disc modes (left panels) and unstable interface modes (right panels) as a function of the disc toroidal field strength. The x-axis specifies the ratio of the Alfve´n velocity to disc rotation velocity at the inner boundary rin. The magnetosphere density is set to zero. The disc parameters are rin = rISCO, m = 2, cs = 0.1rΩ and ρ ∝ r−1. For the interface mode (right panels), the two dotted vertical lines from left to right mark the entering of the outer magnetic resonance and inner magnetic resonance into the flow, respectively.
(iii) As Bφ increases, the real mode frequency decreases. When Bφ exceeds some critical value, corotation resonance appears in the flow (disc), and wave absorption at corotation then enhances the growth rate. But as Bφ increases further, the corotation radius lies at a larger distance from rin, thus the corotational effect becomes less important (because of the large evenacent zone separating rin and the corotation radius) and the mode growth rate decreases again also because of magnetic tension.
Concerning (iii) above, as noted in Fu & Lai (2011a), in the presence of the disc toroidal magnetic field, the corotation resonance (where ω˜ = 0) is split into
104

the inner/outer magnetic resonances, where

ω˜ = ±mωAφ,

(5.37)

where ωAφ = vAφ/r = Bφ/(r 4πρ) is the toroidal Alfve´n frequency of the disc. When these magnetic resonances exist in the flow (disc), wave absorption comes into play in the mode growth. Note that the signs of wave absorption at the two magnetic resonances are different. As Re(ω) decreases (with increasing Bφ), the outer magnetic resonance (where ω˜ = mωAφ) first enters the flow, causing the mode growth rate to increase. When the inner magnetic resonance (where ω˜ = −mωAφ) enters the flow, the wave absorptions at the two magnetic resonances partially cancel, and the mode growth rate starts to decrease again (see the right panels of Fig. 5.7). In this case, the term “corotation resonance” simply refers to the combined effect of two magnetic resonances.

5.4.3 Effects of different inner disc radii
Here we consider the same setup as in Sections 5.4.1-5.4.2 except that the inner disc boundary lies outside the ISCO. The motivation for considering rin > rISCO is that in real accreting NS systems, the magnetosphere radius may well be outside the ISCO (e.g., in accreting millisecond X-ray pulsars with surface magnetic field of 108−9 G, the Alfve´n radius is about 1.5-2 stellar radii). Also, in the “transitional state” (when HFQPOs are observed) of BH X-ray binaries, the thin thermal disc may be truncated at a radius slightly larger than the ISCO (e.g., Done et al. 2007; Oda et al. 2010).
As in Section 5.4.2, we assume that the magnetosphere has a negligible density. Figure 5.8 shows an exmaple (for a non-magnetic disc) of how the (com-
105

Figure 5.8: The real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable inertial-acoustic modes (left panels) and unstable interface modes (right panels) as a function of the inner disc radius rin (in units of rISCO). The magnetosphere inside rin has zero density and the disc toroidal magnetic field is set to zero. The other disc parameters are m = 2, cs = 0.1rΩ and ρ ∝ r−1.
plex) mode frequencies depend on rin. When expressed in units of Ωd, the real frquencies of both interface and inertial-acoustic modes are only modestly affected by rin/rISCO, but the growth rates decrease rapidly with increasing rin/rISCO. For the inertial-acoustic mode, this arises from the reduced vortensity slope at distances further away from the ISCO, which results in smaller wave absorption at corotation. For the interface mode, the larger vorticity just beyond rin leads to a stronger rotational suppression of the RT instability, thus a smaller mode growth rate.
Figure 5.9 shows how the complex frequency of the interface mode depends on the disc magnetic field Bφ when the inner disc radius is set to 1.5rISCO. The behaviour of the mode growth rate as a function of (vAφ/rΩ)in is similar to the
106

Figure 5.9: Real and imaginary parts of the wave frequencies (in units of Ωd, the disc rotation rate at rin) for unstable interface modes as a function of the disc toroidal field strength, with the inner disc boundary located at rin = 1.5rISCO. The x-axis specifies the ratio of the Alfve´n velocity to disc rotation velocity at rin. The magnetosphere density is set to zero. The disc parameters are rin = rISCO, m = 2, cs = 0.1rΩ and ρ ∝ r−1. The two dotted vertical lines from left to right mark the entering of outer magnetic resonance and inner magnetic resonance into the flow, respectively.
right panels of Fig. 5.7. Here, the mode growth rate is small when Bφ = 0. So when the outer magnetic resonance enters the flow, wave absoprtion dramatically increases the mode growth rate.
5.5 Summary
In this chapter we have studied the non-axisymmetric MHD modes and instabilities in a 2D model of magnetosphere-disc systems (see Fig. 1), as may be
107

realized in accreting neutron star or black-hole X-ray binaries (see Section 5.1). We have examined various physical effects and parameters that can influenece the global modes in the system, including the density and magnetic field of the magnetosphere, the velocity contrast across the magnetosphere-disc interface, the rotation profile (Newtonian vs GR), temperature and magnetic field of the disc. We restrict to modes that do not have vertical structure, but otherwise our calculations include all possible instabilities and global oscillations associated with the interface and the disc. We highlight several key findings of this chapter below.
5.5.1 Interface instabilities in a rotating, magnetized system
This chapter includes a comprehensive study of the large-scale Rayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) instabilities associated with the interface of a rotating, magnetized system. RT and KH instabilities have been studied intensively in plane-parallel flows through both theoretical analysis and laboratory experiments (e.g. Chandrasekhar 1961; Drazin & Reid 1981), and have found applications in various astrophysical and space environments. But few papers have focused on rotating systems (e.g., Spruit et al. 1995; Lovelace et al. 2009, 2010). Our study generalizes previous works by Li & Narayan (2004) and Tsang & Lai (2009b) by considering compressible fluid, magnetic field and rotation profile of the disc, as well as generic field (poloidal and toroidal) configuration of the magnetosphere.
As in plane-parallel flows, the interface modes are mainly driven unstable by the RT and KH instabilities. Toroidal magnetic field tends to suppress the
108

instabilities through magnetic tension. The magnetic field also indirectly affects the RT instability by modifying the effective gravity. Except for the m = 1 mode, for which the two opposite effects of Bφ cancel (as conjectured by Li & Narayan 2004), we find that increasing Bφ generally tends to reduce the growth rates of the interface modes. Differential rotation (with finite vorticity) also tends to suppress the interface instability. To overcome this suppression effect, the disc must have sufficiently large temperature (sound speed) (see Tsang & Lai 2009b). General relativity (GR) can significantly affects the growth rate of interface modes because the disc rotation near the ISCO has smaller vorticity in GR than in Newtonian theory.
Another qualitatively new finding of this chapter is that corotation resonance can significantly influence the interface instabilities. As the toroidal field strength in the magnetosphere or in the disc increases, the real frequency of the interface mode falls below mΩd (where Ωd is the disc rotation rate at the interface), and corotation resonance (or its generalization to magnetic resonances) appears in the disc. Wave absorption at corotation can then significantly change the interface mode growth rate (see Figs. 5.4, 5.7, 5.9).
5.5.2 Inertial-Acoustic Modes of relativistic Disc
Our model system (Fig. 1) also acomodates interial-acoustic modes (p-modes) of relativistic discs. These modes are driven unstable primarily by wave absorption at corotation resonance. The magnetosphere-disc interface naturally serves as the inner boundary for the disc. Our study in this chapter complements our previous works (Lai & Tsang 2009; Tsang & Lai 2009c; Fu & Lai 2011a) by
109

properly treating the inner boundary condition for disc oscillations. Our result shows that the magnetosphere behaves as a robust “reflector” for spiral waves in the disc: The p-mode frequency and growth rate do not depend sensitively on the property (density, magnetci field) of the magnetosphere. In agreement with Fu & Lai (2011a), we find that a modest disc toroidal field tends to reduce the growth rate of disc p-modes.
110

CHAPTER 6 PAPALOIZOU-PRINGLE INSTABILITY OF MAGNETIZED ACCRETION
TORI
6.1 Introduction
Hot accretion tori around a compact object are known to be susceptible to a global hydrodynamical instability, the so-called Papaloizou-Pringle (PP) instability, arising from the interaction of non-axisymmetric waves across the corotation radius, where the wave pattern speed matches the fluid rotation rate. However, accretion tori produced in various astrophysical situations (e.g., collapsars and neutron star binary mergers) are likely to be highly magnetized. In this chapter, we study the effect of magnetic fields on the PP instability in incompressible tori with various magnetic strengths and structures. In general, toroidal magnetic fields have significant effects on the PP instability: For thin tori (with the fractional width relative to the outer torus radius much less than unity), the instability is suppressed at large field strengths with the corresponding toroidal Alfve´n speed vAφ ∼> 0.2rΩ (where Ω is the flow rotation rate). For thicker tori (with the fractional width of order 0.4 or larger), which are hydrodynamically stable, the instability sets in for sufficiently strong magnetic fields (with vAφ ∼> 0.2rΩ). Our results suggest that highly magnetized accretion tori may be subjected to global instability even when it is stable against the usual magneto-rotational instability.
111

6.2 Equilibrium Models

As mentioned above, the PP instability operates in modes with no vertical structure (kz = 0). As such, the dynamics can be captured by height-averaged fluid equations. We consider a cylindrical shell (of finite width) of incompressible non-self-gravitating fluid, which is rotating differentially in the external gravitational field produced by a central compact object. We adopt the cylindrical coordinates (r, φ, z) with the z-axis being the rotation axis. The cylindrical shell is assumed to be infinitely long in the z-direction and threaded by magnetic fields. The fluid satisfies the ideal MHD equations:

∂u 1 1 ∂t + (u · ∇)u = −ρ ∇Π − ∇Φ + 4πρ(B · ∇)B,
∂B ∂t = ∇ × (u × B),
∇ · B = 0,

(6.1)
(6.2) (6.3)

∇ · u = 0.

(6.4)

Here, ρ is the constant fluid density, u the fluid velocity, B the magnetic field,

and Π = P + B2/8π the total pressure with P being the gas pressure. The gravi-

tational potential is Φ = −GM/r, where M is the mass of the central object. The background flow velocity u = rΩ(r)φˆ and magnetic field B = Bφ(r)φˆ + Bz(r)zˆ

also depends only on r. For convenience, we assume the flow has a power-law

rotation profile

Ω(r) ∝ rp.

(6.5)

The flow is confined between two boundaries (r1 ≤ r ≤ r2), where the gas pressure vanishes (P|r1, r2 = 0). Outside the fluid zone is a vacuum devoid of matter but maybe permeated with magnetic fields. In the equilibrium state, we

112

assume that the magnetic field is continuous across the fluid boundaries so that there is no surface electric current at r = r1, r2 (however, we allow for surface current to develop when the fluid is perturbed). We will consider two models of magnetic field structure.

6.2.1 Model (a)

In this model, we assume that there is an external current running vertically at small radii much inside r1, giving rise to Bφ(r) ∝ r−1 in the inner region (r < r1). There is no azimuthal current in this region, so Bz is constant. In the fluid zone, we adopt a power-law magnetic field profile

Bφ(r) ∝ rq, Bz(r) ∝ rs,

(6.6)

which means that both the azimuthal and vertical current densities are also of power-law form. Outside the fluid zone (r > r2), there is no current. Hence, Bφ(r) ∝ r−1, Bz(r) = const. The complete magnetic field profile is illustrated in the upper two panels of Fig. 6.1.

Integrating the radial equilibrium equation

1 ρ

dΠ dr

=

GM − r2

+

rΩ2

−

B2φ , 4πρr

gives the gas pressure profile

P ρ

=

GM r

+

r2Ω2 2p + 2

−

1 2

v2Az

−

1 2

1 1+
q

v2Aφ − C,

where

vAz = Bz/ 4πρ, vAφ = Bφ/ 4πρ

(6.7) (6.8) (6.9)

113

Figure 6.1: The two magnetic field profiles adopted in the equilibrium torus model.

are the Alfve´n velocities and C is the integration constant. The location of gas

pressure maximum is determined by

d (P/ρ)
dr

=

GM − r2

−

s v2Az r

−

(1

+

q) v2Aφ r

+

rΩ2

=

0,

(6.10)

which defines a reference radius r0:

GM r02

=

r0Ω20 −

s v2Az(r0) r0

− (1 + q) v2Aφ(r0) . r0

(6.11)

Let C ≡ λGM/r0 with λ being a constant and use Eq. (6.11) to substitute GM in

Eq. (6.8), we can rewrite Eq. (6.8) in the dimensionless form

P ρ

=

1 r

−

λ

+

r2p+2 2p + 2

+

sv2Az0

1 − +λ−

1 r2s

r 2s

+ (1 + q)v2Aφ0

1 − +λ−

1 r2q

r 2q

,

(6.12)

where

vAφ0 = vAφ(r0)/(r0Ω0), vAz0 = vAz(r0)/(r0Ω0).

(6.13)

114

Here and hereafter we will use units such that r0 = Ω0 = 1. Once we specify p, q, s, vAφ0, vAz0 and λ, we can determine the locations of the torus boundary by solving P = 0. However, there are several constraints on these parameters:

(i) dP/dr = 0 only guarantees the extremum of the P(r) profile. To ensure that we find a pressure maximum instead of minimum, we require d2P/dr2 < 0 at r = r0 = 1, which implies

2p + 3 − s(1 + 2s)v2Az0 − (1 + q)(1 + 2q)v2Aφ0 < 0.

(6.14)

This requirement reduces to p < −3/2 in the B = 0 limit.

(ii) Both sides of Eq. (6.11) need to be positive so that the gas pressure maxi-

mum exists. Thus

1 − sv2Az0 − (1 + q)v2Aφ0 > 0.

(6.15)

(iii) The maximum gas pressure Pmax must be positive. Thus, requiring the RHS of Eq. (6.12) to be positive at r = 1 gives

1 λ < 1 − v2Az0 − 2v2Aφ0

2p 2p

+ +

3 2

−

3 2

v2Az0

−

3v2Aφ0

,

provided that Eq. (6.15) is satisfied.

(6.16)

Figure 6.2 illustrates some examples of Model (a). We specify the values of p, q, s, vAφ0 and vAz0, then by varying λ, we obtain solutions for different torus thickness. For a given ∆r/r2, both r1 and r2 change when vAφ0 changes as a result of magnetic support in the torus. In the hydro limit (vAφ0 = vAz0 = 0), r1 approaches 0.5 as ∆r/r2 → 1. This feature is shown analytically in Pringle & King (2007). For a finite field strength, we see that r1 → 0 and r2 → ∞ as ∆r/r2 → 1. Note that for a relatively thin torus (∆r/r2 0.6), the differences of r2 and r1 between different field strengths are quite small. Since Bφ and Bz have similar

115

Figure 6.2: Some examples of Model (a) with a pure toroidal magnetic field (Bφ ∝ r in the fluid and Bz = 0) and constant angular momentum distribution (Ω ∝ r−2). The x-axis is the relative thickness of the torus with ∆r = r2 − r1 being the absolute width, and the y-axis shows the locations of the two boundaries. The different lines represent different values of vAφ0 = vAφ(r0)/(r0Ω0), as indicated. The horizontal line indicates the location of gas pressure maximum (r0 = 1).
116

effects on the equilibrium structure (see Eq. [6.12]), these features also apply to models with a finite vertical field. The special case of s = 0 and q = −1 is worth mentioning: In this case, the magnetic field is force-free and has no effect on the equilibrium structure 1.

6.2.2 Model (b)

The magnetic field profile in this case is shown in the bottom two panels of

Fig. 6.1. Compared with Model (a), the difference is that there is no vertical

current at small radius. Therefore the azimuthal field in the inner region (r < r1) is zero. In the fluid zone, we assume Bz ∝ r and Bφ(r) ∝ r − r12/r such that both azimuthal and vertical current densities are uniform. Following the same

procedure as in Section 6.2.1, we can derive the dimensionless expression for

gas pressure profile:

P ρ

=

1 r

−

λ

+

r2p+2 2p + 2

+

v2Az0

λ − 1 − 1r2 r2

+v2Aφ0

λ− 1 r

1 + r12 1 − r12

−

(1

1 − r12

)2

(r2

−

r12

−

2r12

ln

r)

.

(6.17)

Similarly, for a viable equilibrium model to exist, the model parameters must

satisfy the following requirements:

2p

+

3

−

3v2Az0

−

3 + r14 (1 − r12)2

v2Aφ0

<

0,

(6.18)

1

−

v2Az0

−

1 1

+ −

r12 r12

v2Aφ0

>

0,

1

λ<

1

−

v2Az0

−

v1+r12 2
1−r12 Aφ0

2p 2p

+ +

3 2

−

3 2

v2Az0

−

2 1

+ −

r12 r12

v2Aφ0

.

(6.19) (6.20)

1This is why in Curry & Pudritz (1996) the one-to-one mapping between r2/r1 and (r2 − r0)/r2 remains unchanged for different uniform vertical B field strengths.

117

6.3 MHD Equations for Perturbations

Assuming that the Eulerian perturbation of any physical variable f is of the

form δ f ∝ eimφ−iωt (with no dependance on z), the linearized perturbation equa-

tions are

1 ∂ im r ∂r (rδur) + r δuφ = 0

−iω˜ δur − 2Ωδuφ

=

1 −ρ

∂δΠ ∂r

+

imBφ 4πρr

δBr

−

Bφ 2πρr

δBφ

κ2 −iω˜ δuφ + 2Ω δur

=

im 1 − δΠ +
ρr 4πρ

∂1 +
∂r r

BφδBr

+

imBφ 4πρr

δBφ

−iω˜ δuz

=

imBφ 4πρr

δBz

+

1 4πρ

dBz dr

δBr

−iω˜ δBr

=

im

Bφ r

δur

−iω˜ δBφ

=

imBφ r

δuφ

−

r

d dr

Bφ r

dΩ δur + r dr δBr

−iω˜ δBz

=

imBφ r

δuz

−

dBz dr

δur,

(6.21) (6.22) (6.23) (6.24) (6.25) (6.26) (6.27)

where ω˜ = ω − mΩ is the wave frequency in the co-rotating frame and the radial

epicyclic frequency κ is given by

κ2 = 2Ω d (r2Ω) = 2(p + 2)Ω2. r dr

(6.28)

Using ∆u = δu+ξ·∇u = dξ/dt = −iωξ+(u·∇)ξ, we find that the Eulerian perturba-

tion δu is related to the Lagrangian displacement vector ξ by δu = −iω˜ ξ − rΩ′ξrφˆ

(prime donates radial derivative) and we can further combine Eqs. (6.21)-(6.27)

into two equations for ξr (radial Lagrangian displacement) and δΠ/ρ:

dξr dr

=

A11ξr

+

δΠ A12 ρ ,

d δΠ

δΠ

dr ρ = A21ξr + A22 ρ ,

where

A11

=

−

1 r

ω˜ 2

− 2mω˜ Ω + m2ω2Aφ ω˜ 2 − m2ω2Aφ

,

(6.29) (6.30)
(6.31)

118

m2 A12 = r2 ,

(6.32)

A21

=

ω˜ 2

−

m2ω2Aφ

−

2rΩ

dΩ dr

+

d ln 2

Bφ

−

1

d ln r

ω2Aφ

−

4

(ω˜ Ω + mω2Aφ)2 (ω˜ 2 − m2ω2Aφ)

,

(6.33)

A22

=

2m ω˜ Ω + mω2Aφ , r ω˜ 2 − m2ω2Aφ

(6.34)

and ωAφ ≡ vAφ/r = Bφ/(r 4πρ) is the toroidal Alfve´n frequency. Equations (6.29) and (6.30) are the same as Eqs. (119) and (120) (derived for a pure toroidal mag-

netic field) in Chandrasekhar (1961).

Note that although we start with a mixed magnetic field B = Bφφˆ + Bzˆz, the final Eqs. (6.29) and (6.30) do not contain Bz. The reason is that the z-component only appears in Eqs. (6.24) and (6.27), which in fact can be decoupled from the other five perturbation equations. Indeed, using Eq. (6.25) to replace δBr in Eq. (6.24), and combining with Eq. (6.27), we find

(ω˜ 2 − m2ω2Aφ)

δBz

+

i ω˜

dBz dr

δur

= 0.

(6.35)

In general, ω˜ 2 − m2ω2Aφ 0. Comparing the above equation with Eq. (6.27) we have δuz = 0. This is to be expected since the perturbed quantities are assumed to be independent of z. Also note that when the wave frequency ω is real, the coefficients A11, A21 and A22 are singular at

ω˜ 2 = m2ω2Aφ.

(6.36)

We shall call them the Magnetic Resonances (MRs). Obviously, they reduce to the corotation resonance when Bφ = 0.

119

6.4 Boundary Conditions

In general, the boundary of any magnetized flow should satisfy the following conditions:

ρun = 0,

[n · B] = 0,

P

+

ρu2n

+

B2t 8π

= 0,

ρunut

−

Bn Bt 4π

= 0,

(6.37) (6.38) (6.39)
(6.40)

where n is a unit vector normal to the boundary surface, the subscript n and t denote the normal and tangential components, and the square bracket represents the difference in a quantity across the boundary (e.g., Schmidt 1979; Shu 1992). For the system we study in this chapter, there is no radial background flow, we only need to consider Eqs. (6.38)-(6.40) with un = 0. Obviously, with the magnetic field continuous across the boundaries and with no radial field component, our equilibrium models constructed in Sec. 6.2 already satisfy the boundary conditions. In the perturbed state, the boundary conditions read

∆ [n · B] = 0, ∆ P + B2t = 0,
8π ∆ BnBt = 0.
4π

(6.41) (6.42) (6.43)

Note that ∆[n · B] = [(∆n) · B] + [n · ∆B], and since both n and ∆n are the same
across the boundary 2, Eq. (6.41) simply becomes n · [∆B] = 0, where we have
used [B] = 0 (as assumed in our model setup). Since (∆B)r = (δB + ξ · ∇B)r =
2Note that in the non-axisymmetric case, the perturbed surface normal vector ∆n is not the same as ˆr; see Schmidt (1979) for a derivation of ∆n.

120

δBr − ξφBφ/r, and both ξφ and Bφ are continuous across the boundaries, we find

[δBr] = 0.

(6.44)

On the other hand, Eqs. (6.41) and (6.42) combine to give

[∆Π] = 0.

(6.45)

The condition (6.43) is already satisfied because ∆ [BnBt] = [(∆Bn)Bt] + [Bn∆Bt] = 0 when the background field is continuous. We note that the perturbed magnetic field does not need to be continuous across the boundary. This means that there could be surface current induced by the perturbation.

To implement the two boundary conditions (6.44)-(6.45), we need to calculate the perturbed magnetic field in the vacuum region (r < r1 and r > r2). This can be done by solving

∇ × δB = 0, ∇ · δB = 0.

(6.46)

Clearly, δB is a potential field δB = ∇Ψ with Ψ (also ∝ eimφ−iωt) satisfying

∇2Ψ = 0.

(6.47)

The solution of Eq. (6.47) is

Ψ = C1rm + C2r−m,

(6.48)

where C1 and C2 are integration constants. Requiring δB to be regular at r → 0 and r → ∞, we find

δBr

=

C1mrm−1,

δBφ

=

C1

im r

rm,

for r < r1,

(6.49)

and

δBr

=

−C2mr−m−1,

δBφ

=

C2

im r

r−m

,

for r > r2.

(6.50)

121

The two constants C1 and C2 can be determined by using [δBr] = 0. At r = r1, this implies C1mrm−1 = imBφξr/r (see Eq. [6.25] for δBr inside the fluid zone). Thus,

C1 = iBφξrr−m|r1.

(6.51)

Similarly,

C2 = −iBφξrrm|r2.

(6.52)

Since the detailed realization of the boundary condition [∆Π] = 0 depends on the specific equilibrium model, we address the two models separately.

6.4.1 Model (a)

At the inner boundary r = r1, by using Eq. (6.51), we find that the Lagrangian

perturbation of the total pressure in the vacuum just inside r1 (i.e., r = r1−) is

given by

∆Π|r1−

=

BφδBφ 4π

−

ξr

B2φ 4πr

=

−(m

+

1)ξr

B2φ 4πr

.

(6.53)

In the fluid just outside r1 (i.e., r = r1+), we have

∆Π|r1+

=

δΠ

 + ξr r2p+1

−

v2Aφ r

−

1 r2

 1 − sv2Az0 − (1 + q)v2Aφ0  ρ.

(6.54)

Thus the condition [∆Π] = 0 at r = r1 can be written as

δΠ ρ

+

ξr

 r2p+1

+

mv2Aφ r

−

1 r2

 1 − sv2Az0 − (1 + q)v2Aφ0  = 0,

at r = r1+.

(6.55)

Note that in deriving the above equation, we have implicitly used [Bφ] = 0 and

[ξr] = 0. The same procedure yields the boundary condition at r = r2

δΠ ρ

+

ξr

 r2p+1

−

mv2Aφ r

−

1 r2

 1 − sv2Az0 − (1 + q)v2Aφ0  = 0,

at r = r2−.

(6.56)

122

Figure 6.3: The growth rate of Papaloizou-Pringle instability (in units of Ω0, the fluid rotation rate at the pressure maximum of the torus) as a function of the relative thickness of the torus for different values of m. The rotation profile is Ω ∝ r−2. This figure is similar to Fig. 1 in Blaes & Glatzel (1986) and to Fig. 1 in Abramowicz et al. (1987).

6.4.2 Model (b)

The derivation is similar to Model (a). In this case, the boundary conditions are

δΠ ρ + ξr

r2p+1

−

1 r2

1

−

v2Az0

−

1 1

+ −

r12 r12

v2Aφ0

= 0,

at r = r1+,

(6.57)

δΠ ρ

+

ξr

 r2p+1

−

mv2Aφ0 (1 − r12)2

1

−

r12 r2

1 − r2

1−

v2Az0

−

1+ 1−

r12 r12

v2Aφ0

 

=

0,

at r = r2−. (6.58)

123

Figure 6.4: The instability growth rate as a function of the relative thickness of the torus for a range of toroidal magnetic field strengths. Different lines represent different vAφ0 with the solid line denoting the hydrodynamic case. The upper and bottom panels are for the m = 1 and m = 2 modes, respectively.
6.5 Numerical Results
We employ the standard shooting method (Press et al. 1992) to solve Eqs. (6.29) and (6.30) subjected to the boundary conditions (6.55)-(6.56) [Model (a)] or (6.57)-(6.58) [Model (b)] to obtain the eigenvalue ω = ωr + iωi. For most of our analysis, we set the rotation index p = −2 (i.e., Ω ∝ r−2) such that our results can be directly compared with the original Papaloizou-Pringle instability.
124

Figure 6.5: The instability region in the parameter space defined by the dimensionless torus thickness ∆r/r2 and the toroidal magnetic field strength. The other parameters are fixed to m = 2, Ω ∝ r−2 and Bφ ∝ r. The dotted area denotes the region where a growing mode can be found.
Before discussing our results for finite magnetic fields, we briefly review the main features of the classical (hydrodynamical) PP instability. As seen in Fig. 6.3, the instability growth rate increases with increasing torus thickness for small ∆r/r2 but terminates at some finite thickness. As m increases, the termination point shifts to smaller ∆r/r2, although the peak growth rate remains approximately the same. This means that the PP instability only exists for relatively thin tori, as shown by Blaes & Glatzel (1986) and by Abramowicz et
125

Figure 6.6: The instability growth rate as a function of vAφ0 = vAφ(r0)/(r0Ω0) for tori with different thickness ∆r/r2. The upper and bottom panels depict the cases with m = 1 and m = 2, respectively.
al. (1987). The former also provides an approximate analytical expression for the limiting maximum growth rate as m → ∞. These features can be understood from the fact that the PP instability arises from the interaction of the surface gravity waves at the torus boundary. For a thin torus, the velocity shear across the corotation point is small, and there would not be enough shear rotational energy available to drive the growth. When the torus thickness is too large, the wave amplitudes at the corotation radius (where the waves exchange angular momentum) is too small to allow for adequate interactions. Thus, only for the
126

Figure 6.7: Some special radii for the m = 1 overstable mode in a torus as a function of the dimensionless magnetic field strength vAφ0. The upper and bottom panels correspond to a thin ∆r/r2 = 0.2) and thick (∆r/r2 = 0.7) torus, respectively. The two solid lines show the inner and outer torus boundaries (r1 and r2). The dotted line represents the corotation radius, while the shot-dashed and long-dashed lines show the inner and outer magnetic resonances (IMR and OMR), respectively [see Eq. (6.59)]. The thin solid line is the mode growth rate with the scale shown on the right. The torus rotation and magnetic field profiles are Ω ∝ r−2 and Bφ ∝ r.
127

Figure 6.8: The instability growth rate as a function of vAz0 for a torus with relative thickness ∆r/r2 = 0.2. Different lines represent different values of vAφ0. The other parameters are fixed to m = 2, Ω ∝ r−2, Bφ ∝ r and Bz ∝ r in the fluid.
“intermediate” torus thickness, with m∆r/r ∼< 1, will the instability operates. This explains why for a larger m, the PP instability terminates at a smaller torus thickness.
128

Figure 6.9: The instability growth rate as a function of vAφ0 for accretion tori described by Model (b) with a pure toroidal magnetic field. The solid and dotted lines correspond to two different inner disc boundary radii. The bottom pane shows the corresponding torus relative thickness.
6.5.1 Model (a): Pure Toroidal Field Configuration
In this section, we present the numerical results for an accretion torus with a pure power-law profile toroidal magnetic field. We choose the power-law index q = 1 so that the vertical current density is uniform as in Model (b).
Figure 6.4 shows the growth rate ωi as a function of thickness ∆r/r2 for different vAφ0. The two panels share similar characteristics: (i) For relatively weak
129

Figure 6.10: The instability growth rate as a function of vAz0 for a thin torus (with r1 = 0.9 and ∆r/r2 ≃ 0.2). The different lines represent different values of vAφ0. Unlike Fig. 6.8, here the magnetic field profile is described by Model (b).
B field (vAφ0 0.1), the instability resembles the B = 0 limit in that it always starts from infinitely small thickness and terminates beyond a certain ∆r/r2; (ii) For stronger B fields, the instability starts beyond certain finite ∆r/r2 and then extends all the way to very large thickness (although as ∆r/r2 approaches unity, the growth rate becomes increasingly small). As vAφ0 increases, the critical thickness for the onset of instability also increases.
Figure 6.5 maps out the unstable zone in the thickness – magnetic field
130

strength parameter space. It shows the similar feature as as Fig. 6.4. We can see that the unstable region is mainly located at the lower-left (thin torus with weak B field) and the upper-right (thicker torus with strong B field) corners of the parameter space.

In Fig. 6.6, we present our numerical results in a different way. We fix the dimensionless thickness ∆r/r2 and plot the growth rate as a function of magnetic field strength. For a thin torus, we find that as vAφ0 increases, the growth rate first goes up slightly compared to the B = 0 case, then decreases and becomes completely suppressed when the magnetic field is sufficiently strong (vAφ comparable to rotation velocity). For a thick torus, the instability can survive for a relatively stronger B field and then vanishes beyond a certain vAφ0.

To probe the underlying physics of how magnetic fields affect the PP instability, we show in Fig. 6.7 the locations of several special points in the fluid: The corotation radius rc is where the wave pattern corotates with with the background flow, i.e., ω˜ r = ωr − mΩ = 0. The inner/outer magnetic resonances (IMR/OMR) are defined by [see Eq. (6.36)]

ω˜ r = ωr − mΩ = ±mωAφ.

(6.59)

At the IMR, the wave is trailing the background flow but corotates with the azimuthal Alfve´n wave traveling in the counter-rotational direction (viewed in the corotating frame), while at the OMR, the wave is leading the background flow and corotates with the Alfve´n wave in the rotational direction. Recall that for PP instability to operate in the B = 0 limit, it is essential that the corotation radius lies in between torus boundaries (i.e., r1 < rc < r2). Now, with the inclusion of the magnetic field, we see from Fig. 6.7 that as vAφ0 increases, both rc and rOMR shift beyond the outer boundary of the torus. The IMR radius, rIMR, however,

131

always stays inside the fluid. This suggests that in a magnetic torus, the IMR plays a similar role as the corotation resonance does in a non-magnetic torus.
6.5.2 Other Magnetic Field Configurations
Model (a): Mixed magnetic field
Although the vertical magnetic field Bz does not enter into the perturbation equations, the presence of a finite Bz can affect the mode growth rate through the boundary conditions. Figure 6.8 shows some results for the accretion tori with a mixture of vertical and toroidal magnetic fields. We take the power-law index of Bz to be s = 1 so that the azimuthal electric current density is constant. In Fig. 6.8, we fix the toroidal field and plot the mode growth rate as a function of the vertical field strength. We find that the effect of finite vAz0 is small (note the scale of the y-axis). For vAz = 0, the results agree with what is shown in Fig. 6.6b.
Model (b): Pure toroidal field
In this case, since the background toroidal magnetic field has a profile that depends on the inner boundary radius r1, to solve for the equilibrium structure and the global mode, we must first specify r1. Once we fix r1, we can easily solve for the other boundary radius r2. In Fig. 6.9, we show in the upper panel the mode growth rate as a function of vAφ0. The result is qualitatively similar to Fig. 6.6b. The bottom panel shows that for a fixed inner boundary radius r1, the thickness does not change appreciably as vAφ0 varies. So the two values of r1 we choose adequately depict the thin and thick tori, respectively. Again, we see
132

that for a thin torus, the original Papaloizou-Pringle instability is suppressed by the toroidal field, while for a thick torus, the instability can survive for larger field strengths.
Model (b): Mixed magnetic field
In Fig. 6.10, we show the mode growth rate as a function of vAz0 for different fixed values of vAφ0. Similar to the case shown in Fig. 6.8, we see that the vertical magnetic field has a small effect on the stability property of a magnetized accretion torus. Again, this is understandable given that Bz does not enter into the differential equations for the perturbations, but only affects the modes through boundary conditions.
6.6 Summary
In this chapter, we have studied the effect of magnetic fields on the global nonaxisymmetric instability (the PP instability) in accretion tori. For simplicity, we assume that both the perturbation and the background flow variables have no z-dependance (thus our tori are essentially 2D cylinders). We have explored various possible magnetic configurations in the torus. Although the detailed property of the instability is model-dependent, Figs. 6.4-6.6 illustrate our general findings: (i) For thin tori (with the dimensionless thickness ∆r/r2 ∼< 0.2, where r2 is the outer torus radius), the instability exists for zero and weak magnetic fields, but is suppressed when the toroidal field becomes sufficiently strong (with the corresponding Alfve´n speed vAφ ∼> 0.2rΩ measured at gas pressure maximum); (ii) For thicker tori (∆r/r2 ∼> 0.4), the PP instability does not operate for zero and
133

weak magnetic fields, but becomes active when the field is sufficiently strong (vAφ ∼> 0.2rΩ measured at gas pressure maximum). A vertical magnetic field also influences the PP instability, but its effect is generally smaller that that of the toroidal field.
It is difficult to precisely pin down the physical origin of the magnetic field effect on the PP instability. For example, with a finite toroidal magnetic field, we find that the corotation resonance radius may lie outside the torus body, and yet the torus is still unstable. On the other hand, the inner magnetic resonance radius, where ωr − mΩ = −mωAφ [see Eq. (6.59)], always lie inside the fluid body. Thus we suspect that in a magnetic torus, the inner magnetic resonance plays a similar role as the corotation resonance does in a non-magnetic torus.
Note that boundary conditions play an important role in the PP instability. In the hydrodynamical case, if we had used a rigid boundary condition (ξr = 0) (as adopted in the standard treatment of laboratory Couette flows; e.g., Chandrasekhar 1961), the torus with constant angular momentum profile would be stable to two-dimensional perturbations according to the Rayleigh’s inflexionpoint theorem (see Drazin & Reid 1981). Including a finite magnetic field, the theorem cannot not be applied (our attempt to generalize the theorem did not reveal a simple instability criterion). The instability for thick, magnetized tori found in this chapter may be a result of magnetic shear instability.
We note that the PP instability (or its magnetic generalization) involves wave modes that do not have vertical structure (i.e. kz = 0). Thus it is distinctly different from the usual MRI (Balbus & Hawley 1998). Our finding about the instability of thick tori with strong magnetic fields is particularly interesting: Since the MRI can be suppressed when the magnetic field is too strong, our
134

results suggest that magnetized tori may be subject to the instability even when it is stable against the usual MRI.
135

CHAPTER 7 LOW-T/|W| INSTABILITIES IN DIFFERENTIALLY ROTATING
PROTO-NEUTRON STARS WITH MAGNETIC FIELDS
7.1 Introduction
Recent hydrodynamical simulations have shown that differentially rotating neutron stars formed in core-collapse supernovae may develop global nonaxisymmetric instabilities even when T/|W| (the ratio of the rotational kinetic energy T to the gravitational potential energy |W|) is relatively small (less than 0.1). Such low-T/|W| instability can give rise to efficient gravitational wave emission from the proto-neutron star. In this chapter, we investigate how this instability is affected by magnetic fields using a cylindrical stellar model. Wave absorption at the corotation resonance plays an important role in facilitating the hydrodynamic low-T/|W| instability. In the presence of a toroidal magnetic field, the corotation resonance is split into two magnetic resonances where wave absorptions take place. We show that the toroidal magnetic field suppresses the lowT/|W| instability when the total magnetic energy WB is of order 0.2 T or larger, corresponding to toroidal fields of a few ×1016 G or stronger. Although poloidal magnetic fields do not influence the instability directly, they can affect the instability by generating toroidal fields through linear winding of the initial poloidal field and magneto-rotational instability. We show that an initial poloidal field with strength as small as 1014 G may suppress the low-T/|W| instability.
136

7.2 Equilibrium Model of A Magnetized Rotating Cylinder

We consider a rotating star with purely toroidal magnetic fields and assume a polytropic equation of state

P = KρΓ = Kρ1+1/N ,

(7.1)

where P and ρ are the gas pressure and density, K, Γ and N are constants. Al-

though hydromagnetic stellar equilibrium models can be constructed easily us-

ing the HSCF method (Hachisu 1986; Tomimura & Eriguchi 2005; see Lander &

Jones 2009 for recent works on uniformly rotating stars), linear eigenvalue anal-

ysis for such models is difficult. Thus we follow the setup in Saijo & Yoshida

(2006) by treating the star as an infinite cylinder. We adopt the cylindrical co-

ordinates (r, φ, z). All the background variables are assumed to be functions of

cylindrical radius r only. The equilibrium state of the cylinder is determined by

force balance equation in the radial direction

1 dP = − dΦ + rΩ2 − 1 dPm − B2φ ,

ρ dr dr

ρ dr 4πρr

(7.2)

where Ω is the flow rotation rate, Bφ is the toroidal magnetic field strength, Pm = B2φ/8π is the magnetic pressure, and Φ is the Newtonian gravitational potential which relates to density ρ via Poisson’s equation

∇2Φ = 4πGρ.

(7.3)

Eliminating Φ from Eqs. (7.2) and (7.3) yields

d dr

r dP ρ dr

= −4πGρr + d (r2Ω2) − d dr dr

r dPm ρ dr

−

d dr

 

B2φ 4πρ

 

.

(7.4)

For numerical convenience, we nondimensionalize variables as follows

ρ = ρˆρc = θNρc,

(7.5)

137

r = rˆ (N + 1)Kρ1c/N−1 , 4πG
Ω = Ωˆ 4πGρc,
Bφ = Bˆφ 4π(N + 1)Kρ1c+1/N ,

(7.6) (7.7) (7.8)

where ρc is the central density and the hatted variables denote dimensionless

quantities. Similar dimensionless variables for other quantities can be con-

structed from the list above. We follow Saijo & Yoshida (2006) to adopt the

following rotation profile:

Ωˆ

=

C rˆ2 +

, A

(7.9)

where A and C are constants. The toroidal B field profile we employ is

Bˆφ = brˆ(Rˆ − rˆ)

(7.10)

where Rˆ denotes the dimensionless boundary of the cylinder and the constant b specifies the field strength. For simplicity, we will omit the hats on all variables hereafter unless otherwise noted. The above profile implies that Bφ = 0 at both the center and the surface of the star. For small r, we have Bφ ≃ brR, implying a constant axial current for r → 0. Eq. (7.4) in dimensionless form now reads

d2θ dr2 +

1 − Nθ−N−1b2r(R − r)(R − 2r) r

dθ + θN dr

+

4b2(R − r)(R − 2r) − b2r(3R − 4r)

θ−N

d = 2Ω (rΩ).

dr

(7.11)

For an nonmagnetized star, Eq. (7.11) reduces to Eq. (3.2) in Saijo & Yoshida

(2006). In the limit of zero rotation, it recovers the well-known Lane-Emden

equation in cylindrical geometry.

In the limit of Bφ = 0, we can simply integrate Eq. (7.11) starting from r = 0 using boundary condition that θ = 1 and θ′ = 0 to a point where θ goes to zero,

138

Figure 7.1: Density profiles for hydrodynamic and hydromagnetic equilibria. The polytropic index is N = 1. The solid, short-dashed and long-dashed lines represent different rotation profiles for non-magnetic models, while the dotted line is for a magnetic model with WB/|W| = 0.03.

which defines the cylinder surface R. The hydrodynamic equilibrium can thus

be easily constructed. When Bφ is non-zero, we choose an initial guess for the

surface radius R based on results for the equivalent hydrodynamic model and

integrate Eq. (7.11) imposing the same boundary condition at the center. We

stop the integration at r = R to check the value of θ. We then adjust our guess

for R and go through the same process, until θ|r=R comes close enough to 0. For a

given equilibrium state, the rotational kinetic energy T , gravitational potential

energy W and magnetic energy WB of the cylinder have the following form

1 T=2 W=−

R
ρr2Ω2dV = θNΩ2r3dr,
0

ρr dΦdV =

R2
θNrdr ,

dr 0

139

WB =

B2φ dV

=

R6b2 ,

8π 60

(7.12)

where all the variables are dimensionless and the corresponding physical unit

for energy is (N + 1)2K2ρ2c/N/4. Examples of the equilibrium density profile are

given in Fig. 7.1. We see that the density profile is not always monotonic: for

large C and small A (i.e, large rotation rate and large degree of differential rota-

tion), the density maximum is off-centered.

7.3 Linear Perturbation Analysis

7.3.1 Perturbation equations

The cylindrical flow we are considering satisfies the usual ideal MHD equations

∂ρ + ∇ · (ρv) = 0,
∂t ∂v 1 1 ∂t + (v · ∇)v = −ρ ∇Π − ∇Φ + 4πρ (B · ∇)B,
∂B ∂t = ∇ × (v × B),
∇ · B = 0,

(7.13) (7.14) (7.15) (7.16)

∇2Φ = 4πGρ,

(7.17)

where Π = P + Pm is the total pressure. We apply linear perturbations to the above equations by assuming the perturbation of any physical variable f to

have the form δ f ∝ eimφ−iωt with m being the azimuthal mode number and ω

the wave frequency. The resulting linearized perturbation equations contain

variables δv, δρ, δΠ, δΦ and δB. To simply the algebra, we define a new variable

δh

=

δΠ ρ

=

δP ρ

+

B · δB. 4πρ

140

Using ∆v = δv + ξ · ∇v = dξ/dt = −iωξ + (v · ∇)ξ = −iωξ + Ω∂ξ/∂φ, we find that the Eulerian perturbation δv is related to the Lagrangian displacement vector ξ by δv = −iω˜ ξ − rΩ′ξrφˆ (prime denotes radial derivative). In terms of ξr, δh and δΦ, the MHD perturbation equations (in dimensionless form) can be cast into four first-order differential equations:

dξr dr

=

A11ξr

+ A12δh +

dδΦ A13δΦ + A14 dr ,

(7.18)

dδh dδΦ dr = A21ξr + A22δh + A23δΦ + A24 dr ,

dδΦ dδΦ dr = A31ξr + A32δh + A33δΦ + A34 dr ,

d dδΦ

dδΦ

dr dr = A41ξr + A42δh + A43δΦ + A44 dr ,

(7.19) (7.20) (7.21)

where

A11

=

rω˜ 2 (ω2Aφ − Ω2)ω˜ 2 + ω2Aφω2 (c2s + v2Aφ)(ω˜ 2 − m2ω2Aφ)(ω˜ 2 − ω2s)

+ (c2s

+

gω˜ 2 v2Aφ)(ω˜ 2

−

ω2s )

−

ω˜ 2

+ 2mω˜ Ω + m2ω2Aφ r(ω˜ 2 − m2ω2Aφ)

,

ω˜ 4 m2

A12

=

− (c2s

+

v2Aφ)(ω˜ 2

−

m2ω2Aφ)(ω˜ 2

−

ω2s )

+

r2(ω˜ 2

−

, m2ω2Aφ)

m2

A13

=

, r2ω˜ 2

(7.22) (7.23) (7.24)

A14 = 0,

(7.25)

A21

=

ω˜ 2

−

m2ω2Aφ

−

4(mω2Aφ + ω˜ Ω)2 ω˜ 2 − m2ω2Aφ

+

r

d dr

(ω2Aφ

−

Ω2)

+

(ω2Aφ

−

Ω2) r ρ

dρ dr

+

g ρ

dρ dr

1 + (c2s + v2Aφ)(ω˜ 2 − m2ω2Aφ)(ω˜ 2 − ω2s)

r

(ω2Aφ − Ω2)ω˜ 2 + ω2Aφω2

+ g(ω˜ 2 − m2ω2Aφ) 2 ,

(7.26)

A22

=

−

rω˜ 2 (ω2Aφ − Ω2)ω˜ 2 + ω2Aφω2 (c2s + v2Aφ)(ω˜ 2 − m2ω2Aφ)(ω˜ 2 − ω2s)

− (c2s

+

gω˜ 2 v2Aφ)(ω˜ 2

−

ω2s )

+

2m(mω2Aφ + ω˜ Ω) r(ω˜ 2 − m2ω2Aφ)

−

1 ρ

dρ .
dr

(7.27)

141

2mΩ

A23 =

, rω˜

A24 = −1,

(7.28) (7.29)

A31 = A32 = A33 = 0,

(7.30)

A34 = 1

A41

=

ρ −
r

−

dρ dr

+

ρ

 

m2ω2Aφ ω˜ 2

−

 1

A11

−

ρ

m2ω2Aφ rω˜ 2

−

2ρmΩ ,
rω˜

(7.31) (7.32)

A42

=

ρ

 

m2ω2Aφ ω˜ 2

−

 1 A12

+

m2 ρ,
r2ω˜ 2

(7.33)

A43

=

ρ m4ω2Aφ r2ω˜ 4

+

m2 ,
r2

1 A44 = − r .

(7.34) (7.35)

In the above expressions, ω˜ = ω − mΩ is the wave frequency in the co-rotating

frame, ρ = θN is the dimensionless density, cs =

√ dP/dρ = θ/N is the dimen-

sionless sound speed,

vAφ =

B2φ = br(R − r)θ−N/2 4πρ

(7.36)

is the toroidal Alfve´n velocity, ωAφ = vAφ/r = b(R − r)θ−N/2 is the toroidal

Alfve´n frequency,

ωs =

c2s

c2s + v2Aφ

m2ω2Aφ

is the slow magnetosonic wave frequency for k = (m/r)φˆ, and

(7.37)

g = dΦ = rΩ2 − dθ − b2r(R − r)(2R − 3r) θ−N dr dr

(7.38)

is the gravitational acceleration in radial direction.

142

7.3.2 Boundary conditions

To solve Eqs. (7.18)-(7.21) as an eigenvalue problem, we need four boundary

conditions. The outer boundary conditions are straightforward. From the per-

turbed Poisson equation, we know the perturbed potential outside the star scales as δΦ ∝ r−m. By requiring this potential to match smoothly with the po-

tential inside, we obtain our first outer boundary condition:

dδΦ dr

+

m r

δΦ

=

0

at

r = R.

(7.39)

Requiring the Lagrangian pressure perturbation to vanish at the stellar surface

yields

dθ δh + dr ξr = 0 at r = R.

(7.40)

The inner boundary conditions are more involved. As r → 0, we observe that

g → 0, Ω → constant, ρ → constant, ωAφ ∝ Bφ/r → constant and δρ is finite.

Thus, near the center of the star, Eqs. (7.18)-(7.21) can be simplified as

dξr dr

=

X −

+ mY X

ξr r

+

m2 δh X r2

+

m2 δΦ ω˜ 2 r2

dδh X2 − Y2 mY δh 2mΩ δΦ dδΦ

= dr

X

ξr + X

+ r

ω˜

−, r dr

1d r dr

r dδΦ dr

m2 − r2 δΦ = 0,

(7.41) (7.42) (7.43)

where,

X = ω˜ 2 − m2ω2Aφ, Y = 2(mω2Aφ + ω˜ Ω).

(7.44)

Eq. (7.43) has the following solution

δΦ = rm, at r ∼ 0,

(7.45)

where we have taken the coefficient of the power-law term to be 1 for simplicity. The perturbations ξr and δh generally take the form

ξr = C1rm−1 + C2rm−1 ln r, δh = C3rm + C4rm ln r,

(7.46)

143

where C1, C2, C3, C4 are constants so that perturbations remain regular at the center. Eq. (B.1) represents the leading terms of the Frobenius expansions of these functions.

Substituting solution (B.1) into Eqs. (7.41) and (7.42) leads to two equations

which have the structure a1 + a2 ln r = 0, where a1, a2 are constants that depend

on the values of m, Ω, ω˜ , ωAφ, C1, C2, C3, C4 near the center. Since these equa-

tions should be satisfied everywhere around the center, we demand a1 = a2 = 0.

This yields

C2

+

m(X + X

Y) C1

=

m2 X C3

+

m2 ,
ω˜ 2

(7.47)

mC4 = (X + Y)C2,

m(X − Y)

X2 − Y2

2Ω − ω˜

C4 + X C3 = X C1 + m ω˜ ,

(7.48) (7.49)

from which we can determine three constants

C2

=

m3(m + −
2X

2) ω2Aφ,

(7.50)

C4

=

m2(m + −
2

2)

X

+ X

Y

ω2Aφ,

C1

=

X

m + Y C3

+

m2(m + 2) 2(X + Y)

ω2Aφ

+

m ω˜ 2

X

X +

, Y

(7.51) (7.52)

once we specify C3. When we solve the eigenvalue problem, C3 will be deter-

mined together with the eigenfrequency ω. We see that since for the specific Bφ

profile we are considering (see Eq. [7.10]), ωAφ remains approximately constant

near the center, C2 and C4 are both finite. Therefore the logarithmic parts in

solution (B.1) cannot be neglected.

In the hydrodynamic limit (Bφ = 0, ωAφ = 0), we have C2 = C4 = 0 so that the solutions of ξr and δh take a purely power-law form, and Eq. (7.52) reduces to

m C1 = ω˜ (ω˜ + 2Ω) (C3 + 1),

(7.53)

144

which is equivalent to

m

ξr

=

rω˜ (ω˜

+

(δh 2Ω)

+

δΦ)

at

r ∼ 0.

(7.54)

Again, constant C3 = δh/δΦ will be determined as a part of the eigenvalue problem. Clearly, for the magnetic cases where d ln Bφ/d ln r > 1 at r ∼ 0 so that ωAφ → 0 as r → 0, the above hydrodynamic boundary condition is also valid.

7.3.3 Cowling approximation

In the Cowling approximation, we neglect the gravitational potential perturbation δΦ. The perturbation equations then become

dξr dr

=

A11ξr

+

A12δh,

dδh dr = A21ξr + A22δh,

(7.55) (7.56)

with the four coefficients given by the same equations as before. Similarly, for

r → 0 the simplified version of Eqs. (7.41) and (7.42) are

dξr dr

=

X −

+ mY X

ξr r

+

m2 X

δh r2

dδh X2 − Y2 mY δh

dr =

X

ξr +

X

. r

The outer boundary condition in this case is again given by Eq. (7.40)

(7.57) (7.58)

dθ δh + dr ξr = 0

at

r = R.

(7.59)

The inner boundary condition can be obtained by substituting the power-law solutions ξr ∝ rm−1 and δh ∝ rm into Eqs. (7.57) and (7.58), giving

m ξr = r[ω˜ 2 + 2ω˜ Ω − m(m − 2)ω2Aφ] δh.

(7.60)

145

Figure 7.2: The m = 2 mode frequency as a function of T/|W| with and with out Cowling approximation. The upper and bottom panels show the real and imaginary parts of the frequency, respectively, with Ωc being the rotation frequency at the center. The star has no magnetic field and the polytropic index is N = 1.

Note that for either m = 2 perturbations or a unmagnetized flow, the above inner B.C. reduces to the same form

m

ξr

=

rω˜ (ω˜ +

δh. 2Ω)

(7.61)

7.4 Numerical Results
For most part of this section, we will employ the standard shooting method (Press et al 1992) to solve the two ODEs, Eqs. (7.55) and (7.56), subject to boundary conditions (7.59) and (7.60). We focus on the effects of toroidal magnetic fields on the low-T/|W| instability previously found for purely hydrodynamic stars.
146

Figure 7.3: The m = 2 mode frequency as a function of WB/|W| for stellar models with different rotation profiles (thus different T/|W|’s). The upper and bottom panels show the real and imaginary parts of the frequency, respectively, with Ωc being the rotation frequency at the center. The other parameters are the same as in Fig. 7.2.
Before moving on to our main results, let us first examine the validity of Cowling approximation for low-T/|W| instability. To this end, we compare the eigenfrequency calculation with and without Cowling approximation. In Fig. 7.2, we fix one of the rotation parameters, A, while change the other, C, to obtain different values of T/|W|. We see that for the whole range of T/|W| considered, the real part of the mode frequency does not show much difference between those two cases. The bottom panel of Fig. 7.2 shows that, for the low T/|W| range ( 0.1), the mode growth rate exhibits qualitatively similar behavior
147

Figure 7.4: The m = 2 mode growth rate as a function of T/|W| with and without toroidal magnetic field. The other parameters are A = 0.86 and N = 1.
with an approximate factor of 2 difference between the two cases. For relatively large T/|W| ( 0.2), the growth rate no longer follows the similar trend when T/|W| increases. Shibata et al. (2002) also found from their hydrodynamic simulation of a similar stellar model that the mode growth rate declines beyond certain T/|W|. The solid lines (“no Cowling”) in our Fig. 7.2 agree well with the results depicted in Fig. 4 of Shibata et al.. With Cowling approximation, we find that the growth always increases with increasing T/|W|. Overall, Fig. 7.2 shows that using the Cowling approximation captures the essential feature of the lowT/|W| instability, especially when T/|W| is not much larger than the threshold.
Figs. 7.3 and 7.4 contain the most important results of this chapter. In Fig. 7.3,
148

Figure 7.5: Example wavefunctions of unstable low-T/|W| (≃ 0.095) mode with A = 0.86, C = 0.9, m = 2 and N = 1. The left column shows the radial displacement as a function of radius whereas the right column shows the radial velocity perturbation, with the solid and short-dashed lines representing the real and imaginary parts, respectively. The upper and lower panels are for nonmagnetic and magnetic stellar models, respectively. The dotted lines indicate the location of the corotation resonance (in the nonmagnetic case) or slow magnetosonic resonances (in the magnetic case). The vertical scales of the wavefunctions are arbitrary.
we plot the eigenfrequency of the m = 2 mode as a function of WB/|W| (the ra-
tio of magnetic energy to gravitational energy) for different rotation profiles.
Note that as we change the magnetic field strength, the equilibrium structure,
therefore T/|W| will also change. However, for the range of WB/|W| we consid-
ered, the modification to the equilibrium structure is so small (see dotted line in
Fig. 7.1 for the small modification) that T/|W| is approximately a constant along
the three curves. Fig. 7.3 demonstrates that the low-T/|W| instability can be sup-
pressed by the toroidal magnetic field. The point where the mode growth is
completely suppressed corresponds to WB/|W| ∼ 0.2 T/|W|. In this figure, we
choose those three models so that one can readily see that the instability is more
149

Figure 7.6: Angular momentum carried by the wave as a function of r. The model parameters are the same as in Fig. 7.5. The upper and lower panels show the nonmagnetic and magnetic models, respectively. The locations of the corotation resonance and slow resonances are indicated by the vertical dotted lines.
prominent with larger degree of differential rotation (smaller A) and larger rotation rate (larger C). Fig. 7.4 shows the mode growth rate as a function of T/|W| for stellar models with different WB/|W|. We see that the finite magnetic field shifts the curve towards larger T/|W|. In particular, the magnetic field increases the threshold for the instability from T/|W| ≃ 0.03 for the nonmagnetic model to T/|W| ≃ 0.035 for the WB/|W| ≃ 0.001 model. This finding can be easily understood: increasing rotation drives the instability, whereas magnetic field suppresses the instability. Therefore when a finite B field is included, in order to maintain the instability a larger rotation rate is needed to overcome the sup-
150

pressing effect.
Fig. 7.5 depicts two example wavefunctions of the overstable low-T/|W| mode. In the nonmagnetic case, the perturbation equations are singular at the corotation radius where ω˜ = ω − mΩ = 0. For low-T/|W| modes, the corotation resonance lies inside the star, so both the radial displacement and the gradient of the radial velocity perturbation undergo large variations across the corotation resonance (see the upper panels). In the magnetic case, however, the corotation resonance is no longer a singularity. Instead, the perturbation equations are singular at two slow magnetosonic resonances (Fu & Lai 2011a) where ω˜ = ±ωs, with ωs given by Eq. (7.37) 1. Therefore, the wavefunctions exhibit sudden changes at these two particular locations (see the lower panels). This splitting of corotation resonance into two magnetic slow resonances can also be seen in the angular momentum flux. In Fig. 7.6, we show the angular momentum carried by the wave across the star as a function of radius (see Fu & Lai 2011a for the flux formula). In the upper panel (nonmagnetic case), we see that F(r) experiences a sudden jump at the corotation resonance, whereas in the lower panel (magnetic case), two jumps occur at the two slow resonances and have different signs. This is similar to thin accretion discs studied in Fu & Lai (2011a). However, since the outer boundary condition we employed here (free surface) is totally different from the one used in Fu & Lai (2011a) (outgoing waves), we cannot directly relate the flux jump (or the net jump in the case of two resonances) to the magnitude of the growth rate. In any case, it is clear from
1Note that ω˜ 2 − m2ω2Aφ = 0 is not a singularity even though it appears to be a singular term similar to ω˜ 2 − ω2s in those coefficients of differential equations. This apparent singular term, one can show, will be canceled by some subtle mathematical manipulations. But this cancelation only works for the particular setup we consider here (pure toroidal B field, no vertical structure in perturbations). In general (i.e., with mixed B field or finite kz), equations will be singular at both ω˜ 2 = m2ω2Aφ and ω˜ 2 = ω2s (see Fu & Lai 2011a).
151

Fig. 7.6 that the corotation resonance indeed plays an important role in driving the hydrodynamic low-T/|W| instability and the toroidal magnetic field affects the instability by splitting the corotation resonance into two magnetic slow resonances. The property of the unstable mode in the presence of a magnetic field is determined by the combined effects from both slow resonances.
7.5 Summary
Recent studies of rotating (but nonmagnetic) core-collapse supernovae (e.g., Dimmelmeier et al. 2008) have demonstrated that newly formed neutron stars can develop nonaxisymmetric global instabilities with low T/|W|, and such instabilities lead to significant gravitational wave emission. In this chapter, we have carried out the linear stability analysis of magnetic, differentially rotating stars (modeled as a cylinder) to examine how magnegic fields affect the low-T/|W| rotational instability. We show that the wave absorption at the corotation resonance plays an important role in the instability. In the presence of a toroidal magnetic field, the corotation resonance is split into two magnetic resonances, where wave absorptions of opposite signs take place. Our main result is that toroidal magnetic fields reduce the growth rate of the low-T/|W| instability and increase the threshold T/|W| value above which the instability occurs. To significantly affect the instability, the required WB/|W| (the ratio of the magnetic energy WB to the gravitational potential energy |W|) should be of order 0.2 T/|W| or larger (see Figs. 7.3-7.4). As the critical T/|W| ranges from 0.01 to 0.1, the required WB/|W| lies between 0.002 and 0.02. Using |W| ∼ (3/5)GM2/R and
152

WB ∼ (B2φ/8π)(4πR3/3), we have

WB ∼ 1 |W| 300

Bφ 2 2 × 1016 G

R4 20 km

M −2 .
1.4M⊙

(7.62)

Thus, only toroidal magnetic fields stronger than 2 × 1016 G can significantly

affect the low-T/|W| instability.

In our simple (cylindrical) stellar model, poloidal (vertical) magnetic fields do not directly affect the rotational instability because the unstable modes do not have vertical structure (i.e., the vertical wavenumber is zero). We believe that this also holds for more realistic stellar models (spherical geometry, nonzero vertical wavenumbers, etc.), although more investigations are needed.

Nevertheless, even a relatively weak poloidal magnetic field present in the proto-neutron star may indirectly affect the T/|W| instability. In the core-collapse supernova scenario, differential rotation naturally arises inside the stellar core during the collapse (e.g., Akiyama & Wheeler 2005; Ott et al. 2006). This differential rotation can generate significant toroidal magnetic fields by winding the initial poloidal field and by magneto-rotational instability (Akiyama et al. 2003; Obergaulinger et al. 2009). Consider first the linear winding of the poloidal field Bp. The toroidal field grows in time as Bφ ∼ Bp∆Ωt, where ∆Ω is the difference in the rotation rate across the proto-neutron star. Thus the ratio of magnetic energy WB ∼ B2φR3/6 and the rotational energy T ∼ 0.2MR2(∆Ω)2 increases as WB/T ∼ B2pRt2/M. The time to reach a given WB/T ≡ f is then ttwist ∼ ( f M/B2pR)1/2 = f /3 R/vAp, where vAp = Bp/ 4πρ is the Alfve´n speed associated with Bp. On the other hand, the growth time of the low-T/|W| instability is tgrow ∼ 1/ωi ∼ (R3/GM)1/2/ωˆ i, where ωˆ i is the dimensionless growth rate in units of the Keplerian frequency (ωˆ i is approximately the vertical axis of the bottom panel in Fig. 7.3). For the low-T/|W| instability to operate before being

153

suppressed by the large Bφ (generated by twisting the initial poloidal field Bp),

we require ttwist tgrow, i.e.,

Bp

ωˆ i

GM2 1/2

f R4

≃

8 × 1013 G

f 1/2 ωˆ i 0.2 10−3

M 1.4M⊙

R −2 .
20 km

(7.63)

Since ωˆ i is of order 10−3 or larger, and the toroidal field suppresses the instability

when f = WB/T ∼ 0.2 (see Fig. 7.3), we see that an initial poloidal field strong than 1014 G can lead to the suppression of the instability. In other words, when

the initial poloidal field is less than 1014 G, the toroidal field will not grow fast

enough by linear winding so that the low-T/|W| instability still has a chance to

develop.

The effect of magneto-rotational instability (MRI) is harder to quantify. In the linear regime, MRI operates in modes with vertical structure (i.e., finite vertical wave number), which are independent from the T/|W|-unstable modes studied in this chapter. However, the nonlinear development of MRI may generate significant magnetic fields (both poloidal and toroidal) on a short timescale (of order the rotation period). There have been many MHD simulations of core collapse supernovae (e.g., Ardeljan et al. 2000, 2005; Kotake et al. 2004; Yamada & Sawai 2004; Obergaulinger et al. 2006; Burrows et al. 2007). Most of these simulations cannot resolve the MRI unless they employ drastically strong initial fields. It has been suggested that when MRI saturates, a large fraction of the kinetic energy in the differential rotation is converted to the magnetic energy (Akiyama et al. 2003; Obergaulinger et al. 2009), i.e., WB/T may approach unity on a dynamical time. Our result in this chapter shows that the T/|W| instability is strongly reduced when WB/|T | reaches 0.2. Therefore it would be important to quantify the saturation field of the MRI in proto-neutron stars. In addition, the MRI can lead to efficient angular momentum transport in different region of the star. This may also affect the T/|W| instability. Clearly, these issues must be

154

resolved in order to evaluate whether the low-T/|W| instability can develop in astrophysically realistic proto-neutron stars.
155

CHAPTER 8 CONCLUSION
Quasi-periodic variability has been observed from a number of Galactic compact X-ray binary systems. Of particular interest are several accreting black hole (BH) binaries which show pairs of quasi-periodic oscillations (QPOs) of fixed frequencies having ratios close to 2 : 3. Various theories/models have been proposed to explain the origin of QPOs. Perhaps the theoretically most developed model for QPOs is the relativistic diskoseismic oscillation model, in which general relativistic (GR) effect produces trapped modes/oscillations in the inner region of the disk. C-modes (corrugation modes) are non-axisymmetric and have nodes in vertical direction. They generally have frequencies lower than the observed QPO frequencies and they suffer corotational damping (Tsang & Lai 2009a). G-modes (inertial-gravity modes) used to be the most promising and most studied discoseismic modes as their trapping zone naturally arise due to GR potential and has nothing to do with the disk boundaries. However, recent progress in this field showed that the survival of g-modes faces two main challenges: firstly, the non-axisymmetric g-modes, like c-modes, are also subject to corotation damping (Kato 2003; Li, Goodman & Narayan 2003; Silbergleit & Wagoner 2007); secondly, its trapping zone can be easily destroyed by a weak (subthermal) magnetic field as we demonstrated in Chapter 2.
P-modes (inertial-acoustic modes) have no vertical structure, are trapped between the inner disc edge and the inner Lindblad resonance, and have the character of inertial-fast magnetosonic waves in their propagation zone. Previous works have shown that in a hydrodynamic disc (with B = 0) these modes can be overstable due to wave absorption at the corotation resonance (CR) when
156

the vortensity (or its generalization for non-barotropic flows) of the disc has a positive gradient (TL08, LT09, Tsang & Lai 2009b). In Chapter 3 we found that the inclusion of a finite toroidal field Bφ splits the CR into two magnetic (slow) resonances (MRs), where the wave frequency in the rotating frame of the fluid, ω − mΩ, matches the slow magnetosonic wave frequency. At the inner and outer MRs, the angular momentum flux carried by the wave undergoes significant change, indicating angular momentum (and energy) absorption at both MRs. One of the resonances absorbs positive angular momentum while the other emits positive angular momentum. Independent of the background flow vortensity gradient, the net angular momentum absorption across the resonance region is always reduced or becomes more negative in discs with Bφ 0 compared to the B = 0 case. This leads to a reduced growth rate of the p-modes. Our calculations showed that the hydrodynamically overstable inertial-acoustic modes can be completely stabilized by the toroidal field at relatively small field strengths (see Fig. 3.10)
For discs with a pure vertical field, the corotation resonance persists for the p-modes, with no additional magnetic resonances. This is understandable since p-modes have no vertical structure (kz = 0). With no bending of the field lines, the vertical field simply modifies the background pressure, or equivalently, changes the effective background sound speed. However, when Bz 0, the effective potential for wave propagation contains a second-order singularity term, in addition to the first-order singularity already present in a Bz = 0 disc.
For discs containing mixed (toroidal and vertical) magnetic fields, the corotation resonance is split into four resonances: in addition to the inner/outer MRs (already present in discs with pure toroidal fields), two Alfve´n resonances
157

appear, where ω−mΩ matches the local Alfve´n wave frequency. We showed that the effect of these additional resonances is to further reduce the super-reflection and the growth rate of the disc p-modes (see Fig. 3.11). Overall, the toroidal field has a larger effect on the corotational wave absorption and the mode growth rate than the vertical field.
As we summarized in the end of Chapter 3, our findings imply that in order to sustain the overstable disc p-modes, other mode excitation/destabilization methods are needed. One possibility is the accretion-ejection instability proposed by Tagger et al. Yu & Lai (2012) recently showed that this instability mechanism indeed could save the p-modes from being completely stabilized by disk toroidal B fields. This means that unstable disk p-modes can still be a promising candidate for HFQPO origin. In Chapter 4, we extended our study on disk p-modes into the non-linear regime. We employed numerical code PLUTO to simulate the non-linear evolution of unstable p-modes. We found very satisfactory agreement between our linear mode calculation (e.g. Chapter 3) and non-linear simulation in terms of frequencies of fastest growing modes. We showed that the mode frequencies obtained from linear mode calculation are quite robust. They vary by less than 4% even after the non-linear saturation.
All our analysis in Chapter 3 and Chapter 4 assumed a simple reflective inner disc boundary condition (radial velocity perturbation vanishes) whereas the condition at a real disc’s inner boundary could be much more complicated. In Chapter 5, we made the first attempt to consider more general inner disc boundary conditions. We studied the dynamics of a magnetosphere-disc system and we found that a magnetosphere inner to the disc, a scenario that could be realized in BH and NS X-ray binaries, can serve as the required “reflector” for
158

unstable disc p-modes. We also studied another type unstable modes – interface modes and the effects of magnetic fields, inner disk boundary locations on them. The large-scale (m = 1, 2, 3, . . .), overstable oscillation modes studied in Chapter 5 may provide an explanation for the high-frequency QPOs observed in NS and BH X-ray binaries (see Section 5.1). Obviously, the simplicity of our model setup precludes detailed comparison with the phenomenology of QPOs. But we note the following relevant features of disc inertial-acosutic modes and interface modes.
(1) The disc inertial-acoustic modes have frquencies ω = βdmΩd (with Ωd the disc rotation rate at rin), with βd < 1 (typically ∼ 0.5) depending on model parameters (such as disc sound speed) and boundary conditions. If the inner disk is located at the ISCO, the mode frequencies can be computed ab initio, and the results are generally consistent with the observations of the HFQPOs in blackhole X-ray binaries (see Lai & Tsang 2009; Tsang & Lai 2009c). The problem with these modes is that even a weak (sub-thermal) disc toroidal magnetic field can suppress their instabilities (see Section 5.3 and Fu & Lai 2011a). Although large-scale poloidal fields can enhance the instability under certain conditions (see Tagger & Pallet 1999; Tagger & Varniere 2006), it is not clear at this point which effects (disc toroidal field vs large-scale poloidal field) will dominate.
(2) The interface modes have frqeuencies ω = βimΩd, with βi ∼ 1. If rin ≃ rISCO, the implied QPO frequencies would be too high compared to observations. Of course, as discussed before (see Sections 5.1 and 5.4.3), it is certainly possible that rin is somewhat larger than rISCO in real systems, but then we would lose the predictive power of our calculations (since in our model rin is a free parameter). On the other hand, the interface modes are robustly unstable, driven by the RT
159

and KH instabilities, and by the corotation effect, especially when the GR effect is included. In most cases, the growth rates of the interface modes are much larger than the disc modes.
Overall, the results of Chapter 5 suggest that if the real accreting NS or BH systems can be approximated by our magnetosphere-disc model, the interface oscillations are more likely to provide an explnanation for the observed QPOs. This is reasonable for NS systems, and is consistent with spectral analysis of kHz QPOs in NS X-ray binaries (e.g., Gilfanov et al. 2003). For BH systems, this would require that the inner disc radius in the “transitional state” to be slightly larger than rISCO (e.g., Done et al. 2007; Oda et al. 2010).
Using the similar techniques as we did in studying disc flows, we also investigated the dynamics of other rotating astrophysical flow, such as accreting tori (Chapter 6) and rotating protoneutron stars (Chapter 7). In the case of accreting tori, we found that magnetic fields tend to suppress the well-known Papaloizou-Pringle instability if the torus is relatively thin; while for a relatively thick torus, sufficiently strong magnetic fields could give rise to a new instability, one that is absent when the torus is non-magnetized. In the case of protoneutron stars, our main finding is that toroidal magnetic fields reduce the growth rate of the low-T/|W| instability and increase the threshold T/|W| value above which the instability occurs.
Although our work in this thesis has made some progress in understanding the dynamics of various astrophysical flows, we note that the flow models employed are still too simplified to allow for any direction comparison with observations. For instance, in all chapters, the fluid perturbations have no vertical structures (kz = 0) while in real disks the perturbation very likely has finite
160

vertical wavenumber. How would our conclusions hold up when more general perturbation structures are considered? With finite kz and finite magnetic field, MRI will develop in the flow, how would it interact with the other modes that we studied here? Also, we have always neglecting viscosity, which could play a very important role in determining the dynamics of many fluid systems. Our inner disk boundary condition which does not incorporate inflow is still far from being realistic. Real accreting tori could be compressible while we only considered incompressible tori. Real stars are mostly spherical or near spherical while our stellar model for T/|W| instability study is an idealized cylinder. All these issues certainly need to be addressed in future investigations.
161

APPENDIX A GENERAL MAGNETIC FIELD PROFILES IN THE MAGNETOSPHERE

In the main text of Chapter 5 we adopt a special toroidal magnetic field profile Bφ ∝ r in the magnetosphere so that the perturbation equations have simple analytical solutions. For a more general magnetic field profile Bφ ∝ rq, with q > 1, analytical solutions are no longer attainable, and we must solve numerically the perturbation equations for both the disc (Eqs. [5.22]-[5.23]) and the magnetosphere (Eqs. [5.9]-[5.10]). To this end, we first need to derive the boundary condition at the the inner boundary of the magnetosphere, which is close to the center of the system (r = 0). This is done by requiring the solutions to be regular at small radius. As r → 0, we observe that ωAφ ∝ Bφ/r → 0. In this limit, Eqs. (5.9)-(5.10) then reduce to

dξr dr

=

1 ω˜ − 2mΩ

− r

ω˜ ξr

+

m2 r2 δh,

(A.1)

dδh dr

=

(ω˜ 2

−

rΩ2)ξr

+

2mΩ δh, rω˜

(A.2)

i.e, recovering the hydrodynamic equations1. These simplified equations can be

readily solved, leadimg to a relation between ξr and δh:

mδh

ξr

=

rω˜ (ω˜

+

. 2Ω)

(A.3)

With this regularity condition implemented at some finite yet small inner boundary rc, we carry out the integration in the magnetosphere towards the interface, while at the same time we also integrate the disc equations towards the same interface and employ shooting procedure to solve for the eigenvalues by requiring solutions from both regions satisfy the matching condition at the
1See Fu & Lai (2011c) for the derivation of the regularity condition for a even more general system (differentially rotating, magnetized, self-gravitating, etc.).

162

interface (i.e., continuity of ξr and ∆Π). Note that the eigenvalues of the system now have two components: the wave frequency ω, and the relation between the solutions in two regions, e.g., δh(rc)/δh(rout) or ξr(rc)/ξr(rout).
163

APPENDIX B PLANE-PARALLEL FLOW WITH A COMPRESSIBLE UPPER LAYER AND
A MAGNETIZED LOWER LAYER

Consider a system with two separate fluid layers in a constant gravitational field g = −gzˆ. The upper layer (z > 0) is a compressible fluid of density ρ = ρ+e−z/Hz and constant horizontal velocity u+, and the lower layer (z < 0) is an incompressible fluid of constant density ρ− and horizontal velocity u−. Here, Hz = c2s/g is the density scale height with cs being the sound speed, and both velocities are along the x-axis. This system was studied by Tsang & Lai (2009b), but here in Chapter 5 we we add a uniform horizontal magnetic field Bx in the lower fluid. As in Tsang & Lai (2009b) (correcting typos in their Appendix), we apply perturbations of the form eikx−iωt to both layers and solve for ω by demanding the Lagrangian displacement and Lagrangian total pressure perturbation be continuous across the interface between upper and lower fluids. Denoting

ω˜ + = ω − ku+,

kz2 = (k2 − ω˜ 2+/c2s), k˜ = ( 1 + 4Hz2kz2 − 1)/2Hz,

and ρ˜+ = ρ+k/k˜,

the final solution for ω can be written as

ω = k(ρ˜+u+ + ρ−u−) ± ρ˜+ + ρ−

−

k2(u+ − u−)2ρ˜+ρ− (ρ˜+ + ρ−)2

−

kg(ρ+ − ρ−) ρ˜+ + ρ−

+

k2B2x . 4π(ρ˜+ + ρ−)

(B.1)

The three terms under the square root clearly have different physics mean-

ings. The first term signifies the Kelvin-Helmholtz instability which disappears

164

when velocity shear across the interface vanishes; the second term depicts the Rayleigh-Taylor instability which only exists for non-zero density contrast; the last term, the only difference between Eq. (B.1) and Eq. (A14) in Tsang & Lai (2009b), describes the stabilizing effect of the magnetic field that lies along the direction of the perturbation wave vector. If we take the incompressible and hydrodynamic limit (i.e., cs → ∞, Hz → ∞ so that k˜ → k, ρ˜+ → ρ+, and Bx → 0), then the above equation reduces to Eq. (B12) in Li & Narayan (2004).
165

BIBLIOGRAPHY
[1] Abramowicz M. A., Kluzniak, W. 2001, A&A, 374, L19 [2] Altamirano D., Belloni T., 2012, arXiv:1201.2106 [3] Acheson, D. J. & Gibbsons, M. P., 1978, Phil. Trans. Roy. Soc. Lond. A 289,
459 [4] Akiyama, S., Wheeler, J. C., Meier, D. L., Lichtenstadt, I., 2003, ApJ, 584, 954 [5] Akiyama, S., Wheeler, J. C., 2005, ApJ, 629, 414 [6] Andersson, N., 2003, Class. Quantum Grav., 20, 105 [7] Ardeljan, N. V., Bisnovatyi-Kogan, G. S., Moiseenko, S. G., 2000, A&A, 355,
1181 [8] Ardeljan, N. V., Bisnovatyi-Kogan, G. S., Moiseenko, S. G., 2005, MNRAS,
359, 333 [9] Arras P., Blaes O.M., Turner N. J., 2006, ApJ, 645, L65 [10] Bachetti M., et al., 2010, MNRAS, 403, 1193 [11] Balbus, S. A., 2003, ARAA, 41, 555 [12] Balbus, S. A. & Hawley, J. F., 1991, ApJ, 376, 214 [13] Balbus, S. A. & Hawley, J. F., 1998, Rev. Mod. Phys., 70, 1 [14] Balmforth, N.J., & Korycansky, D.G., 2001, MNRAS, 326, 833 [15] Baruteau, C. & Masset, F., 2008, ApJ, 672, 1054 [16] Beckwith, K., Hawley, J. F., Krolik, J. H., 2008, MNRAS, 390, 21 [17] Beckwith, K., Hawley, J. F., Krolik, J. H., 2009, ApJ, 707, 428 [18] Begelman, M. C., Blandford, R. D. & Rees, M. J., 1984, Rev. Mod. Phys., 56,
255
166

[19] Belloni T. M., Motta S. E., Munoz-Darias T., 2011, BASI, 39, 409 [20] Bisnovatyi-Kogan G. S., Ruzmaikin A. A., 1974, Ap&SS, 28, 45 [21] Bisnovatyi-Kogan G. S., Ruzmaikin A. A., 1976, Ap&SS, 42, 401 [22] Blaes, O. M., 1987, MNRAS, 227, 975 [23] Blaes, O. M. & Glatzel, W., 1986, MNRAS, 220, 253 [24] Blaes, O. M., Sramkova, E., Abramowicz, M. A., Kluzniak, W. & Torkelsson,
U., 2007, ApJ, 665, 642 [25] Blokland, J. W. S., van der Swaluw, E., Keppens, R. & Goedbloed, J. P., 2005,
A&A, 444, 337 [26] Bondeson, A., Iacono, R. & Bhattacharjee, A., 1987, Phys. Fluids, 30, 2167 [27] Brandenburg, A., Nordlund, A., Stein, R. F. & Torkelsson, U., 1995, ApJ,
446, 741 [28] Brown, J. D., 2000, Phys. Rev. D, 62, 084024 [29] Burrows, A. et al., 2007, ApJ, 664, 416 [30] Camarda, K. D., Anninos, P., Fragile, P. C., Font, J. A., 2009, ApJ, 707, 1610 [31] Cazes, J. E., Tohline, J. E., 2000, ApJ, 532, 1051 [32] Centrella, J. M., New, K. C. B., Lowe, L, Brown, J. D., 2001, ApJ, 550, L193 [33] Cerda-Duran, P., Quilis, V., Font, J. A., 2007, Comp. Phys. Comm., 177, 288 [34] Chan, C-K. 2009, ApJ, 704, 68 [35] Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability, Ox-
ford University Press, Oxford [36] Chandrasekhar, S., 1969, Ellipsoidal Figures of Equilibrium, Yale Univer-
sity Press, New York
167

[37] Chandrasekhar, S., 1970, Phys. Rev. Lett., 24, 611 [38] Corvino, G., Rezzolla, L., Bernuzzi, S., De Pietri, R., Giacomazzo, B., 2010,
Class. Quantum Grav., 27, 114104 [39] Curry, C., Pudritz, R. E. & Sutherland, P. G., 1994, ApJ, 434, 206 [40] Curry, C. & Pudritz, R. E., 1995, ApJ, 453, 697 [41] Curry, C. & Pudritz, R. E., 1996, MNRAS, 281, 119 [42] Davis, S. W., Stone, J. M. & Pessah, M. E., 2010, ApJ, 713, 52 [43] de Val-Borro, et al. 2006, MNRAS, 370, 529 [44] De Villiers, J. P., Hawley, J. F., & Krolik, J. H., 2003, ApJ, 599, 1238 [45] Dimmelmeier, H., Ott, C. D., Marek, A., Janka, H.-T., 2008, Phys. Rev. D,
78, 064056 [46] Done C., Gierlinski M., Kubota A., 2007, Astron. Astrophys. Rev., 15, 1 [47] Drazin P. G., Reid W. H., 1981, Hydrodynamic Stability, Cambridge Univ.
Press, Cambridge [48] Duez, M., Foucart, F., Kidder, L. E., Ott, C. D. & Teukolsky, S. A., 2010,
Class. Quantum Grav., 27, 114106 [49] Etienne, Z. B., Liu, Y. T., Shapiro, S. L. & Baumgarte, T. W., 2009, Phys. Rev.
D, 79, 044024 [50] Feng, H., Rao, F. & Kaaret, P., 2010, ApJ, 710, L137 [51] Ferreira B. T., Ogilvie G. I., 2008, MNRAS, 386, 2297 [52] Fragile, P.C., Blaes, O., Anninos, P., Salmonson, J.D. 2007, ApJ, 668, 417 [53] Friedman, J. L., Schutz, B. F., 1978, ApJ, 222, 281 [54] Frieman, E. & Rotenberg, M., 1960, Rev. Modern Phys., 32, 898
168

[55] Fromang, S., 2010, A&A, 514, L5 [56] Fromang, S. & Papaloizou, J., 2007, A&A, 476, 1113 [57] Fu W., Lai D., 2009, ApJ, 690, 1386 [58] Fu W., Lai D., 2011a, MNRAS, 410, 399 [59] Fu W., Lai D., 2011b, MNRAS, 410, 1617 [60] Fu W., Lai D., 2011c, MNRAS, 413, 2207 [61] Fu W., Lai D., 2012, submitted [62] Gardiner, T.A., Stone, J.M., 2005, AIPC, 784, 475 [63] Ghosh P., Lamb, F. K., 1978, ApJ, 223, L83 [64] Gierlinski, M., Middleton, M., Ward, M. & Done, C., 2008, Nature, 455, 369 [65] Gilfanov M., Revnivtsev M., Molkov S., 2003, A&A, 410, 217 [66] Glatzel, W., 1987a, MNRAS, 225, 227 [67] Glatzel, W., 1987b, MNRAS, 228, 77 [68] Goldreich, P., Goodman, J. & Narayan, R., 1986, MNRAS, 221, 339 [69] Goodman, J., Narayan, R. & Goldreich, P., 1987, MNRAS, 225, 695 [70] Goldreich, P. & Tremaine, S., 1979, ApJ, 233, 857 [71] Goodman, J., 1993, ApJ, 406, 596 [72] Goossens, M., Hollweg, J. V. & Sakurai, T., 1992, Solar Physics, 138, 233 [73] Guan, X., Gammie, C. F., Simon, J. B. & Johnson, B. M., 2009, ApJ, 694, 1010 [74] Hachisu, I., 1986, ApJS, 61, 479 [75] Hameiri, E., 1981, J. Math. Phys., 22, 2080
169

[76] Hawley, J. F., Gammie, C. F. & Balbus, S. A., 1996, AJ, 464, 690 [77] Hawley, J.F., Guan, X., & Krolik, J.H., 2011, ApJ, 738, 84 [78] Heger, A., Woosley, S. E., Spruit, H. C., 2005, ApJ, 626, 350 [79] Heinemann, T. & Papaloizou, J. C. B., 2009, MNRAS, 397, 52 [80] Henisey, K. B., Blaes, O. M., Fragile, P. C., Ferreira, B. T., 2009, ApJ, 706, 705 [81] Hirose, S., Blaes, O. & Krolik, J. H., 2009, ApJ, 704, 781 [82] Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042 [83] Ikhsanov N. R., Pustil’nik L. A., 1996, A&A, 312, 338 [84] Ipser, J. R., 1994, ApJ, 435, 767 [85] Kato, S., 1990, PASJ, 42, 99 [86] Kato S., 2001, PASJ, 53, 1 [87] Kato S., 2003a, PASJ, 55, 257 [88] Kato S., 2003b, PASJ, 55, 801 [89] Kato S., 2008, PASJ, 60, 111 [90] Kato S., 2011a, PASJ, 63, 125 [91] Kato S., 2011b, PASJ, 63, 861 [92] Kato, S. & Fukue, J., 1980, PASJ, 32, 377 [93] Kato S., Fukue J., Mineshinge S., 2008, Black-Hole Accretion Disks, Kyoto
Univ. Press, Kyoto [94] Keppens, R., Casse, F. & Goedbloed, J. P., 2002, ApJ, 569, L121 [95] Kotake, K., Sawai, H., Yamada, S., Sato, K., 2004, ApJ, 608, 391
170

[96] Kluzniak, W. & Abramowicz, M. A., 2002, astro-ph/0203314 [97] Knobloch, E., 1992, MNRAS, 255, 25p [98] Kulkarni, A. K., et al., 2011, MNRAS, 414, 1183 [99] Kulkarni A. K., Romaonva M. M., 2008, MNRAS, 386, 673 [100] Kumar, S., Coleman, C. S. & Kley, W., 1994, MNRAS, 266, 379 [101] Lai, D., Shapiro, S. L., 1995, ApJ, 442, 259 [102] Lai D., Tsang D., 2009, MNRAS, 393, 979 [103] Lander, S. K., Jones, D. I., 2009, MNRAS, 395, 2162 [104] Latter, H. N. & Balbus, S. A., 2009, MNRAS, 399, 1058 [105] Lee, W. H., Abramowicz, M. A. & Kluziniak, W., 2004, ApJ, 603, L93 [106] Li L., Goodman J., Narayan R., 2003, ApJ, 593, 980 [107] Li L.-X., Narayan R., 2004, ApJ, 601, 414 [108] Liu, Y. T., Lindblom, L., 2001, MNRAS, 324, 1063 [109] Liu, Y. T., 2002, Phys. Rev. D, 65, 124003 [110] Longaretti, P. -Y., & Lesur, G., 2010, A&A, 561, A51 [111] Lovelace, R. V. E., Li, H., Colgate, S. A. & Nelson, A. F., 1999, ApJ, 513, 805 [112] Lovelace R. V. E., Turner L., Romanova M. M., 2009, ApJ, 701, 225 [113] Lovelace, R. V. E., Rothstein, D. M. & Bisnovatyi-Kogan, G. S., 2009, ApJ,
701, 885 [114] Lovelace R. V. E., Romanova M. M., Newman W. I., 2010, MNRAS, 402,
2575 [115] Lubow S. H., Spruit, H. C., 1995, ApJ, 445, 337
171

[116] Machida, M., Matsumoto, R. 2003, ApJ, 585, 429 [117] Machida, M., Nakamura, K. E., & Matsumoto, R., 2006, PASJ, 58, 193 [118] Meszaros, P., 2006, Rep. Prog. Phys., 69, 2259 [119] Middleton, M. & Done, C., 2010, MNRAS, 403, 9 [120] Middleton, M., Done, C., Ward, M., Gierlinski, M. & Schurch, N., 2009,
MNRAS, 394, 250 [121] Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [122] Mikhailovskii, A. B., Fridman, A. M., Churikov, A. P., Pustovitov, V. D. &
Smolyakov, A. I., 2009, Plasma Phys. Controlled Fusion 51, 045003 [123] Miller, K. A. & Stone, J. M., 2000, ApJ, 534, 398 [124] Miller M. C., Lamb F. K., Psaltis D., 1998, ApJ, 508, 791 [125] Montero, P. J., Zanotti, O., Font, J. A. & Rezzolla, L., 2007, MNRAS, 378,
1101 [126] Montero, P. J., Font, J. A. & Shibata, M., 2010, Phys. Rev. Lett., 104, 191101 [127] Moscibrodzka, M., Gammie, C. F., et al., 2009, ApJ, 706, 497 [128] Muto, T., Machida, M. N. & Inutsuka, Shu-ichiro, 2008, ApJ, 679, 813 [129] Narayan, R., Goldreich, P. & Goodman, J., 1987, MNRAS, 228, 1 [130] Noble, S. C., Krolik, J. H., Hawley, J. F., 2009, ApJ, 692, 411 [131] Noble, S. C., et al., 2011, ApJ, 743, 115 [132] Nowak M. A., Wagoner R. V., 1991, ApJ, 378, 656 [133] Nowak, M. A., Wagoner, R. V., 1992, ApJ, 393, 697 [134] Obergaulinger, M., Aloy, M. A., Dimmelmeier, H., Muller, E., 2006, A&A,
457, 209
172

[135] Obergaulinger, M., Cerda-Duran, P., Muller, E., Aloy, M. A., 2009, A&A, 498, 241
[136] Oda H., et al., 2010, ApJ, 712, 639 [137] Ogilvie, G. I. & Pringle, J. E., 1996, MNRAS, 279, 152 [138] Okazaki A. T., Kato S., Fukue J., 1987, PASJ, 39, 457 [139] O’Neill, S. M., Reynolds, C. S. & Miller, M. C., 2009, ApJ, 693, 1100 [140] O’Neill, S.M., Reynolds, C.S., Miller, C.M. & Sorathia, K. 2011, ApJ, 736,
107 [141] Ortega-Rodriguez, M., Silbergleit, A. S., & Wagoner, R. V. 2008, Geophys.
Astrophys. Fluid Dyn., 102, 75 [142] Ott, C. D., 2009, Class. Quantum Grav., 26, 063001 [143] Ott, C. D., Ou, S., Tohline, J. E., Burrow, A., 2005, ApJ, 625, L119 [144] Ott, C. D., Burrows, A., Thompson, T. A., Livne, E., Walder, R., 2006, ApJS,
164, 130 [145] Ott, C. D., Dimmelmeier, H., Marek, A., Janka, H.-T., Hawke, I., Zink, B.,
Schnetter, E., 2007, Phys. Rev. Lett., 98, 261101 [146] Ou, S., Tohline, J. E., 2006, ApJ, 651, 1068 [147] Paczynski B, Witta P. J., 1980, A&A, 88, 23 [148] Papaloizou, J. C. B. & Pringle, J. E., 1984, MNRAS, 208, 721 [149] Papaloizou, J. C. B. & Pringle, J. E., 1985, MNRAS, 213, 799 [150] Papaloizou, J. C. B. & Pringle, J. E., 1987, MNRAS, 225, 267 [151] Penna, R.F. et al., 2010, MNRAS, 408, 752 [152] Pessah, M. E., Chan, C. & Psaltis, D., 2006, MNRAS, 372, 183
173

[153] Pessah, M. E. & Goodman, J., 2009, ApJ, 608, L72 [154] Pickett, B. K., Durisen, R. H., Davis, G. A., 1996, ApJ, 458, 714 [155] Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numer-
ical Recipes. Cambridge Univ. Press, Cambridge [156] Pringle, J. & King, A., 2007, Astrophysical Flows, Cambridge University
Press, Cambridge [157] Rebusco, P., 2008, New Astron.Rev., 51, 855 [158] Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49 [159] Reynolds C. S., Miller M. C., 2009, ApJ, 692, 869 [160] Rezzolla, L., Baiotti, L., Giacomazzo, B., Link, D. & Font, J. A., 2010, Class.
Quantum Grav., 27, 114105 [161] Rezzolla, L., Yoshida, Si., Maccarone, & Zanotti, O. 2003, MNRAS, 344,
L37 [162] Romanova M. M., Kulkarni A. K., Lovelace R. V. E., 2008, ApJ, 673, L171 [163] Rothstein D. M., Lovelace R. V. E., 2008, ApJ, 677, 1221 [164] Saijo, M., Yoshida, S., 2006, MNRAS, 368, 1429 [165] Saijo, M., Shibata, M., Baumgarte, T. W., Shapiro, S. L., 2001, ApJ, 548, 919 [166] Sakurai, T., Goossens, M. & Hollweg, J. V., 1991, Solar Physics, 133, 227 [167] Schmidt, G., 1979, Physics of High Temperature Plasmas, Academic Press [168] Schnittman, J. D. & Rezzolla, L., 2006, ApJ, 637, L113 [169] Sekiya, M. & Miyama, S. M., 1988, MNRAS, 234, 107 [170] Shafee, R., et al., 2008, ApJ, 687, L25 [171] Shibata, M, Baumgarte, T. W., Shapiro, S. L., 2000, ApJ, 542, 453
174

[172] Shibata, M., Karino, S., Eriguchi, Y., 2002, MNRAS, 334, L27 [173] Shibata, M., Karino, S., Eriguchi, Y., 2003, MNRAS, 343, 619 [174] Shibata, M., Sekiguchi, Y.-I., 2005, Phys. Rev. D, 71, 024014 [175] Shu, F. H., 1992, Gas Dynamics, University Science Books [176] Silbergleit, A. S., & Wagoner, R. V. 2008, ApJ, 680, 1319 [177] Silbergleit, A. S., Wagoner, R. V., & Ortega-Rodriguez, M. 2001, ApJ, 548,
335 [178] Simon, J. B., Hawley, J. F., Beckwith, K., 2009, ApJ, 690, 974 [179] Sorathia, K. A., Reynolds, C. S. & Armitage, P. J., 2010, ApJ, 712, 1241 [180] Spruit H. C., Stehle R., Papaloizou J. C. B., 1995, MNRAS, 275, 1223 [181] Stella L., Vietri M., Morsink S. M. 1999, ApJ, 524, L63 [182] Stergioulas, N., 2003, Liv. Rev. Rel., 6, 3 [183] Stone, J. M., Hawley, J. F., Gammie, C. F. & Balbus, S. A., 1996, ApJ, 463,
656 [184] Strohmayer, T. E. & Mushotzky, R. F., 2009, ApJ, 703, 1386 [185] Sramkova, E., Torkelsson, U. & Abramowicz, M. A., 2007, A&A, 467, 641 [186] Swank J., 1999, Nuclear Phys. B, Proc. Suppl., 69, 12 [187] Tagger M., Pellat R. 1999, A&A, 349, 1003 [188] Tagger M., Varniere P. 2006, ApJ, 652, 1457 [189] Terquem, C., 2003, MNRAS, 341, 1157 [190] Terquem, C. & Papaloizou, J. C. B. 1996, MNRAS, 279, 767 [191] Tohline, J. E., Durisen, R. H., McCollough, M., 1985, ApJ, 298, 220
175

[192] Toman, J., Imamura, J. N., Pickett, B. K., Durisen, R. H., 1998, ApJ, 497, 370
[193] Tomimura, Y., Eriguchi, Y., 2005, MNRAS, 359, 1117 [194] Tsang D., Lai D., 2008, MNRAS, 387, 446 [195] Tsang D., Lai D., 2009a, MNRAS, 393, 992 [196] Tsang D., Lai D., 2009b, MNRAS, 396, 589 [197] Tsang D., Lai D., 2009c, MNRAS, 400, 470 [198] van der Klis M., 2006, in Lewin W. H. G., van der Kils M., eds, Compact
Stellar X-ray sources. Cambridge Univ. Press, Cambridge [199] Varniere, P. & Tagger, M., 2002, A&A, 394, 329 [200] Velikhov, E. P., 1959, Sov. Phys. JETP 36, 995 [201] Wagoner R. V., 1999, Phys. Rep., 311, 259 [202] Watts, A. L, Andersson, N., Jones, D. I., 2005, ApJ, 618, L37 [203] Woosley, S. E., 1993, ApJ, 405, 273 [204] Yamada, S., Sawai, H., 2004, ApJ, 608, 907 [205] Yu, C. & Lai, D., 2012, in prep. [206] Yu, C. & Li, H., 2009, ApJ, 702, 75 [207] Zhang, H. & Lai, D. 2006, MNRAS, 368, 917
[208] Zurek, W. H. & Benz, W., 1986, ApJ, 308, 132
176