A&A 439, 443-459 (2005)
DOI: 10.1051/0004-6361:20041511

Forced oscillations in magnetized accretion disks and QPOs

J. Pétri

Astronomical Institute, University of Utrecht, PO Box 80000, 3508 TA Utrecht, The Netherlands

Received 22 June 2004 / Accepted 12 April 2005

Abstract
Quasi-periodic oscillations (QPOs) have been observed in accretion disks around neutron star, black hole, and white dwarf binaries with frequencies ranging from a few 0.1 Hz up to 1300 Hz. Recently, a correlation between their low- and high-frequency components was discovered and fitted with a single law, irrespective of the nature of the compact object. That such a relation holds over 6 orders of magnitude strongly supports the idea that the physical mechanism responsible for these oscillations should be the same in all binary systems. We propose a new model for these QPOs based on forced oscillations induced in the accretion disk due to the stellar magnetic field. First, it is shown that a magnetized accretion disk evolving in a rotating nonaxisymmetric magnetic field anchored to a neutron star will be subject to three kinds of resonances: a corotation resonance, a Lindblad resonance due to a driving force, and a parametric resonance due to the time varying epicyclic frequencies. The asymmetric part of the field is assumed to contain only one azimuthal mode $m\ge1$. We focus on the m=1 disturbance, which is well studied for an inclined dipolar rotator; but our results are general and easily extend to m>1. However, the radial location of the resonances will be affected by this number m. For instance, with an m=1 asymmetric structure, the resonances reach regions very close to the innermost stable circular orbit (ISCO) and can account for observations of kHz-QPOs at frequencies as high as 1200-1300 Hz. If we replace the dipolar by a higher order multipolar component, the resonance location is shifted to larger radius, implying lower QPO frequencies. To compare the MHD situation with the hydrodynamical case, we also consider an m=2 component in the magnetic perturbation in order to prove that, at least in the linear regime, the conclusions in both cases are the same. In the second part of the paper, we focus on the linear response of a thin accretion disk, developing the density perturbation as the sum of free wave solutions and non-wavelike disturbances. In the last part, we show results of 2D numerical simulations of a simplified version of the accretion disk consisting of a column of plasma threaded by a vertical magnetic field. These simulations are performed for the Newtonian gravitational potential, as well as for a pseudo-general relativistic potential, which enables us to explore the behavior of the resonances around both rotating neutron stars and black holes. We found that the density perturbations are only significant in the region located close to the inner edge of the disk near the ISCO where the magnetic perturbation is maximal. They induce fluctuations in the density which persist over the whole time of the simulations and are closely related to the spin of the magnetic perturbation. It is argued that the nearly periodic motion induced in the disk will produce high quality factor QPOs.

Key words: accretion, accretion disks - magnetohydrodynamics (MHD) - instabilities - methods: analytical - methods: numerical - stars: neutron

1 Introduction

To date, more than 20 accreting neutron stars in Low Mass X-Ray Binaries (LMXBs) are known to exhibit rapid variability in their X-ray fluxes. These high frequency quasi-periodic oscillations (QPOs) show very strong similarities in their shape, as well as in their amplitude, and span frequencies from 300 Hz to about 1300 Hz (see van der Klis 2000, for a review).

The comparison between high frequency QPOs (HFQPOs) in black holes and neutron stars shows a scaling proportional to 1/M_*, where M_* is the mass of the compact object, suggesting that these oscillations are due to the orbital motion in the accretion disk near the innermost stable circular orbit (ISCO) (McClintock & Remillard 2003). Moreover, the QPO frequencies depend on the spectral state of the source and are correlated with the accretion rate. QPOs are observed in the tail of the power law spectrum indicating that they are not due to thermal emission but probably related to Comptonization of soft photons. The saturation of the spectral index was explained by Titarchuk & Fiorito (2004) in the context of coronal heating and phase transition.

Therefore, any QPO model must take not only oscillations in the accretion disk into account but also the way photons propagate to the observer, including photon scattering and general relativistic effects such as beaming, Doppler shift, and aberration. This was investigated by Schnittman & Bertschinger (2004), who used exact integration of the geodesic motion of hot spots in the Kerr spacetime and looked at the harmonic content of the signal received by an observer in the asymptotically flat spacetime.

In black hole binaries, the HFQPOs appear in pairs, as in the neutron star system, but their frequency remains uncorrelated to the accretion rate. They act as a voiceprint at a ratio 3:2 that favors a resonance mechanism. For example a parametric resonance between vertical and horizontal epicyclic frequencies was suggested by Kluzniak et al. (2004). Another mechanism proposed by Lee et al. (2004) is resonance due to an external driven force.

The physics of accretion onto a magnetized neutron star has been described by Ghosh & Lamb (1978, 1979a,b). In this picture close to the star's surface, the magnetic field channels matter to the polar caps and exerts a torque on it, which accelerates or decelerates the neutron star depending on the angular momentum exchange rate. This is the standard model of magnetized accretion onto a neutron star. However, in order to explain pulsating X-ray sources, Anzer & Börner (1980) propose an alternative view in which the magnetic moment of the neutron star and the accretion disk rotation axis are perpendicular. They also studied the influence of the Kelvin-Helmholtz instability on such a configuration (Anzer & Börner 1983).

When the magnetic moment is misaligned with respect to the spin axis of the star, Vietri & Stella (1998) show that some inhomogeneities, treated as diamagnetic blobs, could be lifted off the disk and suffer the Lense-Thirring precession.

Lai (1999) has identified a magnetically driven warping and precession instability that occurs in the inner region of an accretion disk. This occurs by interaction between the surface current density in the disk with the stellar magnetic field and produces motions at frequencies much lower than the orbital one. It can, therefore, account for the low frequency QPOs (LFQPOs), as does the usual Lense-Thirring precession (Stella & Vietri 1998). Magnetic precession, combined with the relativistic dragging effect, was applied to weakly magnetized neutron stars to explain these LFQPOS in LMXBs, (Shirakawa & Lai 2002a). The study of the global magnetic warping/precession process for strongly magnetized neutron stars can also explain the millihertz QPOs in some systems, see for example Shirakawa & Lai (2002b). The Rayleigh-Taylor instability associated with Rossby waves and the rotational splitting gives another self-consistent description of the QPO phenomenon (Titarchuk 2003).

It was realized by Kato (2001a,b) that non-axisymmetric g-mode oscillations can be trapped in a thin relativistic disk. These are excited by the corotation resonance. The effect of non-linear coupling between oscillations in the disk and warp is examined by Kato (2004). However, his study does not include the rotation of the compact object.

Recent observations of accretion disks orbiting white dwarfs, neutron stars, or black holes have shown a strong correlation between the low and high frequency QPOs, LFQPOs and HFQPOs respectively, (Psaltis et al. 1999; Mauche 2002; Belloni et al 2002). This relation holds over more than 6 orders of magnitude and strongly supports the idea that the QPO phenomenon is a universal physical process independent of the nature of the compact object. Indeed, the presence or absence of a solid surface, a magnetic field, or an event horizon plays no relevant role in the production of X-ray variability (Wijnands 2001). This correlation has instead been explained in terms of the centrifugal barrier model of Titarchuk et al. (2002) where the LFQPOs are associated with magneto-acoustic waves, and the HFQPOs correspond to the Keplerian motion in the disk.

In this paper we propose a new resonance in the star-disk system arising from the response of the accretion disk to a non-axisymmetric rotating magnetic field. The paper is organized as follows. In Sect. 2, we describe the initial equilibrium state where the accretion disk is embedded in an axisymmetric magnetic field anchored to the neutron star. Furthermore, we assume that the magnetic moment is aligned with the rotation axis of the star. We then estimate the field perturbation induced by a small inclination of the magnetic moment. In Sect. 3, we look for the resonance conditions and elucidate the meaning of the resonances. For a thin and weakly magnetized accretion disk, assuming no warp and only motion in the orbital plane of the disk, we investigate the precise linear behavior of the disk which is presented in Sect. 3.2. In Sect. 5, we perform 2D numerical MHD simulations of the disk and show that some resonances persist on a long timescale with a narrow radial extension. This is done for Newtonian, as well as for pseudo-general relativistic gravitational potential describing the Kerr geometry in a Newtonian way adapted to a rotating neutron star. Conclusions are drawn in Sect. 6.

   
2 The initial conditions

2.1 The stationary state

In the equilibrium state, the accretion disk evolves in a perfectly spherically symmetric gravitational field and an axisymmetric dipolar magnetic field, both due to the compact object. We assume that the magnetic moment of the star and its rotational axis are aligned. The disk is assumed to lie in the equatorial plane of the star so that it experiences an axisymmetric stationary field. By adopting a cylindrical coordinate system labeled by  $(r,\varphi,z)$, we can therefore describe all physical quantities with only a (r,z)dependence. Neglecting the matter inflow because in the thin disk approximation  $v_r\ll v_\varphi$ and because we only focus on oscillations around a stationary state, the density $\rho$, pressure p, current density $\vec{j}$, velocity $\vec{v}$, and magnetic field $\vec{B}$ at equilibrium are thus given by:

 

$$\rho \rho(r,z)$$ = (1)

p = p(r,z) (2)

$$\vec{v} r ~ \Omega(r,z) ~ \vec{e}_\varphi$$ = (3)

$$\vec{B} B_r(r,z) ~ \vec{e}_r + B_z(r,z) ~ \vec{e}_z$$ = (4)

$$\vec{j} j_\varphi(r,z) ~ \vec{e}_\varphi.$$ = (5)

The magnetic field has two different origins:

• one component $\vec{B}_*$ is attached to the neutron star and corresponds to the rotating stellar dipolar magnetic field;
• the second component  $\vec{B}_{\rm d}$ of the magnetic field is induced by the fluid motion in the accretion disk and follows from Ampere's law  $\bf {rot}~\vec{B}_{\rm d}=\mu_0~\vec{j}$.

The total magnetic field is then simply the sum of these two components  $\vec{B}=\vec{B}_*+\vec{B}_{\rm d}$.

In the radial direction, the gravitational attraction from the compact object $\vec{g}$ is balanced by the centrifugal force, the Lorentz force, and the pressure gradient, whereas in the vertical direction we simply have hydromagnetic equilibrium:

  

$$-\rho ~ \frac{v_\varphi^2}{r} = \rho ~ g_r - \frac{\partial p}{\partial r} + j_\varphi ~ B_z \rho ~ g_r - \frac{\partial}{\partial r} ~ \left( p + \frac{B_z^2}{2~\mu_0} \right)$$ = (6)

$$\rho ~ g_z - \frac{\partial p}{\partial z} - j_\varphi ~ B_r$$ = 0. (7)

To specify a given equilibrium state, we prescribe the initial density and current, or, equivalently, the magnetic field in the disk and the strength of the dipolar magnetic field of the neutron star. Assuming a thin accretion disk and a low $\beta$-plasma parameter, defined as

$$\beta = \frac{p}{B^2/2~\mu_0},$$ (8)

the pressure gradient and Lorentz forces will be negligible compared to the gravitational force. The motion will then only deviate slightly from Keplerian rotation. We then perturb the system by adding a small inclination angle $\chi$ between the magnetic and rotational axis of the star. As a consequence, and due to the magnetic perturbation, the disk will leave its equilibrium state and start to oscillate around the equatorial plane. We next give the expression for the stellar magnetic field perturbation.

2.2 Perturbed dipolar magnetic field of a rotating star

In the case of an aligned rotator, $\chi=0$, the magnetic field is axisymmetric. Therefore, in the disk plane, there is only a vertical component B_z. When the magnetic moment $\vec{\mu}$ is inclined by an angle $\chi\ne0$ relative to rotation axis  $\vec{\Omega}_*$, an additional magnetic component dragged by the star's rotation arises. This perturbation is responsible for the oscillation in the disk. We derive its expression below. The magnetic moment can be developed into a parallel (along $\vec{e}_z$) and a perpendicular (along  $\vec{e}_r^*$) component such that

$$\vec{\mu} = \mu ~ \sin\chi ~ \vec{e}_r^* + \mu ~ \cos\chi ~ \vec{e}_z,$$ (9)

where  $\vec{e}_r^*$ is a radial unitary vector corotating with the star. In the Cartesian coordinates of a distant observer at rest  $(\vec{e}_x, \vec{e}_y, \vec{e}_z)$, it can be expressed as

$$\vec{e}_r^* = \cos(\Omega_*~t) ~ \vec{e}_x + \sin(\Omega_*~t) ~ \vec{e}_y$$ (10)

assuming that at t=0 $\vec{e}_r^* =\vec{e}_x$. We adopt a particular geometric configuration corresponding to the situation where the disk is contained in the star's equatorial plane. We could generalize to the case where magnetic moment, disk, and star rotation axis are all misaligned. However, the rotational frequency of the magnetic perturbation, which is its most important property, would not be affected.

The angular rotation is assumed to be equal to  $\vec{\Omega}_* = \Omega_*~\vec{e}_z$. The stellar dipolar magnetic field has the usual expression derived from the vector potential

$$\vec{A}_* = \frac{\mu_0}{4~\pi} ~ \frac{\vec{\mu}\wedge\vec{r}}{r^3}\cdot$$ (11)

Therefore, the perturbation is  $\delta\vec{A}_* = \frac{\mu_0}{4~\pi} ~ \mu ~ \sin\chi ~ \frac{\vec{e}_r^*\wedge\vec{r}}{r^3}$. It represents a perpendicular rotator with magnetic moment  $\mu_\perp=\mu~\sin\chi$. Expressed in the cylindrical coordinate system, the components of the perturbation of the magnetic field are given by:

  

$$\delta B_*^r \frac{\mu_0}{4 ~ \pi} ~ \frac{\mu ~ \sin \chi}{(r^2+z^2)^{3/2}} ~ \left[ \frac{3~r^2}{r^2+z^2} - 1 \right] ~ \cos (\varphi-\Omega_*~t)$$ = (12)

$$\delta B_*^\varphi \frac{\mu_0}{4 ~ \pi} ~ \frac{\mu ~ \sin \chi}{(r^2+z^2)^{3/2}} ~ \sin (\varphi-\Omega_*~t)$$ = (13)

$$\delta B_*^z \frac{\mu_0}{4 ~ \pi} ~ \frac{\mu ~ \sin \chi}{(r^2+z^2)^{3/2}} ~ \frac{3~r~z}{r^2+z^2} ~ \cos (\varphi-\Omega_*~t).$$ = (14)

We see that in the equatorial plane the perturbation is $\delta B_*^z = ~0$. Each component has a time function that depends on $\{\cos$/ $\sin\}(\varphi-\Omega_*~t)$. The nature of the driving is therefore an m=1 azimuthal mode with an excitation frequency at a rate $\Omega _*$. This has to be contrasted with the m=2 mode and, consequently, excitation frequency of $2\Omega_*$ produced by the quadrupolar gravitational field (the lowest order perturbation in gravity for a hydrodynamical accretion disk for example). This situation will be studied in detail in a forthcoming paper. The nature of the perturbation has some consequences for the derivation of the resonance conditions, as shown in the following Sect. 3.

   
3 Linear analysis

The distorted magnetic field described above will kick the disk out of its equilibrium state. We study its linear response by treating the asymmetric part of the magnetic field as a small perturbation of the equilibrium state prescribed in Sect. 2. We start with the ideal MHD equations of the accretion disk with adiabatic motion given by:

    

$$\frac{\partial\rho}{\partial t} + {\rm div}~(\rho~\vec{v}) = 0$$ (15)

$$\frac{\partial\vec{v}}{\partial t} + (\vec{v}\cdot\vec{\nabla}~) ~ \vec{v} = \vec{g} - \frac{\vec{\nabla}~p}{\rho} + \vec{j} \wedge \vec{B}$$ (16)

$$\frac{D}{Dt} \left( \frac{p}{\rho^\gamma} \right) = 0$$ (17)

$$\frac{\partial\vec{B}}{\partial t} = \bf {rot}~( \vec{v} \wedge \vec{B} )$$ (18)

$$\bf {rot}~\vec{B} = \mu_0 ~ \vec{j}.$$ (19)

All quantities have their usual meaning: $\rho$ density of mass in the disk, $\vec{v}$ the velocity of the disk, p the gaseous pressure, $\gamma$ the adiabatic index, $\vec{g}$ the spherically symmetric gravitational field, and $\vec{B}$ the total magnetic field.

3.1 Lagrangian displacement

We then perturb the equilibrium state from Sect. 2 by introducing the Lagrangian displacement $\vec{\xi}$. Expanding to first order we get for the Eulerian perturbations of the magnetohydrodynamical quantities attached to the disk:

$$\delta\rho = - ~ {\rm div}~(\rho ~ \vec{\xi} ~ )$$ (20)

$$\delta\vec{v} = \frac{\partial\vec{\xi}}{\partial t} + \vec{v}\cdot\vec{\nabla}~\vec{\xi} - \vec{\xi} \cdot \vec{\nabla}~\vec{v}$$ (21)

$$\delta p = - ~ \gamma ~ p ~ {\rm div}~\vec{\xi} - \vec{\xi}\cdot\vec{\nabla}~p$$ (22)

$$\delta\vec{B} = \vec{B} \cdot \vec{\nabla}~\vec{\xi} - \vec{B} ~ {\rm div}~\vec{\xi} - \vec{\xi} \cdot \vec{\nabla}~\vec{B}.$$ (23)

