Unstable Disk Accretion to Magnetized Stars: First Global 3D MHD Simulations
Marina M. Romanova, Akshay K. Kulkarni, Richard V.E. Lovelace
aa r X i v : . [ a s t r o - ph ] D ec Draft version November 8, 2018
Preprint typeset using L A TEX style emulateapj
UNSTABLE DISK ACCRETION TO MAGNETIZED STARS: FIRST GLOBAL 3DMHD SIMULATIONS
Marina M. Romanova
Department of Astronomy, Cornell University, Ithaca, NY 14853-6801; [email protected]
Akshay K. Kulkarni
Department of Astronomy, Cornell University, Ithaca, NY 14853-6801; [email protected]
Richard V.E. Lovelace
Department of Astronomy and Applied and Engineering Physics, Cornell University, Ithaca, NY 14853-6801; [email protected]
Subject headings: accretion, dipole — plasmas — magnetic fields — stars
Draft version November 8, 2018
ABSTRACTWe report on the first global three-dimensional (3D) MHD simulations of disk accretion onto a rotatingmagnetized star through the Rayleigh-Taylor instability. The star has a dipole field misaligned relativeto the rotation axis by a small angle Θ. Simulations show that, depending on the accretion rate, a starmay be in the stable or unstable regime of accretion. In the unstable regime, matter penetrates deepinto the magnetosphere through several elongated “tongues” which deposit matter at random placeson the surface of the star, leading to stochastic light-curves. In the stable regime, matter accretes inordered funnel streams and the light-curves are almost periodic. A star may switch between these tworegimes depending on the accretion rate and may thus show alternate episodes of ordered pulsations andstochastic light-curves. In the intermediate regime, both stochastic and ordered pulsations are observed.For Θ > ◦ , the instability is suppressed and stable accretion through funnel streams dominates. INTRODUCTION
Magnetospheric accretion occurs in different stars, forexample in classical T Tauri stars (CTTSs), which are theprogenitors of Solar-type stars (e.g., Hartmann 1998; Bou-vier et al. 2007), in magnetized white dwarfs in some cat-aclysmic variables (e.g., Warner 1995; Warner & Woudt2002), and in accreting millisecond pulsars, which areweakly magnetized neutron stars (e.g., van der Klis 2000).The matter close to the star may either flow around themagnetosphere, forming ordered funnel streams going to-wards the magnetic poles of the star (e.g., Ghosh & Lamb1978; Camenzind 1990; K¨onigl 1991), or may accrete di-rectly through the magnetosphere due to the Rayleigh-Taylor (RT) instability (e.g., Arons & Lea 1976; Elsner &Lamb 1977; Scharlemann 1978). Earlier 2D and 3D nu-merical simulations have shown accretion through funnelstreams (Koldoba et al. 2002; Romanova et al. 2002, 2003,2004). Recent 3D MHD simulations performed for a wider range of parameters have shown that in many cases thedisk-magnetosphere boundary is RT unstable.Previous theoretical models and non-global simulationsgave useful but restricted analyses of this problem. Aronsand Lea (1976) investigated magnetospheric accretionthrough the RT instability in the non-rotating case with aspherical accretion geometry. Spruit and Taam (1990) in-vestigated the stability of an infinitely thin rigidly rotatingdisk. The RT instability in magnetized disks was studiedby Kaisig et al. (1992); Spruit et al. (1995); Lubow &Spruit (1995). Li and Narayan (2004) investigated disk-magnetosphere interaction in the case of an infinitely thickdisk with a vertical magnetic field. 2D simulations havebeen performed by Wang and Nepveu (1983) and Wang& Robertson (1985), while Rast¨atter and Schindler (1999) performed both 2D and 3D simulations, but in a patch(see also Stone and Gardiner 2007). Such simulations,while shedding light on many important features of theinstability, do not take into account many factors whichare present in global simulations, such as the possibilityof matter flowing through funnel streams to the magneticpoles of the star.In this paper we show results from global 3D MHD simu-lations, where the simulation region includes the disk andthe whole magnetosphere of the star. This paper sum-marizes these new results, while in the following paper(Kulkarni & Romanova 2008, hereafter - KR08) we presentresults for a larger range of parameters. MODEL AND RESULTS OF 3D MHD SIMULATIONS . We use a “cubed sphere” Godunov-typenumerical code (Koldoba et al. 2002) and solve the full setof 3D MHD equations with initial and boundary conditionssimilar to those used in Romanova et al. (2004). Simu-lations were done for non-relativistic as well as relativis-tic neutron stars. We approximate relativistic effects us-ing the Paczynski-Wiita potential (Kulkarni & Romanova2005). Our model includes viscosity to regulate the matterflux ˙ M through the disk. The viscous stress in the disk isproportional to the α -parameter, which we varied in therange α = 0 . − .
3. Most of the results shown here arefor the grid N r × N x × N y = 148 × ×
61 in each of the6 blocks of the “cubed sphere”. Test simulations for twiceas fine a grid and twice as coarse a grid show that thenumber of modes does not depend on the grid resolution.The simulations are done in dimensionless form and areapplicable to stars over a wide range of scales, if the mag-netospheric radius r m is not very large compared to theradius of the star R ∗ , r m = (4 − R ∗ ; r m is determined by1 Fig. 1.— (a) An example of stable accretion. The surface is a constant density surface, and the lines are sample magneticfield lines. The magnetic axis is misaligned relative to the rotational axis by Θ = 15 ◦ . (b) An example of unstableaccretion for Θ = 5 ◦ . The lower panels show light curves from the hot spots. Time is measured in orbital periods at r = 1. Only the inner part of the simulation region is shown.the balance between the magnetospheric and matter pres-sure, so that the modified plasma parameter at the disk-magnetosphere boundary β = ( p + ρv ) / ( B / π ) ≈
1. Thereference units are as follows: the length R = R ∗ / . v = ( GM ∗ /R ) , period P = 2 πR /v , den-sity ρ = B /v , magnetic moment µ = B R , accretionrate ˙ M = ρ v R (see KR08 for a complete descriptionof units). We show results for a star with a dimensionlessmagnetic moment µ = 2. We varied the period P ∗ of thestar so that the corotation radius r cor = ( GM/ Ω ∗ ) / indimensionless units varied in the range r cor = 1 . − r cor ≈ α -parameter of viscos-ity. The disk is sufficiently large ( R d ≈ R ∗ ) to supplymatter during the whole simulation run, which lasted upto 50 periods of rotation P . Most of the runs performedwere for the small misalignment angle Θ = 5 ◦ which helpsexcite the perturbations. Runs for larger Θ have shownthat the instability is present up to Θ ≈ ◦ (see KR08). Inour earlier work we performed 3D MHD simulations ofthe disk-magnetosphere interaction for a variety of mis-alignment angles from Θ = 0 ◦ to Θ = 90 ◦ , and showedthat matter flows from the inner regions of the disk to thestar in symmetric funnel streams, which are quasi-stablefeatures over long times (Romanova et al. 2003, 2004).Figure 1a shows an example of accretion through funnelstreams for Θ = 15 ◦ . In this type of flow, matter is liftedabove the equatorial plane and then flows along field lines.The str eams hit the star near the magnetic poles, formingordered hot spots. Next we chose a case with a small misalignment angle(Θ = 5 ◦ ) and increased the accretion rate in the disk, in-creasing the α -parameter. We observed that at a criticalvalue of α = 0 .
04, the instability appears and the accretionacquires an essentially different pattern: the RT instabilitydevelops at the inner edge of the disk, and matter accretesthrough equatorial “tongues” which penetrate deep intothe magnetosphere. The structure of the tongues is oppo-site to that of the funnel streams: they are narrow and tall(see Figure 1b). Such a tongue shape is determined by thegeometry of the magnetic field and has been predicted byArons and Lea (1976). Simulations performed for differentparameters and grid resolutions have shown that the num-ber of growing modes is always small, m = 2 − α = 0 .
02 which provesthat the matter accretion rate, and not α , is important.During the interchange process, the heavy fluid elementsof the disk (matter dominated) change positions with thelight fluid elements of the magnetosphere (magnetically-dominated). The magnetospheric plasma carries frozen-in magnetic field lines, so that displacing the light fluidleads to the effect of pushing magnetic field lines aside(see Figure 2). Equatorial slices of the density distribu-tion for two moments of time are shown in Figure 3 (seealso http://astro.cornell.edu/us-rus/stereo.htm for anima-tions). from the hot spots on the surface ofthe star were calculated, to study the difference in observa-tional properties between the stable and unstable regimes.We assume that the kinetic energy of the stream is con-verted at the surface of the star into isotropic black-bodyradiation. Integration of radiation in the direction of the Fig. 2.—
Penetration of the tongues through the magnetosphere. The tongues are shown by density contours in theequatorial plane (the highest density is shown in red, and the lowest in dark-blue, with a density contrast of about 300).
Fig. 3.—
Density distribution in the equatorial plane after 20 and 30 periods of rotation at r = 1. The background showsthe density (the colors have the same meaning as in Figure 2). The black line shows the position of β = 1. The plots areshown in the reference frame rotating with the star.observer is performed (see details in R04). In the stableregime, the position of the hot spots is almost constant,so that the light-curve is almost sinusoidal (see Figure 1a,bottom panel). In the strongly unstable regime, sporadi-cally forming tongues hit the star in random places, so thatthe light curve is irregular (see Figure 1b, bottom panel).In the intermediate cases both funnels and tongues arepresent, and mixed light-curves are expected. If a certainnumber of tongues dominates (for example, the m = 2mode dominates in some cases) then quasi-periodic os-cillations are expected. A more detailed analysis of thelight-curves in different cases will be presented in a futurepaper. . We compared our simula-tions with a few relevant theoretical approaches. A heavyfluid supported against gravity by a lighter one is RT un-stable unless there is some opposing force. A homoge-neous vertical field at the disk-magnetosphere boundarydoes not oppose the growth of φ -modes (see e.g. Chan- drasekhar 1961). A more general criterion for differen-tially rotating magnetized accretion disks indicates thatthe disk is unstable to excitation of φ -modes if γ B Σ ≡− g eff d ln(Σ /B z ) /dr > r m d Ω /dr ) ≡ γ (Spruit et al.1995; Spruit & Taam 1990; Kaisig et al. 1992; Lubow &Spruit 1995). Here, − g eff = r (Ω K − Ω ) is the effectivegravity, and Σ = 2 ρh is the surface density. That is, for theinstability to start, the surface density per unit magneticfield strength Σ /B z should drop off fast enough in the di-rection of the effective gravity ( − g eff ), so that the term γ B Σ is larger than the term associated with the shear, γ ,which tends to oppose the instability by smearing out theperturbations. We calculated the φ -averaged values of γ B Σ and γ for the stable and unstable cases (see Figure 4) andcompared the positions of the curves in the inner regionsof the disk where the instability starts, that is at r ∼ > r m .One can see that in the stable case ( α = 0 . γ B Σ ∼ < γ while in the unstable case ( α = 0 . γ B Σ >> γ . Thus thetheory makes a reasonably good prediction for the onset of (α=0.1) γ ΒΣ unstable r γ Ω g eff r m (α=0.02) g eff stable r γ ΒΣ γ Ω r m Fig. 4.—
Parameters corresponding to the stability criterion by Spruit et al. (1995) described in § (α=0.02) stable ρΩΩ Κ r m r (α=0.1) unstable ρΩΩ Κ r m r Fig. 5.—
Radial distribution of the density ρ , angular velocity Ω and Keplerian angular velocity Ω K in the equatorialplane. The dashed line r m marks the position of the magnetospheric radius where β = 1.instability. The term γ B Σ varies strongly mainly becausein the unstable case, a significant amount of matter accu-mulates near r m and the density drop-off is sharper thanin the stable case (see Figure 5).The RT modes observed in simulations are excited bydifferent mechanisms. One of them is the small non-axisymmetry associated with the misalignment of thedipole. Another is the inhomogeneities in the inner regionsof the disk driven by the azimuthal component of the fieldtrapped inside the inner regions of the disk. The Kelvin-Helmhotz instability seems to be less significant (see alsoRast¨atter & Schindler 1999). . According to thetheory of the RT instability, high- m modes should growfaster than low- m ones (e.g. Chandrasekhar 1961). In arotating disk, the shear d Ω /dr may efficiently damp high- m modes (e.g., Lubow & Spruit 1995) and may be themain reason for damping of the high- m modes in our sim-ulations. Initially, in high-resolution simulations we seethe formation of m = 20 −
30 modes, but in one rotationperiod they are damped, and only the low- m modes grow.KR08 analyzed the number of modes based on the Li &Narayan (2004) criterion and reached a similar conclusion.Another possible reason for the suppression of the high- m modes is the presence of the azimuthal component B φ ofthe field at the inner edge of the disk which typically con-stitutes ∼ (5 − B φ is expected from theory (Chandrasekhar1961), and has been noticed by Wang & Nepveu (1983)and Rast¨atter & Schindler (1999). . We performed multiple simulation runs for iden-tical initial density distributions in the disk but different viscosity parameters α and dimensionless periods of thestar P dim ∗ (these bulk simulations were done for a grid72 x ). The symbols in Figure 6a correspond to differ-ent runs. One can see that there is a boundary betweenthe stable and unstable cases with unstable cases corre-sponding to higher α . For similar initial disks and ac-cretion rates through the disks, ˙ M is proportional to α ,and thus the instability occurs at high enough accretionrates through the disk. The boundary is higher at smallerperiods P dim ∗ , because at faster rotation of the star the ef-fective gravity becomes less negative and thus less mattermay accumulate at the inner edge of the disk. That is whya larger α is required to bring more matter to the inneredge of the disk to satisfy the instability criterion. Caseswith periods P dim ∗ ≈ . − P dim ∗ ∼ < .
8, enterthe “propeller” regime, where r cor < r m (e.g. Illarionov& Sunyaev 1975; Lovelace et al. 1999; Rappaport et al.2004). In this regime a star tends to expel a significantpart of incoming matter (Romanova et al. 2005; Ustyu-gova et al. 2006). At high enough α we observed events ofinstability even in this regime (see also Wang & Robertson1985), but the simulation runs were brief, and this regionrequires further analysis.We calculated the accretion rate onto the surface of thestar ˙ M dim ∗ for the above runs. The boundary between sta-ble and unstable accretion ˙ M dim ∗ ≈ . P ∗ , which is probably a result of the combi-nation of the lower disk density due to the propeller effect,and higher α . This boundary is approximate, with inter-mediate cases above and below it in which only marginal unstablemarginalstableboundary α stableunstable r cor P * dim -1 M * p r ope ll e r o b stableunstable r = r cor m dim r cor a P * dim p r ope ll e r Fig. 6.—
Position of stable (circles), unstable (triangles) and marginal (squares) cases in the (a) α - P dim ∗ and (b) ˙ M dim ∗ - P dim ∗ planes. The top horizontal axis shows the corotation radius. The solid lines show an approximate boundary betweenstable (below) and unstable (above) runs. The dash-dotted line marks the approximate boundary of the propeller regime.instability is observed. In this transition region, both fun-nels and tongues are present. Higher up, for ˙ M dim ∗ > . The results of our simulations can be applied to a vari-ety of magnetized stars with small magnetospheres, r m ≈ (4 − R ∗ , including CTTSs, weakly magnetized whitedwarfs (WDs) in cataclysmic variables and neutron stars(NSs) in LMXBs. The dividing line in Figure 6b corre-sponds to ˙ M dim ∗ ≈ .
3. Below we show the accretion ratecorresponding to this boundary in dimensional units fordifferent types of stars:CTTSs : ˙ M cr ∗ ≈ . × − ( B ∗ / G) r m − M ⊙ / yr , WDs : ˙ M cr ∗ ≈ . × − ( B ∗ / G) r m − M ⊙ / yr , NSs : ˙ M cr ∗ ≈ . × − ( B ∗ / G) r m − M ⊙ / yr . The period in dimensional units is: P ∗ = A p P dim ∗ m − r , where scaling factors for different stars are: CTTSs: A p ≈ .
78 days, r = R ∗ / R ⊙ , m = M ∗ / . M ⊙ , WDs: A p ≈ . r = R ∗ / × cm, m = M ∗ / M ⊙ , NSs: A p ≈ .
22 ms, r = R ∗ / cm, m = M ∗ / . M ⊙ .Note that the dimensional ˙ M cr ∗ increases with the mag-netic field as B ∗ so that the boundary will be higher forhigher B ∗ . This reflects the fact that we have fixed the di-mensionless magnetic moment µ , and thus consider starswith an approximately fixed ratio r m /R ∗ ≈ −
5. Weperformed another set of runs for smaller magnetospheres(smaller r m /R ∗ for µ = 0 .
5) and observed that the di-mensionless boundary corresponds to a smaller ˙ M dim ∗ (see KR08). Note also that all the above runs were done forΘ = 5 ◦ . We expect that for larger Θ the boundary willmove up, because funnel accretion is more favorable atlarger Θ. OBSERVATIONAL CONSEQUENCES
The existence of two regimes of accretion changes ourunderstanding of accreting magnetized stars and their ob-servational properties: the lack of a periodic signal in thelight curve does not rule out accretion and an ordered mag-netic field. For example, the random light-curves in manyCTTSs or dwarf novae do not preclude an ordered stellarfield, if the star is in the unstable regime. Depending onthe accretion rate, episodes of stable and unstable accre-tion may alternate, so that periods with pulsations maybe intermittent and may be followed by periods with nopulsations. Recently a number of intermittent accretingmillisecond pulsars have been discovered (e.g., Kaaret etal. 2006; Altamirano et al. 2007; Galloway et al. 2007;Gavriil et al. 2007). Although the reason for this behavioris not yet understood, we note that in, for example, HETEJ1900.1–2455, the pulsations disappeared when the overallX-ray flux increased (Kaaret et al. 2006) which may be anexample of a transition to the unstable accretion regime.In some cases a definite number of tongues may dom-inate, which may lead to quasi-periodic oscillations inthe light curves (Li & Narayan 2004). These oscilla-tions may be candidates for one of the QPOs observedin Type II (accretion-driven) bursts in LMXBs. Comp-tonization of photons by high-energy electrons may lead toonly a small departure of the light-curve from the thermalones obtained using the approximation of isotropic black-body radiation (Poutanen & Gierli´nski 2003), so that theQPO features may survive Comptonization (see, however,Titarchuk et al. 2007).In the case of young stars surrounded by gas/dust diskswhere planets are forming and migrating inward, the un-stable tongues may support inward migration, so that themagnetospheric gap, which can halt migration (e.g. Lin etal. 1996; Romanova & Lovelace 2006; Papaloizou 2007),may form only in the state of stable accretion.We thank the anonymous referee for valuable sugges-tions which improved the paper. We thank Drs. A.Koldoba and G. Ustyugova for their contribution to thecode development and Drs. D. Altamirano, L. Hillen-brand, P. Kaaret, and C. Thompson for discussions.NASA provided high-performance computational facilities for this work. The research was partially supported bythe NSF grants AST-0507760 and AST-0607135, and theNASA grants NNG05GG77G and NNG05GL49G.and may be followed by periods with nopulsations. Recently a number of intermittent accretingmillisecond pulsars have been discovered (e.g., Kaaret etal. 2006; Altamirano et al. 2007; Galloway et al. 2007;Gavriil et al. 2007). Although the reason for this behavioris not yet understood, we note that in, for example, HETEJ1900.1–2455, the pulsations disappeared when the overallX-ray flux increased (Kaaret et al. 2006) which may be anexample of a transition to the unstable accretion regime.In some cases a definite number of tongues may dom-inate, which may lead to quasi-periodic oscillations inthe light curves (Li & Narayan 2004). These oscilla-tions may be candidates for one of the QPOs observedin Type II (accretion-driven) bursts in LMXBs. Comp-tonization of photons by high-energy electrons may lead toonly a small departure of the light-curve from the thermalones obtained using the approximation of isotropic black-body radiation (Poutanen & Gierli´nski 2003), so that theQPO features may survive Comptonization (see, however,Titarchuk et al. 2007).In the case of young stars surrounded by gas/dust diskswhere planets are forming and migrating inward, the un-stable tongues may support inward migration, so that themagnetospheric gap, which can halt migration (e.g. Lin etal. 1996; Romanova & Lovelace 2006; Papaloizou 2007),may form only in the state of stable accretion.We thank the anonymous referee for valuable sugges-tions which improved the paper. We thank Drs. A.Koldoba and G. Ustyugova for their contribution to thecode development and Drs. D. Altamirano, L. Hillen-brand, P. Kaaret, and C. Thompson for discussions.NASA provided high-performance computational facilities for this work. The research was partially supported bythe NSF grants AST-0507760 and AST-0607135, and theNASA grants NNG05GG77G and NNG05GL49G.