Accretion to Magnetized Stars through the Rayleigh-Taylor Instability: Global Three-Dimensional Simulations

We present results of 3D simulations of MHD instabilities at the accretion disk-magnetosphere boundary. The instability is Rayleigh-Taylor, and develops for a fairly broad range of accretion rates and stellar rotation rates and magnetic fields. It manifests itself in the form of tall, thin tongues of plasma that penetrate the magnetosphere in the equatorial plane. The shape and number of the tongues changes with time on the inner-disk dynamical timescale. In contrast with funnel flows, which deposit matter mainly in the polar region, the tongues deposit matter much closer to the stellar equator. The instability appears for relatively small misalignment angles, $\Theta\lesssim30^\circ$, between the star's rotation and magnetic axes, and is associated with higher accretion rates. The hot spots and light curves during accretion through instability are generally much more chaotic than during stable accretion. The unstable state of accretion has possible implications for quasi-periodic oscillations and intermittent pulsations from accreting systems, as well as planet migration.


INTRODUCTION
Magnetospheric accretion occurs in different systems, including classical T Tauri stars (CTTSs), which are the progenitors of solar-type stars (e.g., Hartmann 1998;Bouvier et al. 2007), in magnetized cataclysmic variables, which are accreting white dwarfs (e.g., Warner 1995;Warner & Woudt 2002), and in millisecond pulsars, which are weakly magnetized accreting neutron stars (e.g., van der Klis 2000). The geometry of the accretion flow around magnetized stars is a problem of long-standing interest. It is an important factor in determining the observed spectral and variability properties of the accreting system. An accretion disk around a magnetized central object is truncated at the distance from the central star where the magnetic energy density becomes comparable to the matter energy density. Beyond that point, there are two ways in which the gas can accrete to the star: (1) Through funnel streams, or magnetospheric accretion (e.g., Ghosh & Lamb 1978, 1979; (2) Through plasma instabilities at the disk-magnetosphere interface (e.g., Arons & Lea 1976;Elsner & Lamb 1977). The geometry of the matter flow in these two regimes is expected to be very different. In general, the Rayleigh-Taylor (RT), or interchange, instability is expected to develop at the disk-magnetosphere interface because of E-mail: akshay@astro.cornell.edu † E-mail: romanova@astro.cornell.edu the high-density disk matter being supported against gravity by the low-density magnetospheric plasma. The Kelvin-Helmholtz (KH) instability is also expected to develop because of the discontinuity in the angular velocity of matter at the boundary. The inner disk matter is expected to rotate at the local Keplerian velocity, while the magnetospheric plasma corotates with the star.
A number of theoretical analyses of such instabilities have been performed. Arons & Lea (1976a,b) presented a detailed analysis of magnetospheric accretion through the RT instability for a non-rotating star and a spherical accretion geometry. Anzer (1969) and Anzer & Börner (1980 have analyzed the RT instability for solar prominences and the KH instability at the inner disk edge respectively. Baan (1977Baan ( , 1979 examined the role of the RT instability in X-ray bursts. Lovelace & Scott (1981) investigated the interchange instability in a two-dimensional accretion disk threaded by a vertical magnetic field. Spruit & Taam (1990) investigated the stability of a thin magnetized disk with a rigid rotation profile. They present a solution to one of the problems with the model of Ghosh & Lamb (1979), namely that the incoming disk squeezes the magnetosphere in the equatorial region, and it is energetically impossible for the inner disk matter to follow the resulting magnetospheric field lines. The solution is that the matter at the inner disk edge can be unstable, causing it to move towards the star across the field lines, until it reaches field lines that are energetically possible to follow.
The RT instability in differentially rotating magnetized disks was studied by Spruit, Stehle & Papaloizou (1995) and Lubow & Spruit (1995). Li & Narayan (2004) analyzed RT and KH instabilities at the disk-magnetosphere interface for an infinitely thick disk and a vertical magnetic field.
Theoretical analysis, although useful in understanding the basic features of the instabilities, is limited to the linear regime. To understand the behaviour of the instabilities in the nonlinear regime, numerical simulations are needed. While there have been some numerical studies of such instabilities, they have so far been limited in one way or another, and have mostly focused on the physical characteristics of the instabilities themselves, with relatively less attention being given to their effect on the nature of the accretion flow, and to their dependence on the physical parameters of accreting systems, like the accretion rate and the star's rotation rate and magnetic field. Moreover, realistic accretion geometries can only be obtained from global 3D simulations, which have so far been lacking. 2D simulations have been performed by Wang & Nepveu (1983), Wang & Robertson (1984, 1985 and Kaisig, Tajima & Lovelace (1992). Rastätter & Schindler (1999b) performed 3D simulations in a patch which includes a part of the magnetospheric boundary and the equatorial plane, using semianalytically derived initial conditions (Rastätter & Neukirch 1997;Rastätter & Schindler 1999a). Stone & Gardiner (2007a,b) performed 3D simulations of the magnetic RT instability in a shearing box for the idealized case of two inviscid, perfectly conducting fluids of constant density separated by a contact discontinuity perpendicular to the effective gravity, with a uniform magnetic field parallel to the interface.
Earlier, two-and three-dimensional simulations have shown accretion through funnel streams (Romanova et al. , 2003Kulkarni & Romanova 2005). In this paper we report on accretion through the Rayleigh-Taylor, or interchange, instability in global 3D MHD simulations, where the simulation region includes the disk and the whole magnetosphere of the star. This paper follows up on an earlier, briefer paper (Romanova, Kulkarni & Lovelace 2008), giving a more detailed description of the problem, and showing results for a larger range of star and disk parameters.

THE NUMERICAL MODEL
The model we use is the same as in our earlier 3D MHD simulations Romanova et al. 2003Romanova et al. , 2004. The star has a dipole magnetic field, the axis of which makes an angle Θ with the star's rotation axis. The rotation axes of the star and the accretion disk are aligned. There is a low-density corona around the star and the disk which also rotates about the same axis. To model stationary accretion, the disk is chosen to initially be in a quasi-equilibrium state, where the gravitational, centrifugal and pressure gradient forces are in balance . Simulations were done for both relativistic and non-relativistic stars. General relativistic effects, which are important for neutron stars, are modelled using the Paczyński-Wiita potential Φ(r) = GM * /(r − r g ) (Paczyński & Wiita 1980), where r g ≡ 2GM * /c 2 is the Schwarzschild radius of the star. This potential reproduces some important features of the Schwarzschild spacetime, such as the positions of the innermost stable and marginally bound circular orbits. Viscosity is modelled using the α-model (Shakura & Sunyaev 1973;Novikov & Thorne 1973), and controls the accretion rate through the disk. To model accretion, the ideal MHD equations are solved numerically in three dimensions, using a Godunov-type numerical code, written in a "cubed-sphere" coordinate system rotating with the star Romanova et al. 2003). The boundary conditions at the star's surface amount to assuming that the infalling matter passes through the surface of the star. So the dynamics of the matter after it falls on the star is ignored. The inward motion of the accretion disk is found to be stopped by the star's magnetosphere at the Alfvén or magnetospheric radius, where the magnetic and matter energy densities become equal. The subsequent evolution depends on the parameters of the model.

Reference Values
The simulations are done using dimensionless variables. For every physical quantity q, the dimensionless value is defined as q = q/q 0 , where q 0 is the reference value for q. Because of the use of dimensionless variables, the results are applicable to a wide range of objects and physical conditions, each with its own set of reference values, provided the magnetospheric radius r m is not very large compared with the radius of the star R * , r m = (4 − 6)R * ; r m is determined by the balance of the magnetospheric and matter pressure, so that the modified plasma parameter at the disk-magnetosphere boundary β = (p + ρv 2 )/(B 2 /8π) ≈ 1. To apply our simulation results to a particular situation, we have the freedom to choose three parameters, and all the reference values are calculated from those. We choose the mass, radius and equatorial surface magnetic field of the star as the three independent parameters. Appendix A shows how the reference values are determined, and lists the reference values for three classes of central objects -classical T Tauri stars, white dwarfs and neutron stars.
Subsequently, we drop the primes on the dimensionless variables and show dimensionless values in the figures.

SIMULATION RESULTS
We chose the following parameters for our main case: dipole moment µ = 2, corresponding to an equatorial surface magnetic field of B * = 2, misalignment angle Θ = 5 • , viscosity parameter α = 0.1, stellar rotation period P = 3, initial disk radius = 2. The corotation radius, which is the radius at which the orbital rotation rate of the disk matter equals the star's rotation rate, is r cor = 2. The cubed-sphere grid resolution in most of our cases is N r × N 2 = 72 × 31 2 in each of the six blocks of the grid (which is the resolution in all the figures in this paper, unless otherwise noted), and we have some runs at higher resolutions (100 × 41 2 , 144 × 61 2 , 216 × 91 2 ). The outer radius of our simulation region is r out ≈ 12 ≈ 35R * in most runs, although we only show the region near the star in our plots. Fig. 1 shows two views of the accretion flow around the star through instabilities (top row) and two views of a magnetospheric (funnel) flow (bottom row). The growth of unstable perturbations at the disk-magnetosphere boundary results in penetration of the magnetosphere by the disk matter, in the form of tongues of gas travelling through the equatorial  Left panel: A tongue of gas, shown by density contours in the equatorial plane, pushing aside magnetic field lines on its way to the star. Note that the "hole" in the magnetosphere is not artificially depicted -the field lines start out uniformly spaced on the star's surface, and twist aside around the tongue. plane. Matter energy density dominates inside the tongues. When they come closer to the star, they encounter a stronger magnetic field, which stops their equatorial motion. At this point the tongues turn into miniature funnel-like flows following the field lines, which gives them a characteristic wishbone shape. They deposit matter much closer to the star's equator than true funnel flows do.
The tongues are tall and thin, as opposed to the funnels which are flat and wide. This is because the tongues pene-trate the magnetosphere by prying the field lines aside (Fig.  2), since this is energetically more favourable than bending the field lines inward. This is a standard feature of the interchange instability -if a heavy fluid is supported against gravity by a light fluid, then at the boundary between the fluids, fluid elements on either side displace those on the other side. Since the magnetic field is frozen into the matter, the field lines are also pushed aside in the process. This interchange process continues, producing fingers of each fluid penetrating into the other (see, e.g., Arons & Lea 1976a, Wang & Robertson 1985. The tension of the magnetic field lines suppresses the interchange process in the direction parallel to the field. This increases the characteristic wavelength of deformation of the boundary in that direction (Chandrasekhar 1961). Hence, the tongues are narrow in the direction perpendicular to the field, i.e., the azimuthal direction, and broader parallel to the field, i.e., the vertical direction (Stone & Gardiner 2007a,b). Fig. 2 shows that the magnetosphere of the star is strongly disturbed by the penetrating tongues which push the field lines aside. To conserve angular momentum, the angular velocity of the gas in the tongues increases as it moves inwards, making it rotate faster than the inner disk matter, causing the tongues to curve to their right. This tongue shape allows tongues that move faster to merge with more slowly moving tongues. As a result, the total number of tongues at any given time is of the order of a few. Such merging of Rayleigh-Taylor fingers has been observed in earlier numerical simulations (Wang & Robertson 1985). The merging and subsequent growth of the tongues occurs on the inner-disk dynamical timescale. The number of tongues we see is of the order of a few.
The matter in the tongues is accelerated by the gravitational field of the star, so that its velocity is higher near the star. The density and velocity distribution in the tongues is similar to that in the funnel streams, as Fig. 3 shows. It is also seen that the funnels and tongues are not mutually exclusive. We discuss this last point in more detail in §3.1 and §4. Fig. 4 shows equatorial slices of the circumstellar region at various times. The density enhancements which result in the formation of the tongues can be seen at the bases of the tongues. The number of tongues changes with time of its own accord, without any artifically introduced perturbation.
To estimate the number and positions of the tongues at different times, we plot the plasma density in the equatorial plane along circles of radius R centered at the star, as a function of time (Fig. 5). The density enhancements show the po-sitions of the tongues. The dashed vertical lines show how many tongues cross a certain radius at a specific time. The top panel shows the evolution of the tongues at R = 1.4. This is close to the inner radius of the disk just before the start of the instability. This is the place where the number of tongues is the largest (usually 4-7). The tongues are dense at this radius. The middle panel shows that closer to the star, at R = 1 there are usually 3-4 tongues, and their density is smaller. Even closer to the star, at R = 0.6 (lower panel), the tongues become weaker and their number decreases to 1-3. This is because close to the star, many of the tongues leave the equatorial plane and travel along the magnetic field lines to the star.
This plot also shows us the azimuthal positions φ of the tongues as a function of time, giving a rough idea of their angular velocity. The azimuthal position φ is measured in the coordinate system rotating with the star. We can see that the tongues move and change shape on the dynamical timescale in the inner disk region. The number of tongues varies between about 2 and 7. Fig. 6 shows the hot spots on the star's surface for our main case at different times. We see that the spots are different from pure funnel-flow hot spots , Kulkarni & Romanova 2005, and are significantly different from the simple polar-cap shape that is frequently assumed. Each tongue creates its own hot spots when it reaches the star's surface. Therefore, the shape, intensity, number and position of the spots change on the inner-disk dynamical timescale. As noted earlier in this section, since the density and velocity of the matter in the tongues is comparable to that in the funnel streams ( Fig. 3), the energy flux from the hot spots created by the tongues is also comparable to that from funnel-stream hotspots in other, stable cases. Fig. 7 compares the typical lightcurves during accretion through funnels and through instability. The lightcurve is usually very chaotic during accretion through instability, and shows no definite periodicity. Sometimes, however, a certain number of tongues dominates, and the lightcurve may show quasi-periodicity. Also, as mentioned above, funnels and tongues coexist for certain ranges of parameter values. In this intermediate regime, a mixed lightcurve with both periodic and chaotic components is expected. We plan to discuss this in more detail in a future paper.
In the following subsections, we investigate the dependence of the instability on different parameters, varying one parameter at a time and choosing the same values for the other parameters as in the main case.

Dependence on the Accretion Rate
The accretion rateṀ in our code is regulated by the viscosity coefficient ν which is proportional to the α-parameter (Shakura & Sunyaev 1973). The radial velocity of inward flow in the disk at a distance r from the star is v r ∼ ν/r ∼ αc s h/r, where h is the thickness of the disk and c s is the sound speed. Thus the accretion rate through the disk is approximately proportional to α:Ṁ ≈ 4πrhρv r ∼ α. The accretion rate is also proportional to the fiducial densityρ in the disk, but we keep ρ fixed in all our runs (except in test runs that we describe at the end of this subsection), giving us the same initial density distribution in all runs.
We performed simulations for a wide range of the  parameter α, from very small to relatively large, α = 0.02, 0.03, 0.04, 0.06, 0.08, 0.1, 0.2, 0.3. We find that at very small α ( 0.03), the instability does not appear. At larger α's, the instability appears, and when α is increased, the instability starts earlier and more matter accretes through it. Fig. 8 shows equatorial slices of the plasma density distribution at different α. One can see that there are no tongues at α = 0.02. The tongues are quite weak at α = 0.04, but much stronger at larger α, when more matter comes to the inner region of the disk, and the plasma density in the inner region of the disk is higher than in the low-α cases. This shows that increased accumulation of mass at the inner edge of the disk leads to enhancement of the instability, producing tongues that propagate deeper into the magnetosphere of the star. We should note that in spite of different conditions at the inner region of the disk (much higher density at larger α), the number  and behaviour of the tongues is approximately the same in all cases. Fig. 9 shows the density distribution in the µ − Ω plane. One can see that at α = 0.02, the accretion is entirely through magnetospheric funnel streams. At α = 0.04 and α = 0.06, a significant amount of matter accretes through the funnel streams, though some accretes through instabilities. We call this the intermediate regime of accretion. The bottom row shows the most unstable cases, where most of the matter flows through the tongues. For the unstable cases, what we see inside the magnetosphere in this figure is cross sections of the tongues.
The accretion rate onto the star's surface is higher during accretion through instability, as Fig. 10 shows. We see that the accretion rate increases with increasing α. This is mainly due to increase in the amount of matter transported inwards by the accretion disk. The higher accretion rate is accompanied by a higher angular momentum flux (Fig. 11).
The position and number of the hot spots on the star's surface varies with time. To find the average location of the spots, we calculated the matter flux integrated over magnetic longitude,Ṁ(θ m ), as a function of magnetic latitude θ m .Ṁ(θ m ) is defined such that the total accretion rate onto the star,Ṁ = Ṁ (θ m )dθ m . Fig. 12 shows the evolution ofṀ(θ m ) with time. The top panel shows that at large α, when the instability is strong, the hot spots are located at the mid-latitudes, being brightest at 35 • < θ m < 65 • . The bottom panel shows that at very small α, when matter ac-  cretes through magnetospheric funnels, the spots are located at much higher latitudes, 60 • < θ m < 75 • . In the intermediate case, α = 0.04 (middle panel) both types of accretion are present and the plots reflect hot spots produced both by magnetospheric streams and by instabilities. Fig. 13 shows the θ m -dependence ofṀ(θ m ) at t = 25 for different α. The plot shows that at α = 0.02 and α = 0.03 the maximum of matter flux is located at θ m ≈ 70 • with the half-width of the peak ≈ 75 It is surprising to see that at the largest viscosity coefficients, α = 0.2 and 0.3, the hot spots do not move closer to equator, but have a maximum at θ m ≈ 50 • , like for α = 0.1.
To eliminate the possibility that α-viscosity is directly responsible for the instability, we performed test simulation runs atρ = 2 and α = 0.02, and found that the instability develops in this case, although, as noted earlier, the run atρ = 1 and α = 0.02 is stable. Similarly, in a run withρ = 0.5 and α = 0.04, no instability is observed, although forρ = 1 and α = 0.04, the instability exists. This and a few other test runs have shown that it is really the matter flux, and not α itself, that is responsible for the onset of the instability.
This subsection shows that the accretion rate (controlled by α) is one of the important factors determining whether the matter flow is stable (with magnetospheric funnel flows) or unstable, where most of the matter flows through equatorial instabilities.

Dependence on the magnetic field strength and the star's rotation rate
In qualitative terms, the instability depends on the effective gravitational force (i.e., the difference between the gravitational and centrifugal forces) acting on the inner disk matter. The star's magnetosphere forces the matter at the inner disk boundary to approximately corotate with the star. So if the magnetospheric radius r m is significantly smaller than the corotation radius r cor , then the inner disk matter is appreciably slowed down by the magnetosphere, and the effective gravitational force is large. The effect of this most clearly seen when comparing runs with different magnetic fields or with different stellar rotation rates. First, starting from our main case, with parameters α = 0.1,ρ = 1, P = 3, r cor = 2 and Θ = 5 • , we varied the magnetic moment µ of the star. We observed that the instability appears at a wide range of µ, from µ = 0.2 to µ = 8. Fig. 14 shows that the small µ, r m is very small, and also the typical mode number m = 2. At very large µ (µ = 8), r m ≈ 1.8 − 2.2 becomes comparable to r cor = 2, and the star is in the propeller regime, so that the effective gravity is very small or even radially outwards, and is unable to drive the instability. Thus, increasing the magnetic field leads to suppression of the instability by decreasing the effective gravity.
In another set of runs, we started from our main case and fixed µ = 2 (which fixes r m ), but varied the rotation period P of the star. We observed that at smaller P, the instability becomes weaker and is finally suppressed, again because decreasing the period decreases r cor , bringing it closer to r m , taking the star into the propeller regime and reducing the effective gravity. We discuss the dependence of the instability on the star's magnetic field and rotation rate in more detail in §4. Fig. 15 shows the equatorial matter flow at various misalignment angles Θ, all other parameters being the same as in our main case. We find that the instability shuts off for Θ 30 • . The reason for this is that for large misalignment angles (Θ 30 • ), the magnetic poles are closer to the disk plane. Therefore, the gravitational energy barrier that the gas in the inner disk region has to overcome in order to form funnel flows is reduced, making funnel flows energetically more favourable. The matter seen inside the magnetosphere for Θ = 60 • is part of the warped funnel stream that crosses the disk plane. We also find that for Θ = 30 • , the m = 2 mode (two tongues) usually dominates. Fig. 16 shows the accretion rateṀ(θ) onto the star's surface as a function of rotational latitude θ. We see that when the accretion is through instability (Θ 30 • ), most of the matter accretes onto the mid-latitude (θ ∼ 50 • ) region of the star, independent of the star's misalignment angle.

Dependence on the grid resolution
We found that the azimuthal extent of each tongue is much larger than the size of our grid cells. This indicates that the instability is not an artefact of the coarseness of the grid. Nevertheless, to eliminate that possibility, we performed test simulations at higher grid resolutions. Fig. 17 shows equatorial slices of the region near the star at various grid resolutions. The instability exists at all the resolutions we tested, including a test run at resolution 216×91 2 (not shown here). The number and behaviour of the tongues is similar in all these cases. However, with the finer grid, the tongues are thinner on an average. The accretion rate onto the star at the finest grid is about 20% larger than that with the coarsest grid, which may be related to the smaller numerical diffusivity at finer grids, and the resultant better coupling between the outer magnetosphere and the disk, which leads to a higher accretion rate.

Possible perturbation mechanisms
The Rayleigh-Taylor or Kelvin-Helmholtz instabilities will only manifest themselves in a potentially unstable layer between two liquids if some perturbation of density, pressure or velocity occurs in the layer. There are different mechanisms of perturbation which are expected in real accretion disks. First of all, the matter in the disk is never perfectly homogeneous, and natural density and pressure inhomogeneities may act as perturbations. Also, if the magnetic and rotational axes are not aligned, there will always be some density enhancement near the disk foot-points of the funnel streams (Romanova et al. 2003. This would be a constant source of inhomogeneity in the disk. Our simulations have shown that even at small misalignment angles, Θ ∼ 5 • , the misalignment leads to an inhomogeneous density distribution in the disk, with two oppositely oriented density enhancements or spiral waves. The effect is even stronger at larger misalignment angles, Θ ∼ 10 − 30 • . Another source of perturbation is associated with the magnetic field lines which are trapped inside the inner regions of the disk and are azimuthally wrapped by the disk matter. This leads to increase of magnetic energy in some parts of the disk and to partial expulsion of matter from these regions, and thus to inhomogeneous distribution Figure 15. Equatorial slices of the circumstellar region at various misalignment angles. The colors represent plasma density contours, ranging from red (highest) to blue (lowest). The black line is the β = 1 line. Figure 17. Equatorial slices of the circumstellar region at various grid resolutions and times. The colors represent plasma density contours, ranging from red (highest) to blue (lowest). The black line is the β = 1 line. The grid resolutions are specified as N r × N 2 , where N r and N are the numbers of grid cells in the radial and angular directions respectively in each of the six zones of the cubed-sphere grid. of matter. This mechanism is expected to operate in real astronomical objects as well.
Concerning the role of the grid, it is unlikely, as noted above, that the discrete nature of the grid by itself leads to perturbations. But another perturbing element is the boundary between the sectors of the cubed sphere grid. Four of these boundaries cross the disk. They produce initial density and pressure perturbations at the 5% level near the diskmagnetosphere boundary, and at even larger levels at larger distances from the star, where the grid is coarser. At later times these perturbations become less important. So at early times in the simulations, this boundary effect is the most important contributor to the perturbations. That is why we often see four tongues initially. However, at later times, we often observe anywhere between 2 and 7 tongues, which shows that there is no direct influence of these boundaries on the perturbations at later times.
To check the importance of the boundaries between the grid sectors in the generation of perturbations, we introduced a perturbation in the form of a density enhancement at the inner edge of the disk at the beginning of the simulation. At a certain location in the inner disk region, we chose a sphere of diameter equal to the initial thickness of the disk, and in-creased the plasma density in the sphere by a factor of 10. Fig. 18 shows the temporal evolution of the instability in this case. We see that the perturbation generated by the blob leads to the formation of more than four tongues right from the beginning of the simulation. This shows that the initial perturbations by the boundaries of the grid and by the blob lead to very similar subsequent evolution with similar numbers of unstable tongues.
Another issue is whether the growth of the tongues is the effect of the numerical diffusivity of the code. We do not think this is the case, since the tongue growth occurs on the inner-disk dynamical timescale, which is much shorter than the diffusive timescale at this distance. Also, the numerical diffusivity decreases when the grid resolution is increased, but the instability still exists at higher resolutions.

EMPIRICAL CONDITIONS FOR THE EXISTENCE OF THE INSTABILITY
To investigate in more detail the parameter ranges over which the instability appears, we performed multiple simulation runs for a variety of values of α and P, for two values of the Figure 18. Equatorial slices of the circumstellar region at various times for our main case, with an artificially introduced density enhancement at the inner edge of the disk at t=0. The colors represent plasma density contours, ranging from red (highest) to blue (lowest). The black line is the β = 1 line. magnetic dipole moment, µ = 2 and µ = 0.5, at misalignment angle Θ = 5 • . Fig. 19 shows the resulting regimes of stable and unstable accretion. The top panel shows that, as noted before, a high accretion rate through the disk and slow rotation of the star favour the instability. The bottom panel shows the stable and unstable regimes in theṀ * − P plane, whereṀ * is the accretion rate onto the surface of the star.
Here, the boundary between the regimes has a much weaker dependence on the rotation period. This is probably because, as mentioned in §3.2, increasing the star's rotation rate takes the star closer to the propeller regime, lowering the disk density. The combination of the reduced disk density and higher α possibly produces an almost constant accretion rate onto the star's surface. It is to be noted that, as described in §3.1, the transition between the stable and unstable regimes is not sharp. Near the boundary between the stable and unstable regimes in thė M * − P plane, accretion through both tongues and funnels is seen. This limits the accuracy of the position of the boundary between the stable and unstable regions.
Decreasing the star's magnetic field increases the extent of the unstable region in Fig. 19. This indicates that smaller accretion rates through the disk are sufficient to trigger the instability. The physical effect of changing the dimensionless magnetic moment µ is to change the size of the magnetosphere. For µ = 2, the ratio of the magnetospheric radius to the stellar radius is between 4 and 5, and for µ = 0.5, it is between 2 and 3. Fig. 19 is in dimensionless units, and therefore, as noted in §2.1, can be used for a wide variety of physical situations with appropriately chosen reference values. As an illustration, we show here the critical value of the surface accretion rate separating the stable and unstable regimes for µ = 2, and the reference values of the star's rotation period, for various physical systems. One must be careful when using the reference values in Table A1, because as noted in Appendix A, a dimensionless magnetic moment of µ = 2 corresponds to a  Fig. 19, the dividing line between the stable and unstable regimes for µ = 2 is atṀ crit = 0.3 in dimensionless units, which translates into the following values for classical T Tauri stars (CTTSs), white dwarfs and neutron stars, where the star has mass M, radius R and surface magnetic field B: (i) CTTSs: This is useful because if the mass, radius and magnetic field of the star and the accretion rate are known from observations, one can use equations (1) to (5) to find out if the star is expected to be in the stable or unstable regime of accretion, a fact which is important for variability and spectral modelling. Conversely, if the nature of the accretion flow (stable or unstable) is well constrained by observations, then equations (1) to (5) can be used to constrain the star's magnetic field.
The above results have been obtained from simulations of relativistic stars (using the Paczyński-Wiita potential, as mentioned in §2), for the fiducial neutron-star parameters shown in Table A1, which correspond to r g /R * = 0.4. However, simulations for non-relativistic stars have shown that the unstable accretion has similar features, and a difference of about 20% in quantities such asṀ.

COMPARISON WITH ANALYTICAL STABILITY CRITERIA
In the simple case of a high-density fluid supported against gravity by a low-density fluid with a homogeneous magnetic field, with a plane boundary between them, the development of the Rayleigh-Taylor and Kelvin-Helmholtz instabilities in the direction perpendicular to the field is not affected by the field -all perturbation modes are unstable (Chandrasekhar 1961). This would suggest that for a star with a dipole field, azimuthal perturbations at the inner disk boundary should always be unstable. However, the inner disk usually has a relatively strong azimuthal field component B φ , which develops due to the difference between the rotation rate of the star and the Keplerian rotation speed at the inner edge of the accretion disk. An azimuthal field is expected to suppress short wavelength perturbations, and this has been observed in earlier simulations (Wang & Robertson 1984, 1985  The azimuthal field is usually found to be about (5-30)% as strong as the vertical magnetic field in our simulations. Keeping this in mind, broad conclusions about the instabilities can be drawn from the analytical results for simple cases like the one mentioned above. The Kelvin-Helmholtz instability in the azimuthal direction is suppressed if (Chandrasekhar 1961, §106), where ρ and v are the plasma density and azimuthal velocity, and the subscripts d and m refer to the disk and the magnetosphere, and we have used ρ d >> ρ m . This inequality is applicable only for a plane boundary between the fluids, but is sufficient for the rough estimates we are performing here. Rough values of the above quantities from some of our typical simulations, in dimensionless units, are as follows: We thus see that the inequality (7) is easily satisfied. The Kelvin-Helmholtz instability is completely suppressed. Azimuthal perturbations with wavelength λ are Rayleigh-Taylor unstable if (Chandrasekhar 1961, §97), where g eff ≡ g − Ω 2 r is the effective gravitational acceleration (g and g eff are positive if the acceleration is radially inwards). In the linear regime, the amplitudes of the azimuthal perturbation modes m are proportional to e imφ . Then we have λ = 2πr m /m, where r m is the magnetospheric radius. Thus, unstable modes are those that satisfy From our simulations, we have, roughly, r m ≈ 1.5, ρ d ≈ 0.7, g eff ≈ 0.1, B φ ≈ 0.3, which gives m 7, for both stable and unstable cases. Thus, this criterion successfully predicts the suppression of high-m Rayleigh-Taylor modes, but not the complete suppression of the instability in some cases. Clearly some other mechanism not accounted for here is responsible for suppressing the instability. One important factor missing from the above analysis is the effect of the radial shear of the angular velocity, which can suppress the instability by smearing out the perturbations. Spruit, Stehle & Papaloizou (1995) have performed a more general analysis of disk stability in the thin disk approximation, taking the velocity shear into account. The disk has a surface density Σ and is threaded by a magnetic field with a vertical component B z . Their criterion for the existence of the instability is In other words, Σ/B z should increase with r fast enough to overcome the stabilizing effect of the velocity shear. Fig. 20 shows the relevant quantities for this criterion, azimuthally averaged in the equatorial plane of the disk. The left column shows our main case, which has α = 0.1 and is unstable. The case in the right column differs from the main case only in that it has α = 0.02 and is stable. The radial extent shown spans the inner disk region. Note that for α = 0.1, since the accretion rate through the disk is higher than for α = 0.02, the inner disk is closer to the star. Due to this, the most significant difference between the top two panels is that the departure of the plasma orbital frequency from Keplerian in the inner disk is larger for the α = 0.1 case. The other terms in the criterion (10) do not differ significantly between these two cases. As a result, γ 2 BΣ is much larger for α = 0.1 than for α = 0.02, as seen in the bottom two panels of the figure. If one defines the transition region between the disk and the magnetosphere as that over which the plasma density changes from ρ d to ρ m , then the bottom two panels show that for α = 0.1, γ 2 BΣ > γ 2 Ω over almost the entire transition region, predicting instability, while for α = 0.02, γ 2 BΣ < γ 2 Ω over most of the transition region, predicting stability. The main driver for the instability in the α = 0.1 case is thus seen to be the stronger effective gravitational acceleration. However, it is important to note the role that the shear term γ 2 Ω plays here: although it is approximately the same in both cases, it is larger than γ 2 BΣ for α = 0.02, which suppresses the instability. We performed a similar comparison for stars with different rotation periods P. The right column of Fig. 21 shows the same case as in the right column of Fig. 20, which has P = 3 and is stable, and the case in the left column of Fig.  21 differs from that in the right column only in that it has P = 11 and is unstable. This time, since α is the same in both cases, the accretion rate through the disk is also the same, and the inner disk region is similar in both cases. The angular velocity Ω is almost constant in most of the transition region in both cases. But in the P = 11 case, in the inner part of the transition region, Ω drops off sharply towards the star as magnetic braking tries to force the plasma to corotate with the star. This increases the shear term γ 2 Ω in this region. But precisely because the plasma is slowed down, the effective gravity in the inner transition region is higher, with the result that γ 2 BΣ > γ 2 Ω over almost the entire transition region, once again predicting instability. For the P = 3 case, once again, the shear is almost the same as in the outer transition region for P = 11, but the weaker effective gravity reduces γ 2 BΣ to below γ 2 Ω over almost the entire transition region, predicting stability. Thus, criterion (10) works very well for these cases.
Li & Narayan (2004) have performed linear perturbation analysis on a cylindrically symmetric initial configuration. This analysis is expected to be applicable to the midplane of the disk. They derive analytical stability criteria for this model in a few special cases. One of those cases is where the disk and magnetosphere have constant but different densities, and the angular velocity is constant in the magnetosphere and has a jump at the magnetospheric radius r m . This is the case that is most closely applicable to accretion flows around magnetized stars. The angular velocity profile in the disk is assumed to be Ω = Ω d (r/r m ) −q , with q = 0 or 2. Their criterion for instability is and the vorticity, The term on the first line in criterion (11) is the Rayleigh-Taylor term, and the one on the second line is the Kelvin-Helmholtz term. Since ρ d >> ρ m , µ ≈ 1, and (11) becomes The Kelvin-Helmholtz term is again seen to drop out. For unstable modes, we then have, We see that low vorticity and a high effective gravitational acceleration (measured by Ω eff ) are conducive to the development of the instability. Fig. 22 shows the relevant quantities for this criterion in the inner-disk region for the same two simulation runs as in Fig. 20. The quantities are again in the equatorial plane, and azimuthally averaged. As noted earlier, the effective gravity is stronger for α = 0.1. Also, the vorticity is seen to be slightly smaller. Although this criterion does not work as well as criterion (10), m crit has a definite tendency to be smaller for α = 0.1.

RESULTS AND DISCUSSION
Accretion through instabilities at the disk-magnetosphere interface is found to occur for a wide range of physical parameters. It results in tall, thin tongues of gas penetrating the magnetosphere and travelling in the equatorial plane. The tongues are very transient, and grow and rotate around the star on the inner-disk dynamical timescale. The number of tongues at a given time is of the order of a few, irrespective of the grid resolution. Near the star, the tongues are threaded by the magnetic field lines, and form miniature funnel-like flows, which gives them a characteristic wishbone shape. They deposit matter much closer to the star's equator than true funnel flows do. Each tongue produces its own hot spots on the star's surface, and as a result, the hot spots also change on the inner-disk dynamical timescale. Our simulations were performed for stars with magnetospheric radius 2-5 times the radius of the star. The instability is associated with high accretion rates, and coexists with funnel flows for a relatively broad range of accretion rates. The instability is expected to occur in most accreting systems for typical values of mass, radius, surface magnetic fields and accretion rates.
The instability is suppressed if the misalignment angle Θ between the star's rotation and magnetic axes is large (Θ 30 • ). For Θ 30 • , when the accretion is through instability, the rotational latitude at which most of the accreting matter falls on the star is independent of the misalignment angle.
The number of tongues we see is of the order of a few. The natural question that arises is whether this result is robust with respect to the grid resolution. From theory, we expect small-wavelength Rayleigh-Taylor modes to grow faster, at least in the linear regime (e.g. Chandrasekhar 1961). Solar prominences, which are also subject to the Rayleigh-Taylor instability (Anzer 1969), have a high degree of filamentary structure. If the inner edge of the disk has a similar structure, then increasing the grid resolution would resolve ever thinner tongues, increasing their number. However, we think that this is not the case. Wang & Robertson (1985) investigated the late stages of the Rayleigh-Taylor instability with a plane boundary between the two fluids, in two dimensions. One of their important observations was that although short wavelength modes initially grow faster, their growth is limited by viscous damping. Further evolution is possible only towards larger spatial scales, which is accomplished by the clumping together of the smaller wavelength structures. We see the same behaviour in our simulations. Two other important stabilizing factors for the Rayleigh-Taylor instability are the tension of the azimuthal magnetic field and the radial shear of the angular velocity. The azimuthal magnetic field develops because the stellar magnetic field threads the disk, and is twisted by the inner-disk plasma, which usually rotates faster than the star. The magnetic field tension is expected to suppress short-wavelength perturbations (Chandrasekhar 1961), and this has been observed in earlier simulations (Wang & Robertson 1984, 1985Rastätter & Schindler 1999b;Stone & Gardiner 2007a,b). Our results confirm these observations. The shear of the angular velocity also has a stabilizing effect, since it tends to smear out the perturbations. These, we think, are the reasons why we do not see a large number of tongues. As a final note, the number of grid cells spanning the width of a fully grown tongue in our highest resolution runs is of the order of 10, which indicates that the grid is fine enough to resolve much thinner tongues than those we see.
Although the shear tends to suppress the Rayleigh-Taylor instability, if the shear is large enough, one would expect the Kelvin-Helmholtz instability to develop. However, the angular velocity profile in the inner disk is generally close to being flat, since the stellar magnetic field tries to force the innerdisk plasma closer to corotation. Rough comparisons with the Kelvin-Helmholtz criterion for a plane boundary between two fluids (Chandrasekhar 1961) suggest that the shear that we observe is not strong enough to overcome the stabilizing effect of the azimuthal magnetic field, although it is strong enough to suppress the Rayleigh-Taylor instability in some cases.
Analytical stability criteria have been derived by a number of authors. We compared our simulations with two of them here. Spruit & Taam (1995) derived a stability criterion for thin disks, which is found to work reasonably well. Li & Narayan (2004) derived a stability criterion for a cylindrically symmetric initial equilibrium, which is expected to be applicable in the disk midplane. This criterion works only approximately. This is possibly because the particular angular velocity profiles considered by those authors, and the assumption of a sharp magnetospheric boundary, which are necessary for analytical tractability, do not approximate the conditions in realistic disks as closely as one would wish.
One of the most interesting observational consequences of accretion through instabilities is the effect on the variability. Light curves associated with accretion through funnel streams show clear periodicity, but those associated with accretion through tongues often do not. This of relevance to young stars, whose light curves often lack clear periodicity (e.g. Herbst et al. 2000Herbst et al. , 2002. Such light curves would not preclude the existence of an ordered stellar magnetic field if accretion through instabilities were taken into account. On the other hand, we sometimes see that a certain number of tongues dominates, which may lead to quasi-periodic oscillations in the lightcurves (Li & Narayan 2004). These QPOs may be important for understanding the Type II (accretion-driven) bursts in LMXBs. Comptonization of photons by high-energy electrons may lead to only a small departure of the lightcurve from the thermal ones obtained using the approximation of isotropic black-body radiation (Poutanen & Gierliński 2003) so that the QPO features may survive comptonization (see, however, Titarchuk et al. 2007). From a different angle, the instability found in our simulations may be related to the unstable "corotation" mode found by Lovelace & Romanova (2007) near the disk/magnetosphere boundary by a linear WKB stability analysis.
Due to the stochastic nature of the tongues, it is possible for periods of such quasi-periodic oscillations to alternate with periods of no detectable pulsations. Also, for accretion rates near the critical accretion rate associated with the instability, accretion through both tongues and funnels is observed. This raises the possibility of changes in the accretion rate leading to switching between funnel-and tongue-dominated accretion. These phenomena might be related to the intermittent pulsations observed in some LMXBs (see, e.g., Kaaret et al. 2006;Altamirano et al. 2007;Casella et al. 2007;Galloway et al. 2007;Gavriil et al. 2007). We plan to investigate variability associated with unstable accretion in future work.
The other aspect of variability is spectral variability of CTTSs. The spectral line profiles are dependent upon the density and velocity of the accretion flow (see, e.g., Symington et al. 2005;Kurosawa et al. 2006), and will therefore be affected by instabilities. Recently spectral lines were calculated for CTTSs based on our 3D MHD model of stable accretion through funnel streams using a 3D radiative transfer code (Kurosawa et al. 2007). We plan to undertake similar studies for accretion through instabilities.
In the case of young stars surrounded with gas/dust disks where planets are forming and migrating inward, the tongues may support inward migration, so that the magnetospheric gap, which halts migration (see, e.g., Lin, Bodenheimer & Richardson 1996;Romanova & Lovelace 2006;Papaloizou 2007), may form only in the state of stable accretion.