For a detailed derivation of these perturbed quantities and their associated boundary conditions as well as for the equations described in this section, we refer the reader to Appendix A. By also making allowance for a perturbation in the neutron star magnetic field and following the Frieman-Rotenberg analysis (Frieman & Rotenberg 1960), the Lagrangian displacement satisfies a second order linear partial differential equation given by:

 

$$\rho ~ \frac{D^2\vec{\xi}}{Dt^2} - {\rm div}~\left[ \rho ~ \vec{\... ... ) \vec{v} \right] - \nabla\Pi - \frac{1}{\mu_0} ~ \vec{B} \cdot \nabla \vec{Q}$$

$$\quad- \frac{1}{\mu_0} ~ \vec{Q} \cdot \nabla \vec{B} + {\rm div}~( \rho ~ \vec{\xi} ) ( \vec{g} + \delta\vec{g} ) - \rho ~ \delta\vec{g} =$$

$$\quad \frac{1}{\mu_0} ~ \left[ \bf {rot}~( \vec{B} + \vec{Q} + \d... ...ta \vec{B}_* + \bf {rot}~\delta \vec{B}_* \wedge ( \vec{B} + \vec{Q} ) \right].$$ (24)

We have introduced the convective derivative by  $D/Dt = \partial_t + \Omega~\partial_\varphi$. The perturbation in the magnetic field is represented by the vector  $\vec{Q} = \bf {rot}~(\vec{\xi} \wedge \vec{B})$. The scalar  $\Pi=\gamma~p~{\rm div}~\vec{\xi} + \vec{\xi} \cdot \vec{\nabla} p - \frac{1}{\mu_0} ~ \vec{B} \cdot \vec{Q}$ corresponds to the opposite of the total (gaseous+magnetic) pressure perturbation. Compared to the equation obtained by Frieman & Rotenberg (1960), we have an additional term in the Lorentz force induced by the stellar magnetic perturbation depicted by the expression enclosed by the brackets after the $\frac{1}{\mu_0}$factor, terms containing  $\delta \vec{B}_*$. We also allow for a degree of freedom in the choice of an arbitrary gravitational perturbation induced by the disk  $\delta\vec{g}$. However, for the rest of this paper, we only focus on the consequences of an asymmetry in the magnetic field, leaving the gravity perturbation for future study.

We make a clear distinction between the stellar magnetic perturbation  $\delta \vec{B}_*$ which is known at each time and given by the rotating magnetic dipole Eqs. (12)-(14) or any kind of multipolar fluctuation associated with the star and those caused by the current in the disk  $\delta\vec{B}$ in response to the former perturbation.

Expressing Eq. (24) in cylindrical coordinates for the radial component leads to:

 

$$\rho ~ \frac{D^2\xi_r}{Dt^2} - 2 ~ \rho ~ \Omega ~ \frac{D\xi_\varphi}{Dt} - \frac{\partial\Pi}{\partial r}$$

$$+ \frac{1}{\rho} ~ \frac{{\rm d}p}{{\rm d}r} ~ {\rm div}~(\rho~\vec{\xi}) + \rho ~ r ~ \vec{\xi} \cdot \vec{\nabla}~(\Omega^2)$$

$$- \frac{1}{\mu_0} ~ \left[ B_r ~ \frac{\partial Q_r}{\partial r} ... ...rac{\partial Q_r}{\partial \varphi} - 2 ~ \frac{B_\varphi~Q_\varphi}{r} \right.$$

$$+ \left. B_z ~ \frac{\partial Q_r}{\partial z} + Q_r ~ \frac{\par... ...artial B_r}{\partial \varphi} + Q_z ~ \frac{\partial B_r}{\partial z} \right] =$$

$$\frac{1}{\mu_0} ~ \left[ \delta B_*^z ~ \left\{ \frac{\partial}{\... ...{\partial}{\partial r} \left( B_z + Q_z + \delta B_*^z \right) \right\} \right.$$

$$- \left. \frac{\delta B_*^\varphi}{r} ~ \left\{ \frac{\partial}{\... ... \frac{\partial}{\partial\varphi} ( B_r + Q_r + \delta B_*^r ) \right\} \right.$$

$$+ \left. (B_z + Q_z) ~ \left( \frac{\partial\delta B_*^r}{\partial z} - \frac{\partial\delta B_*^z}{\partial r} \right) \right.$$

$$- \left. \frac{B_\varphi + Q_\varphi}{r} ~ \left( \frac{\partial}... ...lta B_*^\varphi) - \frac{\partial\delta B_*^r}{\partial\varphi} \right) \right]$$ (25)

while for the azimuthal component we get:

 

$$\rho ~ \frac{D^2\xi_\varphi}{Dt^2} + 2 ~ \rho ~ \Omega ~ \frac{D\xi_r}{Dt} - \frac{1}{r} ~ \frac{\partial\Pi}{\partial\varphi}$$

$$- \frac{1}{\mu_0} ~ \left[ B_r ~ \frac{\partial Q_\varphi}{\parti... ...~ \frac{\partial Q_\varphi}{\partial \varphi} + \frac{B_\varphi~Q_r}{r} \right.$$

$$+ \left. B_z ~ \frac{\partial Q_\varphi}{\partial z} + Q_r ~ \fra... ...\frac{Q_\varphi~B_r}{r} + Q_z ~ \frac{\partial B_\varphi}{\partial z} \right] =$$

$$\frac{1}{\mu_0} ~ \left[ \frac{\delta B_*^r}{r} ~ \left\{ \frac{\... ...eft( r ~ ( B_\varphi + Q_\varphi + \delta B_*^\varphi ) \right) \right. \right.$$

$$- \left.\left. \frac{\partial}{\partial\varphi} ( B_r + Q_r + \delta B_*^r ) \right\} \right.$$

$$\left. \delta B_*^z \left\{ \frac{1}{r} \frac{\partial}{\partial\... ...ial}{\partial z} ( B_\varphi + Q_\varphi + \delta B_*^\varphi )\right\} \right.$$ -

$$\times \left.\frac{B_r + Q_r}{r} ~ \left( \frac{\partial}{\partia... ...lta B_*^\varphi) - \frac{\partial\delta B_*^r}{\partial\varphi} \right) \right.$$

$$- \left. (B_z + Q_z) ~ \left( \frac{1}{r} ~ \frac{\partial\delta ... ...partial\varphi} - \frac{\partial\delta B_*^\varphi}{\partial z} \right) \right]$$ (26)

and finally for the vertical oscillations we have:

 

$$\rho ~ \frac{D^2\xi_z}{Dt^2} - \frac{\partial\Pi}{\partial z} - \... ...\partial Q_z}{\partial \varphi} + B_z ~ \frac{\partial Q_z}{\partial z} \right.$$

$$+ \left. Q_r ~ \frac{\partial B_z}{\partial r} + \frac{Q_\varphi}... ... \frac{\partial B_z}{\partial z} \right] = - {\rm div}~(\rho~\vec{\xi}) ) ~ g_z$$

$$+ \frac{1}{\mu_0} ~ \left[ \delta B_*^\varphi \left\{ \frac{1}{r} \frac{\partial}{\partial\varphi} ( B_z + Q_z + \delta B_*^z ) \right. \right.$$

$$- \left. \left. \frac{\partial}{\partial z} ( B_\varphi + Q_\varphi + \delta B_*^\varphi )\right\} \right.$$

$$- \delta B_*^r ~ \left\{ \frac{\partial}{\partial z} \left( B_r +... ... - \frac{\partial}{\partial r} \left( B_z + Q_z + \delta B_*^z \right) \right\}$$

$$+ \left. ( B_\varphi + Q_\varphi ) \left( \frac{1}{r} ~ \frac{\pa... ...partial\varphi} - \frac{\partial\delta B_*^\varphi}{\partial z} \right) \right.$$

$$- \left. (B_r + Q_r) \left( \frac{\partial\delta B_*^r}{\partial z} - \frac{\partial\delta B_*^z}{\partial r} \right) \right] \cdot$$ (27)

We emphasize the fact that these expressions contain terms of the form  $\xi_i ~ \delta B_*^j$, $(i,j=r,\varphi,z)$ that are second order with respect to the perturbation and that, therefore, could be neglected as long as $\delta \vec{B}_*$ remains weak. But in doing so, we suppress the parametric resonance to be studied in more detail below. This resonance becomes relevant when the perturbed stellar magnetic field is on the same order of magnitude as the aligned component. Depending on the magnitude of the perturbation, which we suppose small compared to the background field, this resonance will develop on a timescale closely related to the amplitude of the perturbation so should not be ignored. The terms  $\delta B_*^i~\delta B_*^j$ are also second order when the magnetic asymmetry is small.

Finding an analytical stability criterion for this system is a complicated or perhaps even an impossible task. Furthermore, we cannot apply the classical expansion in plane wave solutions leading to an eigenvalue problem. The presence of some coefficients varying periodically in time, like  $\delta \vec{B}_*$, prevents such a treatment. However, the problem can be cast in a more convenient form if we decouple the oscillations in the orbital plane from those perpendicular to it, as is done in the following subsection.

   
3.2 Detailed linear analysis of the thin disk

If we assume that there is no warping in the disk and neglect the gradient and displacement in the vertical direction, we can perform a detailed and accurate analysis of the motions in the disk. Indeed, by setting $\xi_z=0$ and  $\partial/\partial z =0$ in Eqs. (25) and (26), the linearized MHD equations become much more tractable. Moreover, we first assume a simple geometric structure in which the disk is replaced by a column of plasma threaded by a vertical magnetic field and then allow only magnetic perturbation in this direction. This situation seems far from the real system of an oblique rotator; nevertheless, it gives a good first insight into the oscillations and resonances arising there.

As a result, setting  $B_r=B_\varphi=\delta B_*^r = \delta B_*^\varphi = 0$ we find for the magnetic perturbation

                                   Q_r = 0 (28)

$$Q_\varphi$$ = 0 (29)

$$- \frac{1}{r} ~ \left[ \frac{\partial}{\partial r} \left( r ~ B_z... ...t) + \frac{\partial}{\partial \varphi} \left( B_z ~ \xi_\varphi \right) \right]$$ Q_z = (30)

and for the opposite of the total pressure perturbation:

$$\Pi = \frac{\gamma ~ p}{r} ~ \left[ \frac{\partial}{\partial ... ...r ~ \frac{\partial p}{\partial r} - \frac{B_z~Q_z}{\mu_0}\cdot$$ (31)

We seek solutions by writing each perturbation, such as the components of the perturbed magnetic field  $\delta\vec{B}$ (i.e. those of $\vec{Q}$), those of the Lagrangian displacement $\vec{\xi}$, those of the perturbed velocity  $\delta\vec{v}$, the perturbed density  $\delta \rho$, and the perturbed pressure $\delta p$ as

$$X(r,\varphi,t) = {\rm Re} [\tilde{X}(r) ~ {\rm e}^{{\rm i}(m~\varphi-\sigma~t)}],$$ (32)

where m is the azimuthal wavenumber and $\sigma$ the eigenfrequency of the perturbation related to the speed pattern  $\Omega_{\rm p}$ by $\sigma = m~\Omega_{\rm p}$. Therefore

  

$$\tilde{Q}_z - \frac{1}{r} ~ \frac{\partial}{\partial r} \left( r ~ B_z ~ \tilde{\xi}_r \right) - \frac{{\rm i}~m}{r} ~ B_z ~ \tilde{\xi}_\varphi$$ = (33)

$$\tilde{\Pi} \rho ~ c_{{\rm ma}z}^2 \left[ \frac{1}{r} ~ \frac{\partial}{\part... ...e{\xi}_r ~ \frac{\partial}{\partial r} \left( p + \frac{B_z^2}{2~\mu_0} \right)$$ = (34)

$c_{{\rm ma}z}^2= c_{\rm s}^2 + c_{{\rm a}~z}^2$ is the vertical fast magneto-acoustic wave speed. The sound speed and the Alfven velocity associated with the vertical magnetic field are respectively:

$$c_{\rm s}^2 \frac{\gamma~p}{\rho}$$ (35)

$$c_{{\rm a}~z}^2 \frac{B_z^2}{\mu_0~\rho}\cdot$$ (36)

Neglecting quadratic expressions in the excitation term, like  $Q_z\delta B_*^z$, the azimuthal Lagrangian displacement Eq. (26) is solved by:

 

$$\tilde{\xi}_\varphi = - \frac{\rm i}{\omega_*^2} ~ \left[ - \frac... ...2}{r^2} ~ \frac{\partial}{\partial r} \left( r ~ \tilde{\xi}_r \right) \right].$$ (37)

Due to the differential rotation of the disk, the eigenfrequency $\sigma$ appears Doppler-shifted. It is therefore convenient to introduce the frequency $\omega$ and a second one $\omega_*$ taking the speed of the magneto-sonic wave into account by

$$\omega = \sigma - m ~ \Omega$$ (38)

$$\omega_*^2 = \omega^2 - \frac{m^2 ~ c_{{\rm ma}z}^2}{r^2}\cdot$$ (39)

Strictly speaking, corotation is achieved when  $\omega_*=0$, such that there exist two solutions. When far from the corotation points, because of the weakly magnetized thin accretion disk approximation  $O(c_{\rm s}^2) = O(c_{{\rm a}~z}^2) \ll O(r^2~\omega^2)$ and therefore  $\omega_*^2 \approx \omega^2$. Substituting the expressions for  $\tilde{Q}_z$, Eq. (33), for  $\tilde{\Pi}$, Eq. (34) and for  $\tilde{\xi}_\varphi$, Eq. (37), in the radial equation Eq. (25), we obtain to lowest order:

 

$$c_{{\rm ma}z}^2 ~ \left[ \frac{\partial^2\tilde{\xi}_r}{\partial ... ...{\xi}_r = \frac{B_z}{\mu_0~\rho} ~ \frac{\partial\delta B_*^z}{\partial r}\cdot$$ (40)

The radial epicyclic frequency of a single particle as measured by an observer in the rest frame is

$$\kappa_r^2 = 4~\Omega^2 + r~\frac{\partial\Omega^2}{\partial r}$$ (41)

(see Appendix A). The solutions of this equation consist of free waves corresponding to density perturbations propagating in the disk without magnetic asymmetry related to the homogeneous part and of non-wavelike disturbances due to this asymmetric component related to the inhomogeneous part.

3.2.1 Free wave solutions

When looking for free wave solutions to Eq. (40), we can set  $\delta B_*^z=0$. Thus

$$c_{{\rm ma}z}^2 ~ \left[ \frac{\partial^2\tilde{\xi}_r}{\pa... ... r} \right] + ( \omega^2 - \kappa_r^2) ~ \tilde{\xi}_r = 0 .$$ (42)

This is the generalization to MHD of the sound wave that propagates in the hydrodynamical disk. We just need to replace the sound speed  $c_{\rm s}^2$ by the fast magneto-acoustic speed  $c_{{\rm ma}z}^2$. Introducing a new unknown  $\psi = \tilde{\xi}_r ~ \sqrt{r~\rho~c_{{\rm ma}z}^2}$, it satisfies a kind of Schrödinger equation

$$\psi''(r) + V(r) ~ \psi(r) = 0$$ (43)

where the potential is given by:

$$V(r) = \frac{\omega^2 - \kappa_r^2}{c_{{\rm ma}z}^2}\cdot$$ (44)

A first guess for free wave solutions is given by the WKB expansion as follows:

$$\Psi(r) = \Phi(r) ~ {\rm e}^{{\rm i}\int^r k(s) ~ {\rm d}s}.$$ (45)

Putting this approximation into Eq. (43), the dispersion relation for the fast magneto-acoustic waves propagating in the orbital plane of the accretion disk is given by:

$$\omega^2 = \kappa_r^2 + c_{{\rm ma}z}^2 ~ k^2.$$ (46)

Free waves can only propagate in regions where  $\omega^2 - \kappa_r^2=V(r)~ c_{{\rm ma}z}^2 \ge0$. The boundary between the propagating and damping zones is defined by the inner and outer Lindblad radius  $r_{\rm L}^{\rm in/out}$ defined by  $V(r_{\rm L}^{\rm in/out})=0$. Using the results of Appendix B, we can find a better approximation to the solution of Eq. (43) which is valid even for  $r\approx r_{\rm L}^{\rm in/out}$. For the inner Lindblad resonance of interest here, we introduce the following function $\omega_1$, writing  $r_{\rm L}=r_{\rm L}^{\rm in}$:

$$\omega_1(r) - \left[ - \frac{3}{2} ~ \int_{r_{\rm L}}^r \sqrt{V(s)} ~ {\rm d}s \right]^{2/3} {\rm for\;} r \le r_{\rm L}$$ = (47)

$$\omega_1(r) \left[ \frac{3}{2} ~ \int_{r_{\rm L}}^r \sqrt{-V(s)} ~ {\rm d}s \right]^{2/3} {\rm for\;} r \ge r_{\rm L}.$$ = (48)

The function $\psi$ is then a linear combination of the 2 linearly independent solutions:

 

$$\psi_1(r) \frac{A{\rm i}(\omega_1(r))}{\sqrt{\vert\omega_1'(r)\vert}}$$ = (49)

$$\psi_2(r) \frac{B{\rm i}(\omega_1(r))}{\sqrt{\vert\omega_1'(r)\vert}}\cdot$$ = (50)

Furthermore, we require the solution to remain bounded, which leads to C₂=0. Thus the solution for the Lagrangian displacement is $\xi_r=C_1~\psi_1(r)/\sqrt{r~\rho~c_{{\rm ma}z}^2 }$. At the inner boundary of the accretion disk, in accordance with the simulations performed later on, the density perturbation should vanish. This is expressed as  $\delta\rho=0$. To the lowest order consistent with our approximation, the Lagrangian radial displacement $\xi_r$ must satisfy

$$\tilde{\xi}_r'(R_1) + \left( \left.\frac{\partial\ln(r~\rho... ...~m~\frac{\Omega}{\omega~R_1} \right) ~ \tilde{\xi}_r(R_1) = 0.$$ (51)

This last condition determines the eigenfrequencies $\sigma$ as a function of the azimuthal mode m. For any m, there is an infinite set of eigenvalues. However, the corresponding eigenfunctions become more and more oscillatory, implying larger and larger wave numbers. In the numerical applications, we restrict our attention to the first few that also correspond to the highest possible eigenvalues $\sigma$.

3.2.2 Newtonian vs. general-relativistic disk

We emphasize that the qualitative results presented in this paper are independent of the nature of the spacetime, be it flat as in the Newtonian gravitational potential or be it curved as in the Kerr metric. We make this point clear by showing results according to both spacetime structures.

Let's start by recalling some facts about accretion disks orbiting in the Kerr metric. When the inner edge of the accretion disk reaches a few gravitational radii, general relativistic effects become important. The degeneracy between the three frequencies is lifted, namely the orbital $\Omega$, the radial epicyclic $\kappa_r$, and the vertical epicyclic $\kappa_z$ frequencies. Their value depends not only on the stellar mass M_* but also on the angular momentum of the star a_*. Indeed we distinguish 3 characteristic frequencies in the accretion disk around a Kerr black hole (or equivalently a rotating neutron star):

• the orbital angular velocity:
  

  $$\Omega(r,a_*)=\frac{1}{r^{3/2}+a_*} ;$$ (52)

  
  

• the radial epicyclic frequency:
  

  $$\kappa_r(r,a_*)=\Omega(r,a_*)~\sqrt{1-\frac{6}{r} + 8~\frac{a_*}{r^{3/2}} - 3~\frac{a_*^2}{r^2} } ;$$ (53)

  
  

• the vertical epicyclic frequency:
  

  $$\kappa_z(r,a_*)=\Omega(r,a_*)~\sqrt{1 -4~\frac{a_*}{r^{3/2}} + 3~\frac{a_*^2}{r^2}} \cdot$$ (54)

  
  

The parameter a_* corresponds to the angular momentum of the star in geometrized units. For a neutron star, it is given by $a_*=\frac{c~I_*}{G~M_*^2}~\Omega_*$.

We have the following ordering:

$$\Omega>\kappa_z>\kappa_r {\rm for} \;\; a_*>0$$ (55)

$$\kappa_z>\Omega>\kappa_r {\rm for} \;\; a_*<0.$$ (56)

3.2.3 Results

The eigenvalues for the density waves are shown with decreasing magnitude in Table 1. This holds for a neutron star with angular velocity $\nu_*=\Omega_*/2\pi=100~{\rm Hz}$. We compared the Newtonian case with the Schwarzschild metric. The highest speed pattern, given by $\sigma/m$, never exceeds the orbital frequency at the ISCO. Note that in the m=1 case, the inner Lindblad resonance does not exist for the Newtonian potential, and only a few free wave solutions can be found for the Schwarzschild geometry.

Table 1: The first eight largest eigenvalues $\sigma$ for the free wave solutions of Eq. (43). Values are normalized to the frequency of the ISCO,  $\Omega _{\rm isco}=6^{-3/2}$. Results are given for 3 azimuthal modes m=1,2,5 for the Newtonian, as well as for the Schwarzschild gravitational field. Where no inner Lindblad resonance exist, this is represented by -.

Some examples of the corresponding eigenfunctions for density waves are shown in Fig. 1 with arbitrary normalization. Each of them possesses its own inner Lindblad radius depending on the eigenvalue. They are distinguished by the number of radial nodes they possess, starting from no node at all, which corresponds to the highest speed pattern.

	Figure 1: Density wave perturbation in the disk caused by the free wave propagation for the azimuthal mode m=2. Some examples are shown for different eigenvalues and for Newtonian geometry, on the left, as well as for Schwarzschild geometry, on the right. The vertical bar indicates the location of the inner Lindblad resonance. The normalization of the eigenfunctions is arbitrary.
Open with DEXTER

Because the pattern speed of the first spiral waves is close to the ISCO, we suggest that these free wave solutions can be associated with the highest kHz QPO. More precisely, there exists one eigenfunction which possesses no root, namely the one with the highest eigenvalue. Eigenfunctions with several nodes in the radial distribution will emit no significant radiation since on average the fluctuations will cancel out. Therefore, we expect to see only density fluctuations with no or very few radial nodes. The most interesting candidate is the one with no node, a situation that also corresponds to the highest eigenvalue.

Moreover, while propagating in the perturbed rotating gravitational field of the star, this spiral wave experiences a sinusoidally varying gravity at a rate  $\Omega_{\rm p}{-}\Omega_*$, where $\Omega_{\rm p}$ is the speed pattern of the wave. This induces a modulation in the shape of the wave which will be reflected in the power spectrum density as a kind of beat phenomenon. As a result, in our model, we associate the highest kHz-QPO to $\nu_2=\Omega_{\rm p}$ and the lowest kHz QPO to the beat frequency  $\nu_1=\Omega_{\rm p}-\Omega_*$. The peak separation is then  $\Delta\nu=\Omega_*$.

In a real accretion disk, the precise location of the inner edge does not necessarily reach the ISCO. The magnetic field or the radiation pressure can prevent the disk from doing so. Indeed, the boundary of the inner part of the disk can fluctuate due to the varying accretion rate. For instance when the accretion process is enhanced, the inner edge moves closer to the ISCO. As a result, the highest eigenvalue of the free waves also increases and therefore $\nu_2$ shifts to higher frequencies as well as $\nu_1$.

Examples of this mechanism are shown in Fig. 2. We estimate the highest eigenvalue for the modes m=2,5 in the Newtonian, as well as in the Schwarzschild gravitational potential. The pattern speed of the wave solution  $\Omega_{\rm p}$, normalized to the frequency of the ISCO,  $\Omega_{\rm isco}$, is shown as a function of the location of the inner boundary of the accretion disk  $R_{\rm in}$. In each case, we observe an increase in the eigenvalue when the disk's inner edge moves closer to the compact object, corresponding to an increased accretion rate. When the ISCO is reached, the boundary condition does not change anymore and the eigenvalue saturates to its final value. Thus the kHz-QPO frequencies associated with the speed pattern of the wave saturate. This has been observed in some LMXBs, as reported for instance in a paper by Zhang et al. (1998).

Figure 2: Variation of the highest eigenvalue, corresponding to the eigenfunction having no node, as a function of the location of the inner edge of the disk  $R_{\rm in}$. There is a monotonic increase as the disk approaches the ISCO at  $R_{\rm in}=6R_{\rm g}$. Results are shown for the m=2,5 modes in the Newtonian (N) and Schwarzschild (S) spacetime, red and blue curves respectively. The gravitational radius is defined by  $R_{\rm g}=G~M_*/c^2$.

Open with DEXTER

   
3.2.4 Non-wavelike perturbation

We now solve the full inhomogeneous Eq. (40) in order to find the solution corresponding to the non-wavelike perturbation. In this case, the eigenvalue is known:  $\sigma = m~\Omega_*$. We have to solve the second order ordinary differential equation for $\psi$ with the appropriate boundary conditions Eq. (51). Following the results given in Appendix B, the most general solution is:

 

$$\psi_r(r) = C_1 ~ \psi_1(r) + C_2 ~ \psi_2(r) + \pi ~ {\rm sign}(... ...r \left( \psi_1(r) ~ \psi_2(s) - \psi_1(s) ~ \psi_2(r) \right) F(s) ~ {\rm d}s.$$ (57)

The constant C₂ is chosen such that the solution remains bounded for  $r\gg r_{\rm L}$:

$$C_2 = \lim_{r\to\infty} \pi ~ {\rm sign}(\omega_1') ~ \int_{r_{\rm L}}^r \psi_1(s) ~ F(s) ~ {\rm d}s.$$ (58)

This integral is convergent because the function $\psi_1$ is exponentially decreasing with the radius. We can directly integrate this equation to obtain the shape of the density perturbation between the inner edge of the accretion disk and the inner Lindblad resonance. The results for a particular choice of the density and magnetic profiles (discussed in Sect. 5) are shown in Fig. 3. This can be used later to check the numerical algorithm. It also gives a good estimate of the amplitude of these fluctuations because it is directly related to the strength of the stellar magnetic field asymmetry. The maximum density fluctuation is about 0.1%. The oscillations become evanescent upon reaching the inner Lindblad resonance.

	Figure 3: Non-wavelike disturbance in the inner part of the accretion disk due to the asymmetric component of the magnetic field. The potential is Newtonian and m=2. The numerical values used are described in Sect. 5.
Open with DEXTER

   
3.2.5 Counterrotating disk

The counterrotating disk can be solved in exactly the same way by integrating Eq. (40) to find the solution corresponding to the non-wavelike perturbation. The amplitude of the oscillations is shown in Fig. 4. The Lindblad resonances are no longer in the computational domain. An oscillatory structure appears in the whole disk. The maximum amplitude is about 1% of the average density. This can later be compared with the numerical simulations.

	Figure 4: Non-wavelike disturbance in the inner part of the accretion disk due to the asymmetric component of the magnetic field. The potential is Newtonian, m=2 and the disk is counterrotating.
Open with DEXTER

4 Simplified analysis

To get more insight into the nature of the resonances, we now focus on the displacement of the disk either in the equatorial plane or perpendicular to it. This means that we neglect the coupling between the oscillations occurring in perpendicular directions. This simplification will help us to understand the origin of the instability and to derive some resonance conditions.

To study the motion in the vertical direction, setting  $(\xi_r, \xi_\varphi) = 0$ we find:

 

$$\rho ~ \frac{D^2\xi_z}{Dt^2} - \frac{\partial}{\partial z} \left[... ... \frac{\delta B_*^r }{\mu_0} ~ \frac{\partial B_r}{\partial z} ~ \xi_z \right).$$ (59)

See Appendix A for the details. The Alfven velocity associated with the radial magnetic field and the radial fast magneto-acoustic wave speed are, respectively:

$$c_{{\rm a}r}^2 = \frac{B_r^2}{\mu_0~\rho}$$ (60)

$$c_{{\rm ma}r}^2 = c_s^2 + c_{{\rm a}r}^2.$$ (61)

The same can be done for the orbital motion by setting $\xi_z=0$, which leads to a similar expression:

 

$$\rho ~ \frac{D^2\xi_r}{Dt^2} - \frac{\partial}{\partial r} \left[... ...( \frac{\delta B_*^z}{\mu_0} ~ \frac{\partial B_z}{\partial r} ~ \xi_r \right).$$ (62)

Equations (59) and (62) look very similar, the discrepancy coming only from the difference between the planar and the cylindrical geometry (terms containing r). In the two first terms of Eqs. (59) and (62) we recognize an MHD wave propagation in a tube of spatially varying but time independent cross section (Morse & Feshbach 1953). The first and third terms put together form an harmonic oscillator at the epicyclic frequency. So the three first terms are a generalization of the Klein-Gordon equation and do not give rise to any kind of instability. The interesting parts are those containing the perturbation in the magnetic field  $\delta \vec{B}_*$.

Further, we neglect the magneto-acoustic wave propagation because of the weakly magnetized thin disk approximation, and also assume a homogeneous medium for which the vertical Lagrangian displacement and density are independent of the altitude z in the first case and where the density and radial Lagrangian displacement do not depend on the radius r in the second case. Note also that terms like  ${\delta B_*^i}^2$ are second order with respect to the perturbation and should be neglected. These simplifications make the equations much more easier to manipulate without losing the most interesting physical aspect contained in the resonance terms. We can then simplified Eqs. (59) and (62) to obtain:

 

$$\frac{D^2\xi_z}{Dt^2} + \left[ \kappa_z^2 - \frac{\partial}{\part... ... \frac{\partial}{\partial z} \left( \frac{B_r~\delta B_*^r}{\mu_0~\rho} \right)$$ (63)

$$\frac{D^2\xi_r}{Dt^2} + \left[ \kappa_r^2 - \frac{\partial}{\part... ...c{\partial}{\partial r} \left( \frac{B_z~\delta B_*^z}{\mu_0~\rho} \right)\cdot$$ (64)

Recall that the magnetic perturbation is dragged with the star's rotation. Therefore, when looking at a given fixed point in the accretion disk, the coefficients of the equations containing  $\delta B_*^i$ will vary periodically. Thus, we recognize a kind of Hill equation corresponding to an oscillator with periodically time-varying eigenfrequency and subject to a periodic driven force. It is well known that this type of equation shows what is called a parametric resonance. Moreover, due to the rotation of the star, this perturbation will vary sinusoidally in time. In fact, in the frame locally corotating with the disk, (convective derivative), the modulation follows a harmonic temporal behavior such as:

$$\cos(\varphi - (\Omega-\Omega_*)~t) \quad\sin(\varphi - (\Omega-\Omega_*)~t).$$ , (65)

The Hill equation specializes to Mathieu's equation for which we can derive the resonance conditions analytically. The solutions to Mathieu's equation written in the form  $y''(t) + \omega_0^2 ( 1 + h ~ \cos \gamma~t) ~ y(t) = 0$ will indeed become unbounded if  $\gamma = 2 \frac{\omega_0}{n}$ where n is an integer, (Landau & Lifshitz 1982).

The equation satisfied by each component of the Lagrangian displacement when looking at a given azimuth $\varphi_0$ is therefore in the form:

 

$$\xi''(t) + \left[ \kappa^2 + h ~ \left\{ \begin{array}{c} \sin \\... ...{c} \sin \\ \cos \end{array}\right\} (\varphi_0 - ( \Omega - \Omega_* ) ~ t ),$$ (66)

where $f(\delta B_*^i)$ and  $h(\delta B_*^i)$ are two constants depending on the strength of the stellar magnetic perturbation.

From this analysis we expect three kinds of resonances corresponding to, (see Appendix C for a detailed analysis):

• a corotation resonance at the radius where the angular velocity of the disk equals the rotation speed of the star. This is only possible for prograde disks. The resonance condition to determine the corotating radius is simply
  

  $$\Omega = \Omega_* ;$$ (67)

  
  

• an inner and an outer Lindblad resonance at the radii where the radial or vertical epicyclic frequency equals the frequency of the dipolar magnetic perturbation as measured in the locally comoving frame. Due to the m=1 structure of this perturbation we find the resonance condition to be
  

  $$\vert\Omega-\Omega_*\vert = \kappa_{r/z};$$ (68)

  
  

• a parametric resonance related to the time-varying epicyclic frequency, (Hill equation). The rotation of the star induces a sinusoidally time varying epicyclic frequency leading to the well-known Mathieu equation. The resonance condition is then
  

  $$\vert\Omega-\Omega_*\vert = 2 ~ \frac{\kappa_{r/z}}{n}\cdot$$ (69)

  
  

Note that the driven resonance occurs at the same radius as does the parametric resonance for n=2. However, their growth rates are different, because forcing causes a linear growth in time of the amplitude while parametric resonance causes an exponential growth.

More generally, for a perturbation of arbitrary azimuthal mode m, the resonance conditions are given by:

• for the inner and outer Lindblad resonances: $m~\vert \Omega - \Omega_* \vert = \kappa_{r/z}$
• for the parametric resonance: $m~\vert\Omega - \Omega_*\vert = 2~\frac{\kappa_{r/z}}{n}$.

The conclusion of this study is that the locations of the resonances are determined solely by the rotation rate of the neutron star. Fluctuations of the accretion rate are not included in this picture. However, we will show by a more detailed linear analysis that there are some free-wave solutions for the density perturbation for which the eigenvalue can be correlated to the accretion rate in accordance with observations, see Sect. 3.2.

4.1 Results

4.1.1 Newtonian disk

From the resonance conditions derived previously, Eqs. (67)-(69), we can estimate the radii where each of these resonances occurs. Beginning with the Newtonian potential, it is well known that the angular velocity, the radial and epicyclic frequencies are all equal. For a weakly magnetized thin accretion disk, this statement remains true so that  $\Omega \approx \kappa_r \approx \kappa_z$. Distinguishing between the two signs of the absolute value, we get the following rotation rate for the parametric condition (69):

$$\frac{\Omega}{\Omega_*} = \frac{n}{n+2} = 1/3, 1/2, 3/5, 2/3, 5/7, ...$$ (70)

$$\frac{\Omega}{\Omega_*} = \frac{n}{n-2} = -1 , - , 3, 2, 5/3, ...$$ (71)

As a consequence, each resonance will occur in the frequency range  $[\Omega_*/3, 3~\Omega_*]$. Assuming that the QPOs are related to the orbital motion and taking into account that the spatial structures of the resonances are formed by m=1 patterns, we expect to see QPOs in the same range  $[\nu_*/3, 3~\nu_*]~{\rm Hz}$ where we have introduced the spin frequency by  $\nu_* = \Omega_*/2\pi$. For instance, for a  $300~{\rm Hz}$ spinning neutron star, we can expect QPOs between  $[100,900]~{\rm Hz}$. However, observations have shown that QPOs with frequency up to 1300 Hz exist for this kind of system. We will see in the next subsection that this difficulty can be solved by introducing the general relativistic effects of the innermost stable circular orbit. Table 2 gives an example of the expected orbital frequencies where the resonances are expected to occur for a neutron star spinning at 300 and 600 Hz.

Table 2: Value of the orbital frequencies at the parametric resonance radii for the first five orders n in the case of a Newtonian gravitational potential. The results are given for a  $1.4~{M_\odot}$ neutron star rotating at 300 and  $600~{\rm Hz}$ respectively. The value to the left of the symbol / corresponds to the absolute value sign taken to be - and on the right to be +. The -- sign indicates that no resonance condition can be satisfied.

4.1.2 General relativistic disk

The parametric resonance condition Eq. (69), split into the two cases depending on the sign of the absolute value, becomes:

  

$$\Omega(r,a_*) - 2 ~ \frac{\kappa_{r/z}(r,a_*)}{n} = \Omega_*$$ (72)

$$\Omega(r,a_*) + 2 ~ \frac{\kappa_{r/z}(r,a_*)}{n} = \Omega_*.$$ (73)

For a typical neutron star, we choose:

• mass  $M_*=1.4~{M_\odot}$;
• angular velocity  $\nu_*=\Omega_*/2\pi=300{-}600~{\rm Hz}$;
• moment of inertia  $I_*=10^{38}\;{\rm kg~m^2}$;
• angular momentum  $a_*=\frac{c~I_*}{G~M_*^2}~\Omega_*$.

The angular momentum is then given by  $a_*=5.79*10^{-5}~\Omega_*$. For a given angular momentum a_*, we have to solve Eqs. (72) and (73). Solving them for the radius r and then deducing the orbital frequency at this radius, we get the results shown in Table 3. For the spin rate of the star we find  a_*=0.109-0.218 and so the vertical epicyclic frequency is close to the orbital one  $\kappa_z\approx\Omega$. Thus for the vertical resonance, we are still close to the Newtonian case explored in the previous section.

Table 3: Value of the orbital frequencies at the parametric resonance radii for the first five orders n in the general relativistic case. The results are given for a  $1.4~{M_\odot}$ neutron star rotating at 300 and  $600~{\rm Hz}$ respectively. The value to the left of the symbol / corresponds to the absolute value sign taken to be - and on the right to be +.

We have now overcome the problem faced in the Newtonian approximation. Indeed, even in the case of a 300 Hz spinning neutron star, we can expect frequencies as high as 1332 Hz for the n=1 rank, in accordance with observations. The general-relativistic effect caused by the presence of the ISCO is therefore a crucial characteristic that makes a clear distinction with the Newtonian potential. The ISCO has to be implemented in any realistic QPO model in order to get a good quantitative agreement with observations. Therefore, a fully general-relativistic simulation code should be developed in order to take this ISCO into account in a self-consistent manner. However, we point out that general relativistic effects are not directly responsible for the observed QPOs as claimed in other models (Lee et al. 2004; Stella et al. 1999). They merely lead to quantitative improvements compared to the flat spacetime results.

   
5 Numerical simulations

The linear analysis described above is not the end of the story. It was shown that the disk becomes unstable at some preferred radii where the resonance conditions are satisfied. But this analysis does not tell us what will happen to these resonances on a long time scale. Moreover, as in the hydrodynamical case, the non-linearities in the time evolution of the disk will give rise to a rich variety of physical processes. However, in this simulation, in order to check our numerical algorithm, we start with runs for which the perturbation remains weak, in order to keep a regime close to the linear one described in the previous section. The full non-linear effects arising from a strong magnetic field perturbation will be treated in a forthcoming paper.

To achieve this goal, we performed 2D numerical simulations by solving the magnetohydrodynamical equations for the accretion disk Eqs. (15)-(19). This is done by a pseudo-spectral method for a simplified version of the 2D MHD accretion disk. Indeed, an exact treatment of the magnetic field induced by the motions in the disk would involve solution of an elliptic partial differential equation related to Ampere's law. We don't want to go into such refinement in this paper but only emphasize the role of a rotating asymmetric magnetic field. Nevertheless the physics at hand remains the same.

5.1 Boundary and initial conditions

We expand every quantity in Fourier series in the azimuthal direction and in Tchebyshev series in the radial direction, whereas time integration is performed by an explicit fourth-order Runge-Kutta scheme. In this work we are not interested in the accretion process itself but focus only on the oscillations occurring in this system around its stationary state. We can then impose fixed boundary conditions in the radial direction by setting the radial velocity equal to zero at the two boundaries. Furthermore, to avoid spurious reflections at the edges of the accretion disk, we force the radial velocity to decrease to zero in a narrow layer around these boundaries. We also require the density and the azimuthal velocity to reach their equilibrium. This acts as an ad hoc absorption process damping out the strong oscillations which could arise in these regions. The non-reflecting boundaries conditions are therefore expressed as:

 

$$\rho(r,\varphi) = f(r) ~ ( \tilde{\rho}(r,\varphi) - \rho^{\rm eq}(r,\varphi) ) + \rho^{\rm eq}(r,\varphi)$$ (74)

$$v_r(r,\varphi) = f(r) ~ \tilde{v}_r(r,\varphi)$$ (75)

$$v_\varphi(r,\varphi) = f(r) ~ ( \tilde{v}_\varphi(r,\varphi) - v_\varphi^{\rm eq}(r,\varphi) ) + v_\varphi^{\rm eq}(r,\varphi)$$ (76)

$$B_z(r,\varphi) = f(r) ~ ( \tilde{B}_z(r,\varphi) - B_z^{\rm eq}(r,\varphi) ) + B_z^{\rm eq}(r,\varphi).$$ (77)

The quantities denoted by an $\tilde{X}$ subscript correspond to the physical values before applying the non-reflecting boundary conditions. Those values at equilibrium are denoted by  $X^{\rm eq}$. Function f acts as a filter, as its value is one everywhere except close to the inner and outer boundaries where it decreases to zero in a narrow transition layer.

We adopted the geometry of a cylinder of gas threaded by a vertical magnetic field and a radially directed gravitational field. The problem is invariant under translation in the direction of the magnetic field lines, which is chosen to be the z direction, $\partial/\partial z =0$. Motions of the gas are only allowed in the equatorial plane  $(r,\varphi)$, v_z=0; therefore, we suppress the possibility of a warped disk and a possible precession of its orbital plane, leaving this to future work. The ideal MHD Eqs. (15) and (16) then simplify to:

 

$$\frac{\partial \rho}{\partial t} + \frac{1}{r}~\frac{\partial}{\p... ...(r~\rho~v_r) + \frac{1}{r} ~ \frac{\partial}{\partial \varphi} (\rho~v_\varphi)$$ = 0 (78)

$$\frac{\partial v_r}{\partial t} + v_r ~ \frac{\partial v_r}{\part... ...ac{v_\varphi}{r} ~ \frac{\partial v_r}{\partial\varphi} - \frac{v_\varphi^2}{r} g_r - \frac{1}{\rho} ~ \frac{\partial p}{\partial r} + \frac{j_\varphi ~ B_z}{\rho}$$ = (79)

$$\frac{\partial v_\varphi}{\partial t} + v_r ~ \frac{\partial v_\v... ...hi}{r} ~ \frac{\partial v_\varphi}{\partial\varphi} + \frac{v_r ~ v_\varphi}{r} g_\varphi - \frac{1}{\rho~r} ~ \frac{\partial p}{\partial \varphi} - \frac{j_r ~ B_z}{\rho} \cdot$$ = (80)

Moreover, we adopt a polytropic equation of state for the plasma which is:

$$p = K ~ \rho^\gamma.$$ (81)

The induction Eq. (18) in this particular geometry transforms into an advection equation for the vertical component B_z, similar to the mass conservation. This means that the magnetic flux is conserved during time evolution:

$$\frac{\partial B_z}{\partial t} + \frac{1}{r}~\frac{\partial}{\pa... ...} (r~B_z~v_r) + \frac{1}{r} ~ \frac{\partial}{\partial \varphi} (B_z~v_\varphi)$$ = 0                             (82)

$$\frac{\partial B_r}{\partial t}$$ = 0 (83)

$$\frac{\partial B_\varphi}{\partial t}$$ = 0 . (84)

Inspecting these expressions, we conclude that there is no radial or azimuthal magnetic field generated during the evolution of the disk. Because at the initial conditions  $B_r=B_\varphi=0$, these components will remain equal to zero during the whole simulation. Consequently the divergence free condition of $\vec{B}$ is automatically satisfied by:

$${\rm div}~\vec{B} = \frac{\partial B_z}{\partial z} = 0.$$ (85)

Finally, applying Ampere's law, Eq. (19), we get the current density in the disk by:

$$\frac{1}{\mu_0 ~ r} ~ \frac{\partial B_z}{\partial\varphi}$$           j_r = (86)

$$j_\varphi - \frac{1}{\mu_0} ~ \frac{\partial B_z}{\partial r}$$ = (87)

j_z = 0 . (88)

We expand every quantity in a Fourier series in the azimuthal direction and in a Tchebyshev series in the radial direction. Moreover, to avoid the spectral aliasing that arises from the non-linearities of the equations, we use the 2/3-rule to completely wash out the highest frequencies (Gottlieb & Orszag 1977; Canuto et al. 1990) (desaliasing process).

We use geometrized units in which G=c=1. Distance is measured in units of the gravitational radius given by  $R_{\rm g}=G~M_*/c^2$. Moreover, the simulations are done for a star with M_*=1, so that in the new units we have  $R_{\rm g}=1$. The star is assumed to be perfectly spherical with the three main axes all equal  R_x=R_y=R_z, and the perturbation comes from the magnetic field alone. The adiabatic constant is equal to  $\gamma=5/3$ and typical resolution is  $256\times32$.

We checked that this resolution is sufficient by running a simulation with twice as many points in both directions, namely  $512\times64$. We found no changes in the solutions. Indeed inspecting the Fourier-Tchebyshev coefficients, only half of the radial coefficients and a few of the azimuthal coefficients in the 2D expansions are significantly different from zero. Choosing a higher resolution will not lead to any improvement in the solution but only to excessive computation time.

The initial density profile is prescribed as

$$\rho_0(r) = \frac{10^{-6}}{r}\cdot$$ (89)

Initially the contribution to the magnetic field is given by the star and the disk as follows

 

$$- 5 ~ \frac{10^{-3}}{r^3}$$          B_z^* = (90)

$$\frac{10^{-3}}{r^3},$$ B_z^d = (91)

while the stellar magnetic perturbation is given by

 

$$\delta B_z^* \delta b ~ B_z^* ~ \cos( m~(\varphi - \Omega_* ~ t ))$$ = (92)

$$\delta b$$ = 0.001 . (93)

The r-dependence of the magnetic perturbation Eq. (92) is chosen to mimic the rotation of a misaligned rotator rotating at a rate $\Omega _*$ with a typical decay as r^-3. All other components are equal to zero. The azimuthal mode of the magnetic field perturbation  $\delta B_z^*$ can be chosen arbitrarily. We use the case m=2 to directly compare with the hydrodynamical simulations. However, the m=1 case is better adapted to the true dipolar structure of the stellar magnetic field. We will give some examples for this situation, too.

With these prescribed initial conditions, the thin disk approximation is fairly good. Indeed the disk ratio H/R, H typical height and R radius of the disk, plotted for the Newtonian and Schwarzschild geometry in Fig. 5, never exceed 0.1 in any run. Moreover, the disk is threaded by a weak magnetic field except very close to the inner edge. The plasma $\beta$ parameter is shown in Fig. 6 and remains very high.

	Figure 5: Thin disk approximation in which the ratio H/R never exceeds 0.1, on the left for the Newtonian potential and on the right for the Schwarzschild geometry.
Open with DEXTER

	Figure 6: Plasma $\beta$ parameter for the weakly magnetized disk plotted on a logarithmic scale. The magnetic field dominates only very close to the inner edge, on the left for the Newtonian potential and on the right for the Schwarzschild geometry.
Open with DEXTER

Before the time t=0, the disk stays in its axisymmetric equilibrium state and possesses only an azimuthal motion. At t=0, we switch on the perturbation by adding an asymmetric rotating component to the magnetic field. We then let the system evolve during more than one thousand orbital revolutions of the inner edge of the disk. We performed four sets of simulations. In the first one, the gravitational potential was Newtonian. In the second, we used a pseudo-Newtonian potential in order to take the ISCO into account. This is well-suited to describing the Schwarzschild space-time. In the third one, we took the angular momentum of the star into account by introducing a pseudo-Kerr geometry. And finally in the fourth and last set, we performed simulations with a counter rotating accretion disk evolving in the Newtonian potential described in the first set.

5.2 Numerical test

We have tested the pseudo-spectral method used in this paper on some standard problems. For instance, we solved the Euler equations in cylindrical coordinates for an initial jump in the density profile and retrieved the classical results of the shock problem with the correct shock speed. However, in order to reduce the oscillations (Gibbs phenomenon) due to discontinuity at the shock, we filtered out the high frequency components in Fourier-Tchebyshev space.

In a second set of tests, we launched an m=1 free wave in the magnetized accretion disk, without any perturbation in the magnetic field. The wave is slowly damped due to the absorbing boundary conditions. The linear analysis elaborated in the previous section is in agreement with the non-linear simulations presented below, at least when the perturbations in the variables  $(\rho,r,\vec{v},\vec{B})$remain small. For problems without shock or discontinuity in the flow, this pseudo-spectral code gives very accurate solutions with only a few discretization points.

   
5.3 Newtonian potential

First, we study the behavior of the disk in the Newtonian potential. In this case, the Keplerian rotation rate, the vertical and the radial epicyclic frequencies for a single particle are all identical. However, for a thin gaseous disk, there is a slight difference on the order H/R and we have:

$$\Omega_k \approx \kappa_r \approx \kappa_z.$$ (94)

The normalized rotation rate of the star around the z-axis is  $\Omega _*=0.0043311$. Assuming a  $1.4~{\rm M_\odot}$ neutron star, this corresponds to a spin of $\nu_*=100$ Hz. The disk's inner boundary is located at R₁=6 and the outer boundary at R₂=60. The orbital angular motion at the inner edge of the disk is  $\Omega_{\rm in}=R_1^{-3/2}=0.0680$. We normalize the time by dividing it by the spin period of the star  $T_* = \frac{2~\pi}{\Omega_*} = 1450.7$.

	Figure 7: Final snapshot of the density perturbation  $\delta \rho /\rho _0$ in the accretion disk evolving in a quadrupolar perturbed Newtonian potential. The disk extends from R1=6.0 to R2=60.0. The rotation rate of the star is  $\Omega _*=0.0043311$. The time is normalized to the spin period  T*=1450.7. The m=2 structure emerges in relation to the m=2 quadrupolar magnetic perturbation.
Open with DEXTER

The final snapshot of the density perturbation in the disk calculated as  $\Delta\rho/\rho_0=\rho/\rho_0-1$ is shown in Fig. 7. We clearly recognize the spatial separation between the propagating and non-propagating regions delimited by the inner Lindblad resonance located at  $R_{\rm L}^{\rm in}=23.7$. However, the outer Lindblad resonance at  $R_{\rm L}^{\rm out}=49.3$ is not seen in this run because the magnetic perturbation close to the outer boundary is very weak compared to its amplitude close to the inner edge. The corotation resonance, located in the neighborhood of r=40.0, is also not apparent. The reason is that its growth rate is linear and very weak; an estimate is found by taking the typical driven resonance for a harmonic oscillator as

$$\frac{{\rm d}^2 x}{{\rm d}t^2} + \omega^2 ~ x = f ~ \cos (\omega~t).$$ (95)

which has a solution given by:

$$x(t) = \frac{f}{2~\omega} ~ t ~ \sin(\omega~t).$$ (96)

The numerical value from the simulation is  $f/2~\omega \approx 10^{-10}$. In order to get an appreciable amplitude at the corotation point of 10^-3 of the average density, we would have to wait a time  $\Delta t \ge 1500~T_*$, which is much more than the total duration of our run (about 70 T_*). This explains the missing corotation resonance at this stage. However, the most interesting features are the very pronounced density fluctuations in the innermost part of the accretion disk. They correspond to the non-wavelike disturbance of the density related to the inhomogeneous solution to the Schrödinger type Eq. (43). They keep their own identity even after 1000 revolutions of the innermost orbit and rotate at the stellar spin rate. After a few hundred orbital revolutions, the disk settles down to a new quasi-stationary state in which these disturbances persist. This is confirmed by inspection of Fig. 8 in which we have plotted a cross section of the density perturbation  $\delta \rho /\rho _0$ for a given azimuthal angle, namely $\varphi=0$. As expected, the density perturbations vanish at the disk edges. They behave as predicted by the linear analysis in Fig. 3 when allowance is made for the different boundary conditions.

	Figure 8: Cross section of the density perturbation  $\delta \rho /\rho _0$ in the Newtonian disk at the final time of the simulation. The inner and outer Lindblad resonances are shown by vertical lines at  $r_{\rm L}^{\rm in/out}=23.7/49.3$. Nevertheless the disk feels only the inner resonance, the outer one being much too weak because of the very low perturbation amplitude in this region.
Open with DEXTER

	Figure 9: Amplitude of the Fourier components of the density perturbation. The odd modes are numerically zero. Due to small nonlinearities, the even modes are apparent but with weak amplitude. The components $m\ge 9$ are set to zero because of the desaliasing process.
Open with DEXTER

The nonlinearities are weak throughout the simulation. Indeed, looking at the Fourier spectrum of the density in Fig. 9, the amplitude of each component, plotted vs. the mode m on a logarithmic scale, is small except for that corresponding to the excitation mode. The odd modes are not present. However, the weak nonlinearities create a cascade to high even modes starting with m=2. The largest asymmetric expansion coefficient C_m is m=2, while the next even coefficients roughly follow a geometric series with a factor q=10^-3, so we can write for all m even, $C_m \approx q^{m/2-1} ~ C_2$, until they reach values less than 10^-20, which can be interpreted as zero from a numerical point of view. As deviation from the stationary state is weak, the amplitudes of these even modes decay compared to the previous one, the highest being of course m=2. Thus, even in the full non-linear simulation, the regime remains quasi-linear. Consequently, the parametric resonance phenomenon discussed in the previous section is irrelevant at this stage of our work. The effects of strong nonlinearity will be studied in a forthcoming paper. Note that due to the desaliasing process, the modes $m\ge 9$ are all set to zero. Note also that the free wave solutions leave the computational domain and do not return. Only the non-wavelike disturbance produces significant changes in the density profile.

	Figure 10: Same as in Fig. 7 but the Newtonian potential is replaced by a pseudo-Schwarzschild one.
Open with DEXTER

   
5.4 Pseudo-Schwarzschild potential

In order to take into account the modification of the radial epicyclic frequency due to the curved space-time around a Schwarzschild black hole, we replaced the Newtonian potential by the logarithmically modified potential proposed by Mukhopadhyay (2003). This potential is well-suited to approximating the angular and epicyclic frequencies in accretion disks around a rotating black hole. The radial gravitational field derived from this potential is given by

 

$$g_r = - \frac{G~M_*}{r^2} \bigg[ 1 + R_{ms} \bigg( \frac{9}{20} ~... ...ac{R_{ms}-1}{r} - \frac{3}{2r} ~ \ln \frac{r}{(3r-R_{ms})^{2/9}} \bigg) \bigg],$$ (97)

where R_ms is the last stable circular orbit:

$$3+Z_2 \pm \sqrt{(3-Z_1)(3+Z_1+2Z_2)}$$                                  R_ms = (98)

Z₁ = 1 + (1-a²)^1/3   [ (1+a)^1/3 + (1-a)^1/3 ] (99)

$$\sqrt{3~a^2 + Z_1^2},$$ Z₂ = (100)

where a is related to the angular momentum of the star.

The important feature of this potential is that the radial epicyclic frequency for a single particle vanishes at the innermost stable circular orbit (ISCO), which implies that the parametric resonance condition can be fulfilled in the very inner part of the accretion disk for different values of the integer n but still with a very small radial extension. However, due to the small perturbation we remain in the linear regime, and the parametric resonance is only of second order. The disk boundaries, the star rotation rate, and the normalization are the same as in the previous case.

The final snapshot of the density perturbation in the disk is shown in Fig. 10. The corotation resonance is now located at r=43.7, which differs slightly from the previous simulation because the orbital velocity is no longer Keplerian due to the pseudo-Newtonian force Eq. (97). The non-wavelike disturbance exists between the inner edge of the disk and the inner Lindblad resonance. By viewing a cross section, the Lindblad radii appear clearly as seen in Fig. 11. Inspection of the Fourier coefficients of the density Fig. 12 proves that the linear regime is preserved. The level of excitation of the even modes $m\ne2$ is very weak.

	Figure 11: Same as in Fig. 8 but for the pseudo-Schwarzschild potential. The Lindblad resonances are located at  $r_{\rm L}^{\rm in/out}=21.6/45.5$.
Open with DEXTER

	Figure 12: Same as Fig. 9 but for the pseudo-Schwarzschild potential. The Fourier coefficients form a decaying geometric series.
Open with DEXTER

	Figure 13: Same as in Fig. 10 but with an azimuthal mode m=1.
Open with DEXTER

	Figure 14: Same as in Fig. 11 but with an azimuthal mode m=1. The Lindblad resonances are located at  $r_{\rm L}^{\rm in/out}=12.9/56.0$.
Open with DEXTER

	Figure 15: Same as Fig. 12 but with an azimuthal mode m=1. The Fourier coefficients form a decaying geometric series.
Open with DEXTER

The results of replacing the m=2 mode by a more realistic m=1mode, in accordance with the dipolar structure of the magnetic field, is shown in Fig. 13. The inner Lindblad resonance moves closer to the inner edge and is located at  $R_{\rm L}^{\rm in}=12.9$. There is less space left for the non-wavelike disturbance between the inner edge and the inner Lindblad resonance as can be seen in Fig. 14. Also, all the azimuthal modes are now excited starting with m=1 and with decreasing magnitude, Fig. 15. From a qualitative point of view, we can draw the same conclusions as in the case of the Newtonian disk. The absence or presence of the ISCO does not affect the physics of the resonances, but only has an effect on their amplitude and their location.

   
5.5 Pseudo-Kerr potential

We also ran a simulation in which the rotation of the star is taken into account. This shifts the location of the ISCO closer to the surface of the neutron star.

We chose a star with an angular momentum of a_*=0.5. Therefore, the disk inner boundary corresponding to the marginally stable circular orbit will be located at R₁=4.24, while the outer boundary is chosen to be at R₂=42.4. The star magnetic perturbation in this run is ten times larger than in the previous case, $\delta b = 0.01$. Therefore the computation starts to leave the linear regime and reaches a weak non-linear state.

The corotation resonance at r=23.8 is not reached, as can be seen from Fig. 16. However, the wave and non-wavelike disturbances are clearly identified between the inner edge and the inner Lindblad resonance. They persist during the whole simulation at the low level predicted by the linear analysis. The precise location of the Lindblad resonance is determined by the cross section view shown in Fig. 17. Looking at the Fourier components in Fig. 18, all the even modes are excited at a substantial level. The stronger stellar perturbation starts to excite modes which are multiples of the injection scale given by m=2.

	Figure 16: Final snapshot of the density perturbation in the accretion disk evolving in a perturbed pseudo-Kerr potential with a=0.5. The disk extends from R1=4.24 to R2=42.4. The outer Lindblad resonance is not on the grid.
Open with DEXTER

	Figure 17: Same as in Fig. 8 but for the pseudo-Kerr potential. The inner Lindblad resonance appears clearly at  $r_{\rm L}^{\rm in}=19.2$ while the outer one at  $r_{\rm L}^{\rm out}=45.8$ lies outside the computational grid.
Open with DEXTER

	Figure 18: Same as Fig. 9 but for the pseudo-Kerr potential.
Open with DEXTER

5.6 Retrograde disk

Finally we checked that the Lindblad resonances disappear for a retrograde Newtonian disk. The numerical values are the same as in Sect. 5.3, except that the sign of the stellar spin is taken to be  $\Omega_*=-0.0043311$ and the magnetic perturbation  $\delta b = 0.01$.

Looking at Fig. 19, the evolution of the disk does not show any Lindblad resonance as expected. The mode m=2 with two spiral arms propagates in the whole computational domain. The disk reaches a quasi-stationary state in which the density fluctuations are well established, Fig. 20. They fit the shape found in the linear analysis, Fig. 4. Again, the stronger stellar magnetic perturbation affects the high order azimuthal modes, which are excited at the same amplitude as the injection mode m=2, Fig. 21.

	Figure 19: Final snapshot of the density perturbation in the counterrotating accretion disk evolving in a perturbed Newtonian potential. Same values as in Fig. 7 apply except for the sign of $\Omega _*$. A trailing spiral density wave of m=2 propagates in the entire disk.
Open with DEXTER

	Figure 20: Same as Fig. 8 but for the retrograde disk. The Lindblad resonances are no longer present.
Open with DEXTER

	Figure 21: Same as Fig. 9 but with a retrograde disk.
Open with DEXTER

Here also we can replace the m=2 mode by a more realistic m=1 one, in accordance with the dipolar structure of the magnetic field and with, for example, a smaller magnetic perturbation taken again to be  $\delta b = 0.001$. The results for the density perturbation are shown in Fig. 22. A spiral wave is formed and propagates in the whole disk with a rotation rate equal to the neutron star spin. Now all the Fourier components are excited, the strongest one corresponds to m=1, and the lighter order ones have decreasing amplitudes, as shown in Fig. 24.

	Figure 22: Final snapshot of the density perturbation in the counterrotating accretion disk evolving in a perturbed Newtonian potential. Same values as in the Fig. 7 caption excepted for the sign of $\Omega _*$. A trailing spiral density wave of m=1 is propagating in the whole disk.
Open with DEXTER

	Figure 23: Same as Fig. 8 but with a retrograde disk. The Lindblad resonances have disappeared.
Open with DEXTER

	Figure 24: Same as Fig. 9 but with a retrograde disk.
Open with DEXTER

To summarize, in all of the results obtained by these runs, the exact structure of the gravitational potential, be it Newtonian or pseudo-Newtonian, does not affect the physics of the resonances. The main features of these resonances are their location very close to the inner boundary of the accretion disk. When the applied magnetic asymmetry is small enough, the perturbations remain in the linear regime.

   
6 Conclusion

In this paper, we studied the effects of an asymmetric rotating dipolar or multipolar magnetic field on a magnetized accretion disk. This situation is well-suited to LMXBs containing a neutron star, since they are believed to possess a magnetic field that more or less channels the accreting matter to the polar caps of the neutron star, depending on its strength. The same conclusions can be drawn in the case of a hydrodynamical disk evolving in a quadrupolar gravitational potential perturbation. Although the flows differ, the consequences of a rotating asymmetric field are qualitatively the same. The disk is subject to the corotation resonance, the Lindblad driven resonances, and if the perturbation is strong enough, the parametric instability. We derived an explicit formulation for the regions where the blobs should form, given by the resonance conditions Eqs. (67)-(69). On the one hand, these blobs can persist during thousands of orbital revolutions and can account for the high quality factor Q seen in some high frequency-QPOs. On the other, the pattern speed of the free wave solutions is closely related to the location of the inner boundary of the accretion disk. The observed shift in the HFQPOs in response to the accretion rate can be explained by these solutions.

This study should be extended to include viscosity of the gas flowing in the disk, which should set a lower limit to the magnetic field perturbation amplitude in order to have parametric resonance. In addition, the loss of energy by radiation should also damp these oscillations and select only a few of them, which will also permit estimation of the light curves of such systems. The relation between the simulations presented in this paper and observations will be investigated in a forthcoming paper, in which we show how to calculate the power spectrum density (PSD) of the accretion disk. The prediction of the PSD will be crucial for validation of the model.

We would like to emphasize that this work is only a preliminary study undertaken to prove that a rotating asymmetric magnetic field leads to some interesting resonance mechanisms. The toy model used here cannot be applied directly to observations. Indeed, many aspects, such as the warping/precession of the disk expected in 3D simulations, as well as viscosity and radiation, should be taken into account in a more realistic model.

Nevertheless, we believe that this idea can also be applied to the LMXBs containing a black hole. Indeed, it seems that the existence of QPOs in black hole candidates is correlated with the detection of a jet emanating from these systems, (McClintock & Remillard 2003). The Blandford-Znajek process (Blandford & Znajek 1977) converting the rotational energy of the black hole into the launching of the jet should also be responsible for a connection between the accretion disk and the black hole via the magnetic field. Therefore, the model presented in this paper has the potential to explain the QPOs seen in the BHCs (Wang et al. 2003).

Acknowledgements

This research was carried out in an FOM projectruimte on "Magnetoseismology of accretion disks'', a collaborative project between R. Keppens (FOM Institute Rijnhuizen, Nieuwegein) and N. Langer (Astronomical Institute Utrecht). This work is part of the research programme of the "Stichting voor Fundamenteel Onderzoek der Materie (FOM)'', which is financially supported by the "Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)''.
I am grateful to the referee Lev Titarchuk for his valuable comments and remarks as well as to Jean Heyvaerts and John Kirk for carefully reading the manuscript.

References

• Alpar, M. A., & Shaham, J. 1985, Nature, 316, 239 [NASA ADS] [CrossRef]
• Anzer, U., & Börner, G. 1980, A&A, 83, 133 [NASA ADS] (In the text)
• Anzer, U., & Börner, G. 1983, A&A, 122, 73 [NASA ADS] (In the text)
• Blandford, R. D., & Znajek, R. L., 1977, MNRAS, 179, 433 [NASA ADS] (In the text)
• Belloni, T., Psaltis, D., & van der Klis, M. 2002, ApJ, 572, 392 [NASA ADS] [CrossRef] (In the text)
• Canuto, C., Hussaini, M., Quarteroni, A., & Zang, T. 1990, Spectral Methods in Fluid Dynamics (Berlin and New York: Springer Verlag) (In the text)
• Lee, W. H., Abramowicz, M. A., & Kluzniak, W. 2004, ApJ, 603, L93 [NASA ADS] [CrossRef]
• Frieman, E., & Rotenberg, M. 1960, Rev. Mod. Phys., 32, 898 [NASA ADS] [CrossRef] (In the text)
• Ghosh, P., Pethick, C. J., & Lamb, F. K. 1977, ApJ, 217, 578 [NASA ADS] [CrossRef]
• Ghosh, P., & Lamb, F. K. 1978, ApJ, 223, L83 [NASA ADS] [CrossRef] (In the text)
• Ghosh, P., & Lamb, F. K. 1979a, ApJ, 232, 259 [NASA ADS] [CrossRef] (In the text)
• Ghosh, P., & Lamb, F. K. 1979b, ApJ, 234, 296 [NASA ADS] [CrossRef]
• Goldreich, P., Tremaine, S. 1979, ApJ, 233, 857 [NASA ADS] [CrossRef] (In the text)
• Gottlieb, D., & Orszag, S. 1977, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM (In the text)
• Kato, S. 2001a, PASJ, 53, 1 [NASA ADS] (In the text)
• Kato, S. 2001b, PASJ, 53, L37 [NASA ADS]
• Kato, S. 2004, PASJ, 56, 905 [NASA ADS] (In the text)
• Kluzniak, W., Abramowicz, M. A., Kato, S., et al. 2004, ApJ, 603, L89 [NASA ADS] [CrossRef] (In the text)
• Lai, D. 1999, ApJ, 524, 1030 [NASA ADS] [CrossRef] (In the text)
• Landau, L., & Lifshitz, E. 1982, Cours de physique théorique, Tome I: mécanique (Moscou: Éditions Mir) (In the text)
• Lee, W. H., Abramowicz, M. A., & Kluzniak, W. 2004, ApJ, 603, L93 [NASA ADS] [CrossRef] (In the text)
• Lindblad, P. O. 1974, IAU Symp. 58, 399 (In the text)
• Mauche, C. W. 2002, ApJ, 580, 423 [NASA ADS] [CrossRef] (In the text)
• McClintock, J. E., & Remillard, R. A. 2003,[arXiv:astroph/0306213] (In the text)
• Miller, M. C., & Lamb, F. K., & Psaltis, D. 1998, ApJ, 508, 791 [NASA ADS] [CrossRef]
• Morse, & Feshbach, 1953, Methods of Theoretical Physics (New-York: McGraw Hill) (In the text)
• Mukhopadhyay, B., & Misra, R. 2003, ApJ, 582, 347 [NASA ADS] [CrossRef] (In the text)
• Psaltis, D., Belloni, T., & van der Klis, M. 1999, ApJ, 520, 262 [NASA ADS] [CrossRef] (In the text)
• Schnittman, J. D., & Bertschinger, E. 2004, ApJ, 606, 1098 [NASA ADS] [CrossRef] (In the text)
• Shirakawa, A., & Lai, D. 2002a, ApJ, 564, 361 [NASA ADS] [CrossRef] (In the text)
• Shirakawa, A., & Lai, D. 2002b, ApJ, 565, 1134 [NASA ADS] [CrossRef] (In the text)
• Smirnov, V. 1989, Cours de mathématiques supérieures, Tome 3, Vol. 2 (Moscou: Éditions Mir) (In the text)
• Stella, L., & Vietri, M. 1998, ApJ, 492, L59 [NASA ADS] [CrossRef] (In the text)
• Stella, L., & Vietri, M. 1999, Phys. Rev. Lett., 82, 17 [NASA ADS] [CrossRef]
• Stella, L., Vietri, M., & Morsink, S. 1999, ApJ, 524, L63 [NASA ADS] [CrossRef] (In the text)
• Titarchuk, L. 2003, ApJ, 591, 354 [NASA ADS] [CrossRef] (In the text)
• Titarchuk, L., & Fiorito, R. 2004, ApJ, 612, 988 [NASA ADS] [CrossRef] (In the text)
• Titarchuk, L., & Wood, K. 2002, ApJ, 577, L23 [NASA ADS] [CrossRef] (In the text)
• Titarchuk, L., Lapidus, I., & Muslimov, A. 1998, ApJ, 499, 315 [NASA ADS] [CrossRef]
• van der Klis, M. 2000, ARA&A, 38, 717 [NASA ADS] (In the text)
• Vietri, M., & Stella, L. 1998, ApJ, 503, 350 [NASA ADS] [CrossRef] (In the text)
• Wang, D., Ma, R., Lei, W. & Yao, G. 2003, MNRAS, 344, 473 [NASA ADS] [CrossRef] (In the text)
• Wijnands, R. 2001, Adv. Space Res., 28, 469 [NASA ADS] (In the text)
• Zhang, W., Smale, A. P., Strohmayer, T. E., & Swank, J. H. 1998, ApJ, 500, L171 [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

   
Appendix A: Eigenvalue problem for the Lagrangian displacement

A.1 The Lagrangian displacement $\vec{\xi}$

The linearization of the MHD equations is done by introducing the Lagrangian displacement vector $\vec{\xi}$, which is related to the Eulerian and Lagrangian perturbation of the fluid physical quantities, such as the velocity, the magnetic field, and the density. These perturbations are denoted respectively by $\delta$ and $\Delta$. By definition, the Lagrangian perturbation of the velocity is

$$\Delta\vec{v} = \vec{v} - \vec{v}_0 = \frac{{\rm d}\vec{\xi}}... ...al\vec{\xi}}{\partial t} + \vec{v}\cdot\vec{\nabla}~\vec{\xi}.$$ (A.1)

Moreover, the Eulerian perturbation is related to the Lagrangian one by

$$\delta = \Delta - \vec{\xi}\cdot\nabla.$$ (A.2)

Therefore, the Eulerian perturbation of the velocity reads

$$\delta\vec{v} = \vec{v}_1 = \frac{\partial\vec{\xi}}{\parti... ...vec{\nabla}~\vec{\xi} - \vec{\xi} \cdot \vec{\nabla}~\vec{v}.$$ (A.3)

The boundary conditions at the disk-vacuum interface are

 

$$\left\langle p + \frac{\vec{B}^2}{2~\mu_0} \right\rangle$$ (A.4)

$$\vec{n} \cdot \vec{B}$$ (A.5)

$$\vec{n} \wedge \left\langle\vec{B} \right\rangle$$ (A.6)

where  $\left\langle X\right\rangle=X_{\rm vac}-X_{\rm disk}$ denotes the jump at the interface and $\vec{n}$ is a unit vector normal to the plasma interface.

Let's start by perturbing the mass conservation equation by using the equilibrium condition  ${\rm div}~(\rho_0\vec{v}_0)=0$

$$\frac{\partial\rho_1}{\partial t} + \vec{v}_0 \cdot \nabla\rho_1 ... ...{v}\cdot\vec{\nabla}~\vec{\xi} - \vec{\xi} \cdot \vec{\nabla}~\vec{v}) \right).$$ (A.7)

From the tensorial calculus formalism, the right hand side can be transformed after tedious algebra into

$$\frac{\partial\rho_1}{\partial t} + \vec{v}_0 \cdot \nabla\rho_1 ... ... div}~(\rho_0~\vec{\xi})) - {\rm div}~(\rho_0~\vec{\xi}) ~ {\rm div}~\vec{v}_0.$$ (A.8)

Because the Lagrangian displacement is arbitrary, we make the following identification

$$\rho_1 = -{\rm div}~(\rho_0~\vec{\xi}).$$ (A.9)

To continue by perturbing the induction equation, we get

$$\frac{\partial\vec{B}_1}{\partial t} + \vec{v}_0 \cdot \nabla \ve... ...\vec{v}_1 - \vec{B}_0 ~ {\rm div}~\vec{v}_1 - \vec{v}_1 \cdot \nabla \vec{B}_0.$$ (A.10)

We transform each term on the right hand side from the equilibrium state $\bf {rot}~(\vec{v}_0\wedge\vec{B}_0)=0$. The unperturbed magnetic field is divergenceless neither is the perturbed one, ${\rm div}~\vec{B}_0={\rm div}~\vec{B}_1=0$. Thus, we obtain

$$\vec{B}_0 \cdot \nabla \vec{v}_1 \frac{\partial}{\partial t} (\vec{B}_0\cdot\nabla\vec{\xi}) + \vec{v}_0\cdot\nabla(\vec{B}_0\cdot\nabla\vec{\xi})$$ =

$$+ {\rm div}~\vec{v}_0 ~ \vec{B}_0\cdot\nabla\vec{\xi} - {\rm div}~(\vec{B}_0\vec{\xi}\cdot\nabla\vec{v}_0 )$$ (A.11)

$$- \vec{B}_0 ~ {\rm div}~\vec{v}_1 - \frac{\partial}{\partial t}(\vec{B}_0{\rm div}~\vec{\xi}) - \vec{v}_0\cdot\nabla(\vec{B}_0{\rm div}~\vec{\xi})$$ =

$$+ (\vec{v}_0\cdot\nabla\vec{B}_0){\rm div}~\vec{\xi} + \vec{B}_0 ~\vec{\xi}\nabla {\rm div}~\vec{v}_0$$ (A.12)

$$- \vec{v}_1 \cdot \nabla \vec{B}_0 - \frac{\partial}{\partial t} (\vec{\xi}\cdot\nabla\vec{B}_0) - \vec{v}_0\cdot\nabla\vec{\xi}\cdot\nabla\vec{B}_0$$ =

$$+ \vec{\xi}\cdot\nabla\left[ \vec{v}_0\cdot\nabla\vec{B}_0\right].$$ (A.13)

By introducing the vector

$$\vec{Q}=\bf {rot}~(\vec{\xi}\wedge\vec{B}_0) ,$$ (A.14)

the sum of the three terms can be cast after tedious tensorial calculus into the form

$$\frac{\partial\vec{Q}}{\partial t} + \vec{v}_0\cdot\nabla\vec... ...vec{Q} ~ {\rm div}~\vec{v}_0 - \vec{Q} \cdot \nabla \vec{v}_0.$$ (A.15)

So we make the identification  $\vec{B}_1=\vec{Q}$. This is the perturbation of the magnetic field induced by motion in the disk and does not include the magnetic field perturbation emanating from the star.

Finally for the pressure perturbation, we develop the adiabatic condition to first order. We have

$$\frac{\partial p_1}{\partial t} + \vec{v}_0 \cdot \nabla p_1 ... ...ma ~ p_0 ~ {\rm div}~\vec{v}_1 - \vec{v}_1 \cdot \nabla p_0.$$ (A.16)

Introducing $\pi=-\gamma~p_0~{\rm div}~\vec{\xi}-\vec{\xi}\cdot\vec{\nabla}~p_0$, the right hand side is

$$\frac{\partial \pi}{\partial t} + \vec{v}_0 \cdot \nabla \pi + \gamma ~ \pi ~ {\rm div}~\vec{v}_0$$ (A.17)

and thus the pressure perturbation of the gas is

$$p_1=\delta p = -\gamma~p_0~{\rm div}~\vec{\xi}-\vec{\xi}\cdot\vec{\nabla}~p_0.$$ (A.18)

Equations (A.3, A.14, A.18) are the perturbed quantities induced by the motion in the fluid. However, the magnetic field perturbation possesses another variation that is not imposed by the fluid but is introduced independently by the neutron star. This known perturbation  $\delta \vec{B}_*$ has to be added to the momentum equation. In the realistic case of a dipolar magnetic field, we should expect  $\bf {rot}~\vec{B}_*=\vec{0}$, as well as  $\bf {rot}~\delta\vec{B}_*=\vec{0}$. Nevertheless, in our 2D numerical simulations, we do not keep this property for the stellar magnetic field because in this particular geometry, we impose straight field lines aligned along the rotational z-axis. Therefore, in order to treat the problem in a self-consistent manner we keep this non-vanishing curl term in the linear analysis already in order to facilitate the connection linear analysis-numerical simulations. To get the eigenvalues system satisfied by $\vec{\xi}$ we put expressions (A.14, A.18) into the momentum Eq. (16).

A.2 PDE for $\vec{\xi}$

Replacing all the perturbed quantities by their expression involving the Lagrangian displacement and keeping only the first order term, $\vec{\xi}$, satisfies a second order linear PDE given by

$$\rho_0 ~ \frac{\partial\vec{v}_1}{\partial t} + \rho_0 ~ \vec{v}_... ... \vec{v}_1 \cdot \nabla \vec{v}_0 + \rho_1 ~ \vec{v}_0 \cdot \nabla \vec{v}_0 =$$

$$\rho_0 \delta\vec{g} + \rho_1 ~ ( \vec{g}_0 + \delta\vec{g} ) - \nabla p_1$$

$$+ \frac{1}{\mu_0} ~ \left( \bf {rot}~\vec{B}_0 \wedge (\vec{B}_1 + \delta\vec{B}_*) + \bf {rot}~\vec{B}_1 \wedge (\vec{B}_0 \right.$$

$$+ \left. \delta\vec{B}_*) + \bf {rot}~\delta\vec{B}_* \wedge (\vec{B}_0 + \vec{B}_1 + \delta\vec{B}_*) \right).$$ (A.19)

On the left hand side, we have for each term

$$\rho_0 ~ \frac{\partial\vec{v}_1}{\partial t} \rho_0 ~ \frac{\partial^2\vec{\xi}}{\partial t^2} + \rho_0 ~ \vec... ...ial t} - \rho_0 ~ \frac{\partial\vec{\xi}}{\partial t} ~ \cdot \nabla \vec{v}_0$$  =   (A.20)

$$\rho_0 ~ \vec{v}_0 \cdot \nabla \vec{v}_1 \rho_0 ~ \vec{v}_0 \cdot \nabla \frac{\partial\vec{\xi}}{\partial... ... \nabla ( \vec{v}_0 \cdot \nabla \vec{\xi} - \vec{\xi} \cdot \nabla \vec{v}_0 )$$  =   (A.21)

$$\rho_0 ~ \vec{v}_1 \cdot \nabla \vec{v}_0 \rho_0 ~ \frac{\partial\vec{\xi}}{\partial t} ~ \cdot \nabla \vec... ...\rho_0 ~ \left( \vec{v}_0 \cdot \nabla \vec{\xi} \cdot \nabla \right) \vec{v}_0$$  =  

$$- \rho_0 ~ \left( \vec{\xi} \cdot \nabla \vec{v}_0 \cdot \nabla \right) \vec{v}_0$$ (A.22)

$$\rho_1 ~ \vec{v}_0 \cdot \nabla \vec{v}_0 - {\rm div}~(\rho_0~\vec{\xi}) ~ \vec{v}_0 \cdot \nabla \vec{v}_0.$$ = (A.23)

Finally, summing these four terms gives

$$\rho_0 ~ \frac{\partial\vec{v}_1}{\partial t} + \rho_0 ~ \vec{v}_... ... \vec{v}_1 \cdot \nabla \vec{v}_0 + \rho_1 ~ \vec{v}_0 \cdot \nabla \vec{v}_0 =$$

$$\rho_0 ~ \frac{\partial^2\vec{\xi}}{\partial t^2} + 2 ~ \rho_0 ~ ... ... ) \vec{\xi} - \rho_0 ~ \vec{\xi} ( \vec{v}_0 \cdot \nabla ) \vec{v}_0 \right].$$ (A.24)

Recalling that  $\vec{Q}=\vec{B}_1$ and $p_1=\pi$, we rewrite

$$\bf {rot}~\vec{B}_0 \wedge \vec{Q} + \bf {rot}~\vec{Q} \wedge \ve... ...0\cdot\vec{Q}) + \vec{B}_0\cdot \nabla \vec{Q} + \vec{Q} \cdot \nabla \vec{B}_0$$ (A.25)

Introducing the opposite of the total pressure variation, gaseous and magnetic, by  $\Pi = \gamma ~ p_0 ~ {\rm div}~\vec{\xi} + \vec{\xi} \cdot \vec{\nabla}~p_0 - \frac{1}{\mu_0} ~ \vec{B}_0 \cdot \vec{Q}$, the Lagrangian displacement for the fluid motion in the accretion disk satisfies the following second order linear partial differential equation

 

$$\rho_0 ~ \frac{\partial^2\vec{\xi}}{\partial t^2} + 2 ~ \rho_0 ~ ... ...a ) \vec{\xi} - \rho_0 ~ \vec{\xi} ( \vec{v}_0 \cdot \nabla ) \vec{v}_0 \right]$$

$$- \nabla\Pi - \frac{1}{\mu_0} ~ \vec{B}_0 \cdot \nabla \vec{Q} - \frac{1}{\mu_0} ~ \vec{Q} \cdot \nabla \vec{B}_0$$

$$+ {\rm div}~( \rho_0 ~ \vec{\xi} ) ( \vec{g}_0 + \delta\vec{g} ) - \rho_0 ~ \delta\vec{g} =$$

$$\frac{1}{\mu_0} ~ \left[ \bf {rot}~( \vec{B}_0 + \vec{Q} + \delta... ... \vec{B}_* + \bf {rot}~\delta \vec{B}_* \wedge ( \vec{B}_0 + \vec{Q} ) \right].$$ (A.26)

In the expression above, we allow for a perturbation in the gravitational field induced by the star and denoted by  $\delta\vec{g}$. However, in the following analysis, we remove this freedom and only keep a perturbation in the magnetic field.

A.3 Projection onto the cylindrical coordinate system

In the stationary state, the disk orbits around the compact object without mass in- or outflow from the edges and without precession. From now on, we drop the subscript 0 for the stationary quantities, which should not lead to any confusion. Thus the velocity in the axisymmetric stationary state is simply

$$\vec{v} = v_\varphi ~ \vec{e}_\varphi = r ~ \Omega ~ \vec{e}_\varphi.$$ (A.27)

We find the following expressions

$$( \vec{v} \cdot \nabla ) \vec{v} - r ~ \Omega^2 ~ \vec{e}_r$$  =   (A.28)

$$( \vec{v} \cdot \nabla ) \vec{\xi} \Omega ~ \left[ \left( \frac{\partial\xi_r}{\partial\varphi} - \x... ...) ~ \vec{e}_\varphi + \frac{\partial\xi_z}{\partial\varphi} ~ \vec{e}_z \right]$$  =   (A.29)

$$( \vec{\xi} \cdot \nabla ) \vec{v} - \Omega ~ \xi_\varphi ~ \vec{e}_r + \vec{\xi} \cdot \vec{\nabla}~(r~\Omega) ~ \vec{e}_\varphi.$$   =  (A.30)

We also give the explicit formulas for each term in Eq. (A.26).

$$\rho ~ \frac{\partial^2\vec{\xi}}{\partial t^2} \rho ~ \frac{\partial^2\xi_r}{\partial t^2} ~ \vec{e}_r + \rho ~ ... ...^2} ~ \vec{e}_\varphi + \rho ~ \frac{\partial^2\xi_z}{\partial t^2} ~ \vec{e}_z$$  =  (A.31)

$$2 ~ \rho ~ \vec{v} \cdot \nabla \frac{\partial\vec{\xi}}{\partial... ...\frac{\partial\xi_r}{\partial\varphi} - \xi_\varphi \right) ~ \vec{e}_r \right.$$

$$+ \left. \left( \frac{\partial\xi_\varphi}{\partial\varphi} + \xi... ...) ~ \vec{e}_\varphi + \frac{\partial\xi_z}{\partial\varphi} ~ \vec{e}_z \right]$$ (A.32)

$$\vec{\nabla}~\Pi \frac{\partial\Pi}{\partial r} ~ \vec{e}_r + \frac{1}{r} ~ \frac{... ...vec{e}_z - {\rm div}~( \rho ~ \vec{\xi} ~ (\vec{v} \cdot \vec{\nabla}) \vec{v})$$ = (A.33)

$$\left[ r ~ \Omega^2 ~ {\rm div}~(\rho~\vec{\xi}) + \rho ~ \vec{\x... ...Omega^2 ) \right] ~ \vec{e}_r + \rho ~ \xi_\varphi ~ \Omega^2 ~ \vec{e}_\varphi$$  =  (A.34)

$${\rm div}~(\rho ~ \vec{v} ~ ( \vec{v} \cdot \vec{\nabla} ) \vec{\xi}) \rho ~ \Omega^2 ~ \frac{\partial^2 \vec{\xi}}{\partial\varphi^2}$$  =  (A.35)

$$\rho ~ \Omega^2 ~ \left[ \left( \frac{\partial^2 \xi_r}{\partial\... ...r - 2 ~ \frac{\partial\xi_\varphi}{\partial\varphi} \right) ~ \vec{e}_r \right.$$  = 

$$+ \left. \left( \frac{\partial^2 \xi_\varphi}{\partial\varphi^2} ... ...vec{e}_\varphi + \frac{\partial^2 \xi_z}{\partial\varphi^2} ~ \vec{e}_z \right]$$

$$\vec{B} \cdot \vec{\nabla}~\vec{Q} \left[ B_r ~ \frac{\partial Q_r}{\partial r} + \frac{B_\varphi}{r... ...arphi~Q_\varphi}{r} + B_z ~ \frac{\partial Q_r}{\partial z} \right] ~ \vec{e}_r$$  = 

$$\quad + \left[ B_r ~ \frac{\partial Q_\varphi}{\partial r} + \fra... ...Q_r}{r} + B_z ~ \frac{\partial Q_\varphi}{\partial z} \right] ~ \vec{e}_\varphi$$

$$\quad+ \left[ B_r ~ \frac{\partial Q_z}{\partial r} + \frac{B_\va... ...}{\partial \varphi} + B_z ~ \frac{\partial Q_z}{\partial z} \right] ~ \vec{e}_z$$ (A.36)

$$\vec{Q} \cdot \vec{\nabla}~\vec{B} \left[ Q_r ~ \frac{\partial B_r}{\partial r} + \frac{Q_\varphi}{r... ...arphi~B_\varphi}{r} + Q_z ~ \frac{\partial B_r}{\partial z} \right] ~ \vec{e}_r$$  = 

$$\quad+ \left[ Q_r ~ \frac{\partial B_\varphi}{\partial r} + \frac... ...B_r}{r} + Q_z ~ \frac{\partial B_\varphi}{\partial z} \right] ~ \vec{e}_\varphi$$

$$\quad+ \left[ Q_r ~ \frac{\partial B_z}{\partial r} + \frac{Q_\va... ...}{\partial \varphi} + Q_z ~ \frac{\partial B_z}{\partial z} \right] ~ \vec{e}_z$$ (A.37)

$$\rho ~ \delta\vec{g} ~=~ \rho ~ ( \delta g_r ~ \vec{e}_r + \delta g_\varphi ~ \vec{e}_\varphi + \delta g_z ~ \vec{e}_z )$$ (A.38)

$$- {\rm div}~(\rho~\vec{\xi} ) ~ ( \vec{g} + \delta\vec{g} ) =$$

$$- {\rm div}~(\rho~\vec{\xi} ) ~ [ ( g_r + \delta g_r ) ~ \vec{e}_r + \delta g_\varphi ~ \vec{e}_\varphi + ( g_z + \delta g_z ) ~ \vec{e}_z ]$$ (A.39)

$$\bf {rot}~( \vec{B} + \vec{Q} + \delta \vec{B}_* ) \wedge \delta ... ...c{\partial}{\partial z} \left( B_r + Q_r + \delta B_*^r \right) \right. \right.$$

$$- \left.\left. \frac{\partial}{\partial r} \left( B_z + Q_z + \delta B_*^z \right) \right\} \right.$$

$$- \left. \frac{\delta B_*^\varphi}{r} ~ \left\{ \frac{\partial}{\... ...eft( r ~ ( B_\varphi + Q_\varphi + \delta B_*^\varphi ) \right) \right. \right.$$

$$- \left.\left. \frac{\partial}{\partial\varphi} ( B_r + Q_r + \delta B_*^r ) \right\} \right] ~ \vec{e}_r$$

$$\left[ \frac{\delta B_*^r}{r} ~ \left\{ \frac{\partial}{\partial ... ...eft( r ~ ( B_\varphi + Q_\varphi + \delta B_*^\varphi ) \right) \right. \right.$$

$$- \left.\left. \frac{\partial}{\partial\varphi} ( B_r + Q_r + \delta B_*^r ) \right\} \right.$$

$$- \left. \delta B_*^z \left\{ \frac{1}{r} \frac{\partial}{\partial\varphi} ( B_z + Q_z + \delta B_*^z ) \right. \right.$$

$$- \left.\left. \frac{\partial}{\partial z} ( B_\varphi + Q_\varphi + \delta B_*^\varphi )\right\} \right] ~ \vec{e}_\varphi$$

$$\left[ \delta B_*^\varphi \left\{ \frac{1}{r} \frac{\partial}{\partial\varphi} ( B_z + Q_z + \delta B_*^z ) \right. \right.$$

$$- \left. \left. \frac{\partial}{\partial z} ( B_\varphi + Q_\varphi + \delta B_*^\varphi )\right\} \right.$$

$$- \left. \delta B_*^r ~ \left\{ \frac{\partial}{\partial z} \left( B_r + Q_r + \delta B_*^r \right) \right. \right.$$

$$- \left.\left. \frac{\partial}{\partial r} \left( B_z + Q_z + \delta B_*^z \right) \right\} \right] ~ \vec{e}_z$$ (A.40)

$$\bf {rot}~\delta\vec{B}_* \wedge ( \vec{B} + \vec{Q} ) = \left[ (... ...ta B_*^r}{\partial z} - \frac{\partial\delta B_*^z}{\partial r} \right) \right.$$

$$- \left. \frac{B_\varphi + Q_\varphi}{r} ~ \left( \frac{\partial}... ...phi) - \frac{\partial\delta B_*^r}{\partial\varphi} \right) \right] ~ \vec{e}_r$$

$$\left[ \frac{B_r + Q_r}{r} ~ \left( \frac{\partial}{\partial r}(r... ...lta B_*^\varphi) - \frac{\partial\delta B_*^r}{\partial\varphi} \right) \right.$$

$$- \left. (B_z + Q_z) ~ \left( \frac{1}{r} ~ \frac{\partial\delta ... ...\frac{\partial\delta B_*^\varphi}{\partial z} \right) \right] ~ \vec{e}_\varphi$$

$$\left[ ( B_\varphi + Q_\varphi ) \left( \frac{1}{r} ~ \frac{\part... ...partial\varphi} - \frac{\partial\delta B_*^\varphi}{\partial z} \right) \right.$$

$$- \left. (B_r + Q_r) \left( \frac{\partial\delta B_*^r}{\partial z} - \frac{\partial\delta B_*^z}{\partial r} \right) \right] ~ \vec{e}_z.$$ (A.41)

The perturbations in the disk magnetic field induced by the fluid are in cylindrical coordinates

$$\frac{1}{r} ~ \frac{\partial}{\partial\varphi} \left( \xi_r ~ B_\... ... \right) - \frac{\partial}{\partial z} \left( \xi_z ~ B_r - \xi_r ~ B_z \right)$$                                 Q_r  =  (A.42)

$$Q_\varphi \frac{\partial}{\partial z} \left( \xi_\varphi ~ B_z - \xi_z ~ B_... ...frac{\partial}{\partial r} \left( \xi_r ~ B_\varphi - \xi_\varphi ~ B_r \right)$$   =  (A.43)

$$\frac{1}{r} ~ \frac{\partial}{\partial r} \left( r ~ ( \xi_z ~ B_... ...partial}{\partial\varphi} \left( \xi_\varphi ~ B_z - \xi_z ~ B_\varphi \right).$$ Q_z  =  (A.44)

Projecting on the 3 axis, and taking the stationary condition into account

$$\rho ~ \left ( r ~ \Omega^2 + g_r \right) = \frac{\partial}{\partial r} ~ \left( p + \frac{B_z^2}{2~\mu_0} \right).$$ (A.45)

The radial component of the Lagrangian displacement satisfies

 

$$\rho ~ \frac{D^2\xi_r}{Dt^2} - 2 ~ \rho ~ \Omega ~ \frac{D\xi_\va... ...al r} + \frac{1}{\rho} ~ \frac{{\rm d}p}{{\rm d}r} ~ {\rm div}~(\rho~\vec{\xi})$$

$$\quad + \rho ~ r ~ \vec{\xi} \cdot \vec{\nabla}~(\Omega^2) - \fra... ...artial r} + \frac{B_\varphi}{r} ~ \frac{\partial Q_r}{\partial \varphi} \right.$$

$$\quad- \left. 2 ~ \frac{B_\varphi~Q_\varphi}{r} + B_z ~ \frac{\pa... ...artial B_r}{\partial \varphi} + Q_z ~ \frac{\partial B_r}{\partial z} \right] =$$

$$\quad ( \rho - {\rm div}~(\rho~\vec{\xi}) ) ~ \delta g_r + \frac{... ...c{\partial}{\partial z} \left( B_r + Q_r + \delta B_*^r \right) \right. \right.$$

$$\quad- \left.\left. \frac{\partial}{\partial r} \left( B_z + Q_z + \delta B_*^z \right) \right\} \right.$$

$$\quad- \left. \frac{\delta B_*^\varphi}{r} ~ \left\{ \frac{\parti... ... \frac{\partial}{\partial\varphi} ( B_r + Q_r + \delta B_*^r ) \right\} \right.$$

$$\quad+ \left. (B_z + Q_z) ~ \left( \frac{\partial\delta B_*^r}{\partial z} - \frac{\partial\delta B_*^z}{\partial r} \right) \right.$$

$$\quad- \left. \frac{B_\varphi + Q_\varphi}{r} ~ \left( \frac{\par... ...*^\varphi) - \frac{\partial\delta B_*^r}{\partial\varphi} \right) \right] \cdot$$ (A.46)

The azimuthal component is given by

 

$$\rho ~ \frac{D^2\xi_\varphi}{Dt^2} + 2 ~ \rho ~ \Omega ~ \frac{D\xi_r}{Dt} - \frac{1}{r} ~ \frac{\partial\Pi}{\partial\varphi}$$

$$- \frac{1}{\mu_0} ~ \left[ B_r ~ \frac{\partial Q_\varphi}{\parti... ... + \frac{B_\varphi~Q_r}{r} +B_z ~ \frac{\partial Q_\varphi}{\partial z} \right.$$

$$+ \left. Q_r ~ \frac{\partial B_\varphi}{\partial r} + \frac{Q_\v... ...\frac{Q_\varphi~B_r}{r} + Q_z ~ \frac{\partial B_\varphi}{\partial z} \right] =$$

$$( \rho - {\rm div}~(\rho~\vec{\xi}) ) ~ \delta g_\varphi + \frac{... ...eft( r ~ ( B_\varphi + Q_\varphi + \delta B_*^\varphi ) \right) \right. \right.$$

$$- \left. \left. \frac{\partial}{\partial\varphi} ( B_r + Q_r + \d... ...} \frac{\partial}{\partial\varphi} ( B_z + Q_z + \delta B_*^z ) \right. \right.$$

$$- \left. \left. \frac{\partial}{\partial z} ( B_\varphi + Q_\varphi + \delta B_*^\varphi )\right\} \right.$$

$$\left. \frac{B_r + Q_r}{r} ~ \left( \frac{\partial}{\partial r}(r... ...lta B_*^\varphi) - \frac{\partial\delta B_*^r}{\partial\varphi} \right) \right.$$

$$- \left. (B_z + Q_z) ~ \left( \frac{1}{r} ~ \frac{\partial\delta ... ...partial\varphi} - \frac{\partial\delta B_*^\varphi}{\partial z} \right) \right]$$ (A.47)

and finally the vertical component satisfies

 

$$\rho ~ \frac{D^2\xi_z}{Dt^2} - \frac{\partial\Pi}{\partial z} - \... ...\partial Q_z}{\partial \varphi} + B_z ~ \frac{\partial Q_z}{\partial z} \right.$$

$$+ \left. Q_r ~ \frac{\partial B_z}{\partial r} + \frac{Q_\varphi}... ...artial B_z}{\partial \varphi} + Q_z ~ \frac{\partial B_z}{\partial z} \right] =$$

$$( \rho - {\rm div}~(\rho~\vec{\xi}) ) ~ \delta g_z - {\rm div}~(\rho~\vec{\xi}) ~ g_z$$

$$+ \frac{1}{\mu_0} ~ \left[ \delta B_*^\varphi \left\{ \frac{1}{r}... ...ial}{\partial z} ( B_\varphi + Q_\varphi + \delta B_*^\varphi )\right\} \right.$$

$$- \delta B_*^r ~ \left\{ \frac{\partial}{\partial z} \left( B_r +... ... - \frac{\partial}{\partial r} \left( B_z + Q_z + \delta B_*^z \right) \right\}$$

$$+ \left. ( B_\varphi + Q_\varphi ) \left( \frac{1}{r} ~ \frac{\pa... ...partial\varphi} - \frac{\partial\delta B_*^\varphi}{\partial z} \right) \right.$$

$$- \left. (B_r + Q_r) \left( \frac{\partial\delta B_*^r}{\partial z} - \frac{\partial\delta B_*^z}{\partial r} \right) \right] \cdot$$ (A.48)

Equations (A.46)-(A.48) are the fundamental equations governing the perturbed motion in the disk due to a small gravitational and magnetic field.

A.4 Orbital and vertical decoupled oscillations

To elucidate the meaning of these equations, we focus on the disk response when motions in vertical direction and orbital plane are decoupled. Let's start with the Lagrangian displacement contained in the orbital plane of the disk by setting $\frac{\partial}{\partial z} = 0$ and $\xi_z=0$. We also assume that perturbations are axisymmetric in order to avoid meaningless complications so that $\frac{\partial}{\partial\varphi}=0$, corresponding to m=0 modes. Then the magnetic perturbations read

                Q_r = 0 (A.49)

$$Q_\varphi$$ = 0 (A.50)

$$- \frac{1}{r} ~ \frac{\partial}{\partial r} \left( r ~ B_z ~ \xi_r \right)$$ Q_z = (A.51)

and the opposite of the total pressure perturbations

$$\Pi = \gamma ~ p ~ \frac{1}{r} ~ \frac{\partial}{\partial r} ... ... \frac{\partial p}{\partial r} - \frac{B_z ~ Q_z}{\mu_0}\cdot$$ (A.52)

Keeping only the leading terms in Eq. (A.47) we find a simple relation between radial and azimuthal displacement, namely

$$\frac{D^2\xi_\varphi}{Dt^2} + 2 ~ \Omega ~ \frac{D\xi_r}{Dt} \approx 0.$$ (A.53)

This can be integrated directly. Dropping the constant of integration we find

$$\frac{D\xi_\varphi}{Dt} \approx - 2 ~ \Omega ~ \xi_r\cdot$$ (A.54)

After straightforward algebra and making use of the equilibrium condition Eq. (6) and the relation (A.54), the radial Lagrangian displacement (A.46) satisfies the following second order partial differential equation

 

$$\rho ~ \frac{D^2\xi_r}{Dt^2} - \frac{\partial}{\partial r} \left[... ...t) ~ \frac{1}{r} ~ \frac{\partial}{\partial r} \left( r ~ \xi_r \right) \right]$$

$$\quad+\rho ~ \kappa_r^2 ~ \xi_r = - \frac{\partial}{\partial r} \... ...( \frac{\delta B_*^z}{\mu_0} ~ \frac{\partial B_z}{\partial r} ~ \xi_r \right).$$ (A.55)

$\kappa_r^2 = 4 ~ \Omega^2 + r ~ \frac{\partial\Omega^2}{\partial r}$ is the radial epicyclic frequency of the oscillations. The sound speed and the Alfven velocity associated with the vertical magnetic field are respectively

$$c_{\rm s}^2 = \frac{\gamma~p}{\rho}$$ (A.56)

$$c_{{\rm a}z}^2 = \frac{B_z^2}{\mu_0~\rho}\cdot$$ (A.57)

In the vertical direction, we neglect the variation in the  $(r,\varphi)$ plane implying that  $\frac{\partial}{\partial r} = \frac{\partial}{\partial\varphi}=0$. We set also  $(\xi_r, \xi_\varphi) = 0$ such that the magnetic perturbations are

$$- \frac{\partial}{\partial z} ~ \left( B_r ~ \xi_z \right)$$                  Q_r = (A.58)

$$Q_\varphi$$ = 0 (A.59)

Q_z = 0 (A.60)

and the opposite of the total pressure perturbations

$$\Pi = \gamma ~ p ~ \frac{\partial\xi_z}{\partial z} + \xi_z ~ \frac{\partial p}{\partial z} - \frac{B_r ~ Q_r}{\mu_0} \cdot$$ (A.61)

After straightforward algebra and making use of the equilibrium condition Eq. (7), the vertical Lagrangian displacement (A.48) satisfies the following second order partial differential equation

 

$$\rho ~ \frac{D^2\xi_z}{Dt^2} - \frac{\partial}{\partial z} \left[... ...0} \right) \frac{\partial\xi_z}{\partial z} \right] + \rho ~ \kappa_z^2 ~ \xi_z$$

$$\frac{\partial}{\partial z} \left( \frac{{\delta B_*^r}^2+{\delta... ...c{\delta B_*^r }{\mu_0} ~ \frac{\partial B_r}{\partial z} ~ \xi_z \right) \cdot$$ (A.62)

The Alfven velocity associated with the radial magnetic field  $c_{{\rm a}r}$is

$$c_{{\rm a}r}^2 = \frac{B_r^2}{\mu_0~\rho}\cdot$$ (A.63)

A.5 Derivation of the inner boundary condition Eq. (51)

In the numerical runs presented in this paper, we use non reflecting boundary conditions which leads to a vanishing density perturbation at the inner edge of the disk, namely  $\delta\rho=0$. In order to compare the linear analysis with simulations, we adopt the same conditions in the both cases. Neglecting the displacement in the vertical direction, density perturbations are related to the Lagrangian displacement by

$$\delta\rho = - {\rm div}~(\rho~\vec{\xi}) = - \frac{1}{r} ~ ... ... ~ \rho ~ \xi_r) - {\rm i} ~ \frac{m}{r} ~ \rho ~ \xi_\varphi.$$ (A.64)

To the lowest approximation given by Eq. (A.54), the radial and azimuthal displacement expanded by Eq. (32) are related by

$$\tilde{\xi}_\varphi \approx - 2 ~ {\rm i} ~ \frac{\Omega}{\omega} ~ \tilde{\xi}_r.$$ (A.65)

Injecting into the density perturbation, we have at the inner edge of the disk r=R₁

$$- \delta\rho = \frac{\partial}{\partial r}( r ~ \rho ~ \til... ...) + 2 ~ m ~ \frac{\Omega}{\omega} ~ \rho ~ \tilde{\xi}_r = 0.$$ (A.66)

We recast the terms in such a way to obtain

$$\frac{\partial\tilde{\xi}_r}{\partial r} + \left( \frac{\part... ...m ~ \frac{\Omega}{r~\omega} ~ \rho \right) ~ \tilde{\xi}_r = 0$$ (A.67)

evaluated at r=R₁. This is the inner boundary condition given by Eq. (51).

   
Appendix B: Approximate solution to the Schrödinger type equation

In this appendix, we recall a method to find approximate solutions to the following Schrödinger type equation with a source term f(x)

 

y''(x) + p(x)   y(x) = f(x). (B.1)

The functions f and p are continuous in the interval [a,b], which includes the origin x=0, (a<0,b>0). Moreover, p has one pole in this interval corresponding to the point x=0 where it changes sign, p(0)=0 and $p'(0)\ne0$.

B.1 Homogeneous solution

We start to look for solutions to the homogeneous part of Eq. (B.1) when f(x)=0. The reference equation is Airy's equation given by

 

w''(x) - x   w(x) = 0. (B.2)

We then look for approximate solutions to Eq. (B.1) which can be cast into the following form, see for example Smirnov (1989)

$$y(x) = A(x) ~ w(\omega_1(x))$$ (B.3)

w is an exact solution of Airy's equation. The functions A and $\omega_1$ are chosen such as to satisfy the homogeneous part of Eq. (B.1) within a prescribed accuracy. The first and second derivatives of y(x) are

 

$$A'(x) ~ w(\omega_1(x)) + A(x) ~ \omega_1'(x) ~ w'(\omega_1(x))$$                                y'(x) = (B.4)

$$A''(x) ~ w(\omega_1(x)) + 2 ~ A'(x) ~ \omega_1'(x) ~ w'(\omega_1(x))$$ y''(x) =

$$+ A(x) ~ \left( \omega_1''(x) ~ w'(\omega_1(x)) (\omega_1')^2(x) ~ w''(\omega_1(x)) \right).$$ (B.5)

However, due to Eq. (B.2), $w''(\omega_1(x)) = \omega_1(x) ~ w(\omega_1(x))$ and therefore the second order derivative can be transform into

 

$$y'' = ( A'' + A ~ (\omega_1')^2 ~ \omega_1 ) ~ w(\omega_1) + ( 2 ~ A' ~ \omega_1' + A ~ \omega_1'' ) ~ w'(\omega_1(x)) .$$ (B.6)

Inserting this expression into Eq. (B.1), (remember that we look for homogeneous solutions with f(x)=0), we get

 

$$( A'' + A ~ (\omega_1')^2 ~ \omega_1 + p ~ A ) ~ w(\omega_1)$$

$$\quad+ ( 2 ~ A' ~ \omega_1' + A ~ \omega_1'' ) ~ w'(\omega_1(x)) = 0 .$$ (B.7)

The Airy functions w and w' being linearly independent, A and $\omega_1$ satisfy the following system

  

$$A'' + A ~ (\omega_1')^2 ~ \omega_1 + p ~ A$$ (B.8)

$$2 ~ A' ~ \omega_1' + A ~ \omega_1''$$ (B.9)

Assuming that the amplitude A is positive, Eq. (B.9) is solved by

 

$$A(x) = \frac{K}{\sqrt{\vert\omega_1'(x)\vert}}$$ (B.10)

where K is constant. Neglecting the second derivative A'' which is one order of magnitude less than the other terms in Eq. (B.8), it reduces to

$$(\omega_1')^2 ~ \omega_1 + p = 0.$$ (B.11)

This relation implies that the function $\omega_1$ has to satisfy the property  $p(x)~\omega_1(x)<0$. We have to distinguish the following two cases.

${ p(x)\ge 0}$

It follows that  $\omega_1(x)\le0$. The solution of Eq. (B.11) is

$$\omega_1(x) = - \left[ \pm \frac{3}{2} ~ \int_0^x \sqrt{p(s)} ~ {\rm d}s \right]^{2/3}.$$ (B.12)

The plus sign applies if ${ p(x)\ge 0}$ for $x\ge0$, otherwise the minus sign applies.

${ p(x)\le 0}$

Then  $\omega_1(x)\ge0$. The solution of Eq. (B.11) is

$$\omega_1(x) = \left[ \pm \frac{3}{2} ~ \int_0^x \sqrt{-p(s)} ~ {\rm d}s \right]^{2/3}.$$ (B.13)

The minus sign applies if ${ p(x)\le 0}$ for $x\le0$, otherwise the plus sign applies.

Note that the function $\omega_1$ is strictly monotonic. $\omega_1'$never changes sign and never vanishes. Consequently, the amplitude Agiven by Eq. (B.10) remains finite in the entire interval [a,b].

B.1.3 Conclusion

Finally, two linearly independent solutions of the homogeneous part of Eq. (B.1) can be expressed with help on the Airy functions of the first and the second kind, respectively Ai and Bi

 

$$\frac{1}{\sqrt{\vert\omega_1'(x)\vert}} ~ {\rm Ai}(\omega_1(x))$$                          y₁(x) = (B.14)

$$\frac{1}{\sqrt{\vert\omega_1'(x)\vert}} ~ {\rm Bi}(\omega_1(x)) .$$ y₂(x) = (B.15)

The associated Wronskian is given by

$$W(y_1,y_2) = y_1 ~ y_2' - y_1' ~ y_2 = \frac{{\rm sign}(\omega_1')}{\pi}\cdot$$ (B.16)

It never vanish, proving that the two solutions remain linearly independent in the entire interval [a,b].

B.2 Inhomogeneous solutions

Knowing the general solutions to the homogeneous part, a peculiar solution to the full Eq. (B.1) with source term is given by the formula

$$y_{\rm p}(x) = \int^x \frac{y_1(x) ~ y_2(t) - y_1(t) ~ y_2(x)}{W(y_1,y_2)} ~ f(t) ~ {\rm d}t.$$ (B.17)

Consequently, the most general approximate analytical solution to our initial problem Eq. (B.1) to the lowest order is

 

$$+ {\rm sign}(\omega_1') ~ \pi \int^x (y_1(x) ~ y_2(t) - y_1(t) ~ y_2(x) ) ~ f(t) ~ {\rm d}t$$ (B.18)

As usual, the constants C₁,C₂ are defined by properly fitting the imposed boundary conditions at x=a and x=b. Applied to the linearized MHD accretion disk, we look for solutions which remains bounded for large radii. This gives the first boundary condition. The second boundary condition is derived by an appropriate density or pressure perturbation at the inner edge. (See end of the appendix)

   
Appendix C: Derivation of the resonance conditions

C.1 Corotation resonance

The corotation resonance is defined by the radius $r_{\rm c}$where  $\omega_*(r_{\rm c})=0$. Actually, this equation possesses two solutions corresponding to  $\omega(r_{\rm c})=\pm\frac{m~c(r_{\rm c})}{r_{\rm c}}$. The width of this region is of the order of the disk height O(H). For very thin disks, this discrepancy can be neglected and the two solutions merge together in an unique corotation radius given by  $\omega(r_{\rm c})=0$. In other words, we assume in this case that  $\omega=\omega_*$. However, in our numerical application, the separation between the two corotation radii is large enough to be resolved. Nevertheless, in the derivation of the corotation resonance below, we will assume that corotation is achieved when  $\omega(r_{\rm c})=m~(\Omega_*-\Omega(r_{\rm c}))=0$. For the detailed study of both corotation, it can be shown that the leading term in the equation satisfied by the radial displacement Eq. (25) is given by

 

$$c_{{\rm ma}z}^2 ~ \left[ 1 + \frac{m^2~c_{{\rm ma}z}^2}{r^2~\omega_*^2} \right] ~ \frac{\partial^2\xi_r}{\partial r^2}$$

$$+ \frac{c_{{\rm ma}z}^2}{r} ~ \left[ \frac{\partial \ln~(r~\rho~c... ...\partial}{\partial r}(\omega_*^{-2}) \right] ~ \frac{\partial\xi_r}{\partial r}$$

$$+ \left[ \omega^2 - \kappa_r^2 + 4 ~ \Omega^2~ \left( 1 - \frac{\... ..._{{\rm ma}z}^2}{r} ~ \frac{\partial}{\partial r}(\omega_*^{-2}) \right] ~ \xi_r$$

$$= \frac{m^2~c_{{\rm ma}z}^2}{\rho~r^2} ~ \frac{B_z~\delta B_*^z}{\mu_0} ~ \frac{\partial}{\partial r}\left( \omega_*^{-2} \right) .$$ (C.1)

Indeed, it is obtained from Eq. (25) by replacing the convective derivative D/Dt by $-i~\omega$, $\xi_\varphi$ by Eq. (37) and neglecting second order terms. Keeping only the leading divergent terms in the coefficients of the ordinary differential Eq. (C.1), we obtain

 

$$\frac{m^2~c_{{\rm ma}z}^4}{r^2~\omega_*^2} ~ \xi_r''(r) + \frac{m... ...4}{r^2} ~ \frac{{\rm d}}{{\rm d}r}\left(\frac{1}{\omega_*^2}\right) ~ \xi_r'(r)$$

$$+ 2~ \frac{m~c_{{\rm ma}z}^2~\Omega~\omega}{r} ~ \frac{{\rm d}}{{... ...*^z}{\mu_0} ~ \frac{{\rm d}}{{\rm d}r}\left(\frac{1}{\omega_*^2}\right)\! \cdot$$ (C.2)

We introduce the new independent variable

$$x = \frac{r-r_{\rm c}}{r_{\rm c}}\cdot$$ (C.3)

Expanding $\omega_*$ to the first order around the corotation radius $r_{\rm c}$ we have

$$\omega_*(r) = \omega_*(r_{\rm c}) + (r-r_{\rm c})~\omega_*'(r_{\r... ... + o(r-r_{\rm c}) = x~r_{\rm c}~\omega_*'(r_{\rm c}) + o(x) \approx \alpha ~ x.$$ (C.4)

To this approximation, we have to solve

 

$$\xi_r''(x) - \frac{2}{x} ~ \xi_r'(x) - 4 ~ \frac{\Omega~\omega~r_... ...c_{{\rm ma}z}^2~x} ~ \frac{B_z(r_{\rm c})~\delta B_*^z(r_{\rm c})}{\mu_0} \cdot$$ (C.5)

This is of the form

$$y''(x) - \frac{2}{x} ~ y'(x) - \frac{b}{x} ~ y(x) = \frac{d}{x}$$ (C.6)

with  $b = 4 ~ \frac{\Omega~\omega~r_{\rm c}^2}{m~c_{{\rm ma}z}^2~x}\ge0$ because $\omega~x = m(\Omega_*-\Omega)~(r-r_{\rm c})/r_{\rm c} \ge0$ and  $d = -2 ~ \frac{r_{\rm c}}{\rho~c_{{\rm ma}z}^2~x} ~ \frac{B_z(r_{\rm c})~\delta B_*^z(r_{\rm c})}{\mu_0}$. Making the change of variable  $t=2~\sqrt{b~x}$ and introducing the new unknown v(t)by  y(t)=t³ v(t), v(t) satisfies the modified Bessel equation of order 3

$$v''(t) + \frac{1}{t} ~ v'(t) - \left( 1 + \frac{9}{t^2} \right) ~ v(t) = 0.$$ (C.7)

This is solved by

 

v(t) = C₁   I₃(t) + C₂   K₃(t). (C.8)

Thus, the complete most general solution to Eq. (C.5) for which a particular solution is easily found to be a constant equal to  $\xi_r^{\rm p}(r) = \frac{m}{2~\Omega~\omega~\rho} ~ \frac{B_z(r_{\rm c})~\delta B_*^z(r_{\rm c})}{\mu_0}$, is

 

$$\xi_r(x) = C_1 ~ x^{3/2} ~ I_3(2~\sqrt{b~x}) + C_2 ~ x^{3/2} ~ K_... ...\Omega~\omega~\rho} ~ \frac{B_z(r_{\rm c})~\delta B_*^z(r_{\rm c})}{\mu_0}\cdot$$ (C.9)

Finally, near the corotation radius, the Lagrangian displacement which remains bounded needs C₁=0

$$\xi_r(x) = C_2 ~ x^{3/2} ~ K_3(2~\sqrt{b~x}) + \frac{m}{2... ...o} ~ \frac{B_z(r_{\rm c})~\delta B_*^z(r_{\rm c})}{\mu_0}\cdot$$ (C.10)

The density disturbance induced in the disk by the rotating gravitational perturbation is then to the lowest order

 

$$\frac{\delta\rho}{\rho} = - \frac{{\rm div}~(\rho~\vec{\xi})}{\rho} = - \frac{1}{\rho~r}~\frac{d}{dr}(r~\rho~\xi_r) - \frac{m}{r~\omega_*^2} ~$$

$$\left( - \frac{m}{r~\rho} ~ \frac{B_z(r_{\rm c})~\delta B_*^z(r_{\rm c})}{\mu_0} + 2~\Omega~\omega~\xi_r \right).$$ (C.11)

The displacement Eq. (C.10) is continuous and differentiable everywhere. Thus, the first term on the right hand side has a finite value. The second term on the RHS needs a special treatment. Indeed, when r approaches $r_{\rm c}$ the numerator $( - \frac{m}{r~\rho} ~ \frac{B_z(r_{\rm c})~\delta B_*^z(r_{\rm c})}{\mu_0} + 2~\Omega~\omega~\xi_r )$ and the denominator  $\omega_*^2$ vanish as well, leading to an undetermined expression of the form 0/0. To find the behavior near $r_{\rm c}$ we introduce the function  $f(r) = - \frac{m}{r~\rho} ~ \frac{B_z(r)~\delta B_*^z(r)}{\mu_0}$. We note that near the corotation, $(f(r) + 2~\Omega~\omega~\xi_r)$ behaves as $f(r)-f(r_{\rm c}) = f'(r_{\rm c})~(r-r_{\rm c})$ with  $f'(r_{\rm c})\ne0$. Thus we conclude that it approaches zero as x. Therefore the density perturbations in the disk diverge as

$$\frac{\delta\rho}{\rho} \approx - ~ \frac{m~f'(r_{\rm c})} {\omega_*'(r_{\rm c})^2 ~ r_{\rm c}^2 ~ x}\cdot$$ (C.12)

The divergent term in the density perturbation Eq. (C.12) tends to infinity as  $\frac{1}{x}$. This result is consistent with the conclusions drawn by Goldreich & Tremaine (1979) for a hydrodynamic disk without self-gravity. We just need to replace the sound speed  $c_{\rm s}^2$ by the fast magnetoacoustic wave speed  $c_{{\rm ma}z}^2$ and the gaseous pressure by the total (gaseous+magnetic) pressure.

C.2 Lindblad and parametric resonance

By choosing an appropriate origin of time t₀, the generalized Mathieu Eq. (66) can be recast in the more convenient form as follows

$$\xi''(t) + \left[ \kappa^2 + h ~ \cos ( \gamma ~ t ) \right] ~ \xi(t) = f ~ \cos( \gamma ~ t )$$ (C.13)

where we have introduced the comoving excitation frequency  $\gamma = \Omega - \Omega_*$. It is a second order linear ordinary differential inhomogeneous equation for the variable $\xi$. The homogeneous part corresponds exactly to Mathieu equation. The related solutions known as Mathieu functions are well studied in mathematical physics textbooks, see for example Morse & Feshbach (1953). These authors show that the solution becomes unbounded when the excitation frequency $\gamma$ lies in a close neighborhood to a multiple of twice the harmonic frequency $\kappa$, namely  $n~\vert\gamma\vert \approx 2~\kappa$ where n is an integer. The size of the neighborhood depends on the amplitude of the modulation h. Because the largest growth rate is reached when the relations holds exactly, we only keep this fastest exponential growing solution. This gives the parametric resonance condition as exposed in Eq. (69).

When the amplitude of the modulation is weak compared to the driven force $h\ll f$, the parametric growth rate can be neglected. Indeed, to this approximation, Eq. (C.13) reduce to the well known driven harmonic oscillator

$$\xi''(t) + \kappa^2 ~ \xi(t) = f ~ \cos( \gamma ~ t ) .$$ (C.14)

Resonance occurs when  $\vert\gamma\vert=\kappa$. In that case the most general solution reads

$$\xi(t) = A ~ \cos(\kappa~t) + B ~ \sin(\kappa~t) + \frac{f~t}{2~\kappa} ~ \sin(\kappa~t).$$ (C.15)

While in the parametric resonance the growth is exponential, here it is only linear with time. Therefore, in the early stage of the full solution of Eq. (C.13), this latter form is a good approximation. On a very long timescale, it converges asymptotically the Mathieu function with exponential growth rate.

Note that the driven resonance (called Lindblad resonance in the main text to refer to the same physical process arising in galactic dynamic and described by Lindblad 1974) is a special case of the parametric resonance for n=2. As a consequence, we could enclose the driven resonance condition into the parametric one. Nevertheless, in order to keep track of the difference between linear and exponential growth rate, this distinction will remain throughout the paper. Moreover, in the numerical application shown in the simulations, the parametric resonance does not appear because of the weak magnetic field disturbance. The timescale of the growing parametric mode excess by several orders of magnitude the duration of the simulations. Only the Lindblad resonance is relevant in this case.