Global 3 D MHD Simulations of Waves in Accretion Discs

We discuss results of the first global 3D MHD simulations of warp and density waves in accretion disks excited by a rotating star with a misaligned dipole magnetic field. A wide range of cases are considered. We find for example that if the star’s magnetosphere corotates approximately with the inner disk, then a strong one-arm bending wave or warp forms. The warp corotates with the star and has a maximum amplitude (|zw|/r ∼ 0.3) between the corotation radius and the radius of the vertical resonance. If the magnetosphere rotates more slowly than the inner disk, then a bending wave is excited at the disk-magnetosphere boundary, but it does not form a large-scale warp. In this case the angular rotation of the disk [Ω(r)] has a maximum as a function of r so that there is an inner region where dΩ/dr > 0. In this region we observe radially trapped density waves in approximate agreement with the theoretical prediction of a Rossby wave instability in this region.


Introduction
Different types of accreting stars have dynamically important magnetic fields, for example, young T Tauri stars [1], accreting millisecond pulsars [2], and accreting white dwarfs [3].Doppler tomography of classical T Tauri stars (CTTSs) shows that in spite of the complexity of the magnetic field at the surface of the star, the dipole component is likely the dominant component at the disk-magnetosphere boundary [4].
A rotating star with a tilted dipole (or other non-axisymmetric) magnetic field exerts a periodic force on the inner part of the accretion disk which acts to excite different types of waves that propagate to larger distances.This force can excite density waves in the disk [5] and/or vertical bending waves in the disk [6].These waves can have important observational consequences such as periodically obscuring the light from the star (e.g., [7]) and periodically modulating the radiation from the inner disk and the star possibly giving rise to observed quasi-periodic oscillations (QPOs) in low-mass X-ray binary systems [2].Theoretical investigations of the small amplitude linearized waves in disks has led to understanding that both types of waves may be strongly enhanced or damped at radii corresponding to resonances (e.g., [8]).There are two types of resonances.One is a horizontal resonance (Lindblad resonance) and the other is a vertical resonance.Another important resonance is associated with corotation radius r cr where the angular velocity of the disk matches the angular velocity of the external force.The linear theory is of course not applicable when the amplitude of the waves is large.For this limit simulations are needed.Here, we summarize results from the first global 3D magnetohydrodynamic simulations of the waves excited in an accretion disk by a rotating star with a dipole magnetic field misaligned with the star's rotation axis (which is aligned with the disk's rotation axis at large distances) [9].In Sect. 2 we outline the theory.Sect. 3 gives sample results from our MHD simulations and Sect. 4 briefly discusses the Rossby waves found in our simulations.

Linearized Modes of Thin Disks
The small amplitude waves in a hydrodynamic disk are characterized by the perturbations of the different fluid variables, Q = p, ρ, v r , .., which can be expanded as (e.g., [10]): Here, ω is the angular frequency of the wave, m = 0, 1, 2... is the number or arms in the azimuthal direction, H n , n ≥ 0 are Hermite polynomials (with n = 0, 1, 2, ..), and h(r) the half-thickness of the disk at the radius r.We restrict our attention to the low order modes m = 1, 2 and n = 0, 2.

03004-p.2
Instabilities and Structures in Proto-Planetary Disks Linearization of the equations of motion [8] leads to where ψ = δp/ρ 0 is the perturbation of the enthalpy.For the bending modes ψ is the displacement of the disk center of mass from the equatorial plane, z w .Also in this equation, ω ≡ ω − mΩ is the Doppler-shifted angular frequency of the wave seen by observer orbiting at the disk's angular velocity Ω(r), Ω ⊥ is the frequency of oscillations perpendicular to the disk, κ 2 = r −3 d(r 4 Ω 2 )/dr is the square of the radial epicyclic frequency, k r is the radial wavenumber, and c s is the sound speed in the disk.
Considering the local free wave solution in the form in the WKBJ approximation [8], for (k r r) 2  1, (k r h) 2 1.Thus d 2 ψ/dr 2 = −k 2 r ψ in equation ( 2) and one gets the dispersion relation Propagation of the perturbations is possible in regions where ( ω2 Part of the range where ω2 < κ 2 describes density waves (inertial-acoustic) waves and part of that with ω2 > nΩ 2 ⊥ describes bending waves.The WKBJ approximation breaks down in the vicinity of points where the coefficient in front of ψ in equation ( 2) changes sign or has a singularity.At the points where ω = ±κ there is an outer (+sign) or an inner (-sign) Lindblad resonance.At points where ω = ±nΩ ⊥ , there are vertical resonances.For ω = 0 there is a corotation resonance.The behavior of the different modes is shown in Figure ??.
An important class of waves in the disk are the waves 'driven' by the distributed force from the rotating tilted magnetosphere of the star.The interaction of the magnetosphere with the disk is complex, non-linear and non-stationary process.This force can be written as where Ω * is the angular velocity of the star.In linear approximation, we suggest that the m−harmonic of this force excites the m−arm wave in the disk.The frequency of this force is mΩ * .Then, for a Keplerian disk, the condition for the Lindblad resonances ω = ±Ω is mΩ * − mΩ = ±Ω, or, r LR = r cr (1 ± 1/m) 2/3 , where r cr = (GM/Ω 2 * ) 1/3 is the corotation radius.For bending waves, n = 1, the condition for vertical resonances ω = ±Ω is the same, and hence vertical resonances are located at the same radii.

3D MHD Simulations
We have solved the 3D MHD equations with a Godunov-type code in a reference frame rotating with the star using the "cubed sphere" grid.The model has been described earlier in a series of previous papers ( [11], [9]).Hence, we describe it only briefly here.For initial conditions we assume that the rotating, perfectly conducting, magnetized star is surrounded by an accretion disk and a corona.The disk is cold and dense, while the corona is hot and rarefied, and at the reference point (the inner edge of the disk in the disk plane), T c = 100T d , and ρ c = 0.01ρ d , where the subscripts 'd' and 'c' denote the disk and the corona.The disk and corona are initially in rotational hydrodynamic equilibrium [12].The disk is relatively thin with the half-thickness to radius ratio h/r ≈ 0.1.The grid consists of N r concentric spheres, where each sphere represents an inflated cube.Each of the six sides of the inflated cube has an N × N curvilinear Cartesian grid.The whole grid consists of 6 × N r × N 2 cells.The typical grid used in simulations has N 2 = 61 2 angular cells in each block.The typical number of grid cells in the radial direction is N r = 140.Figure 2 shows sample results for the bending and density waves found in our 3D MHD simulations [9].
Figure 3 shows the power spectral densities (PSD) for the m = 1 warp wave (z w ) and the m = 2 density wave (σ).The left-hand panel shows that warp wave rotates at the angular frequency of the star and that it extends radially outward and inward from the corotation radius.It is in approximate agreement with the theory of [13].The right-hand panel shows the PSD of the density waves.Note in particular the prominent feature on the inner part of the rotation curve where dΩ/dr < 0. We identify this feature as a Rossby wave as discussed below.

Radially Trapped Rossby Modes
A radially localized instability of thin disks know as the 'Rossby wave instability' (RWI) may occur near regions where there is an extremum (as a function of r) of the potential vorticity ẑ•(∇×v)/Σ, where v is the disk velocity and Σ is the disk's surface mass density ( [14], [15], [16]).Of particular interest here is the RWI in the non-Keplerian region of a disk [17].which occurs when a disk encounters the magnetosphere of a slowly rotating star.An example of such a rotation curve derived from simulations [17] is shown in Figure 4.For given profiles of v φ , ρ, etc. and a given m, the value of the growth rate is found by solving for the 'ground state' solution of Schrödinger-like equation for the enthalpy where U(r|ω) is an effective potential shown for a sample case in Figure 4 from [17].The ground state is found by varying both the real and imaginary parts of the mode frequency ω so as to give the ψ(r) most deeply trapped in the potential well.For typical profiles and m = 1, the growth rate is ω i ∼ 0.1(v φ /r) R and the real part of the frequency is ω r ≈ (v φ /r) R , where the R−subscript indicates evaluation at the resonant radius r R .The radial width of the mode is ∆r/r R ∼ 0.05.The linear growth of the Rossby wave is predicted to saturate when the azimuthal frequency of trapping of a fluid particle in the trough of the wave ω T grows to a value equal to the growth rate ω i as discussed by [17].Recently, the saturation of the RW growth has been studied in a large set simulations and the saturation in most cases is found to correspond to ω T ≈ ω i [18].

Figure 1 .
Figure 1.Panels a and b show the radial dependence of the radial wavenumber for the one-armed (m = 1) and two-armed (m = 2) bending waves (n = 1).Here, r cr is the corotation radius, CR indicates the position of the corotation resonance, OVR the outer vertical resonance, and IVR the inner vertical resonance.Panels c and d show the radial wavenumber for the one-armed and two-armed in plane density waves (n = 0), where ILR and OLR are the inner and outer Lindblad resonances.This figure assumes a thin Keplerian disk with Ω = GM/r 3 and disk half-thickness (h) has h/r = c s /(rΩ) = 0.05.

Figure 2 .
Figure 2. Views of the bending and density waves shown at time t = 80 measured in orbits at the inner disk radius.The left two panels show 3D views of the isodensity surface of a warp from the side and face-on orientations.The right two panels show the bending (where z w is the vertical displacement of the center of mass of the disk) and the density wave (where σ is the disk's surface mass density).Also shown are the locations of resonances: r CR is the corotation resonance; r OVR is the outer vertical resonance; and r ILR and r OLR are the inner and outer Lindblad resonances.

Figure 3 .
Figure 3. Power spectral densities (PSD) for the m = 1 warp wave and the m = 2 in plane density wave (σ).The thick red curves show the average angular velocity Ω of the disk matter.At small radii the angular velocity decreases owing to the slow rotation rate of the star.This give dΩ/dr < 0 at small r.

Figure 4 .
Figure 4.The left-hand panel shows the midplane radial dependencies of the azimuthal disk velocity v φ , magnetic field B z , and density ρ, where v K is the Keplerian velocity.The vertical arrow indicating r R is the location of the bottom of a potential well U(r) shown in the right-hand panel from[17]