Magnetic shear-flow instability in thin accretion discs

Abstract

The possibility that the magnetic shear-flow instability (also known as the ‘Balbus-Hawley’ instability) might give rise to turbulence in a thin accretion disc is investigated through numerical simulations. The study is linear and the fluid disc is supposed to be incompressible and differentially rotating with a simple velocity profile with Ω∝R^−q. The simplicity of the model is counterbalanced by the fact that the study is fully global in all three spatial directions with boundaries on each side; finite diffusivities are also allowed. The investigation is also carried out for several values of the azimuthal wavenumber of the perturbations in order to analyse whether non-axisymmetric modes might be preferred, which may produce, in a non-linear extension of the study, a self-sustained magnetic field.

We find the final pattern steady, with similar kinetic and magnetic energies and the angular momentum always transported outwards. Despite the differential rotation, there are only small differences for the eigenvalues for various non-axisymmetric eigensolutions. Axisymmetric instabilities are by no means preferred; in fact for Prandtl numbers between 0.1 and 1, the azimuthal wavenumbers m=0,1,2(10¹⁶ g s⁻¹). All three quantities appear to be equally readily excited. The equatorial symmetry is quadrupolar for the magnetic field and dipolar for the flow field system. The maximal magnetic field strength required to cause the instability is almost independent of the magnetic Prandtl number. With typical white dwarf values, a magnetic amplitude of 10⁵ G is estimated.

accretion, accretion discs, instabilities, MHD

1 Introduction

The long-standing problem of the generation of turbulence in various astrophysical situations, in which ordinary microscopic viscosities seem to be unable to produce the observed instabilities, has perhaps found a solution in recent years in the so-called ‘Balbus-Hawley instability’, in which the presence of a magnetic field has a destabilizing effect on a differentially rotating flow, provided that the angular velocity decreases outwards with the radius (Velikhov 1959; Chandrasekhar 1961). After the first work by Balbus & Hawley (1991), several other studies have been carried out to study the non-linear evolution of the instability (Hawley & Balbus 1991; Brandenburg et al. 1995; Hawley, Gammie & Balbus 1995; Matsumoto & Tajima 1995), but because of the complexity of non-linearities, all of these studies were aimed at understanding the ‘local’ properties of the instability. Recent investigations, however, have pointed out the importance of boundaries on the instability growth rates (Curry, Pudritz & Sutherland 1994; Curry & Pudritz 1995, 1996; Kitchatinov & Rüdiger 1997; Kitchatinov & Mazur 1997). All these studies were linear, a fully global non-linear approach being in its infancy (Drecker et al. 1997).

We investigate here the effect of the instability in a simple numerical model for a slim accretion disc in which a magnetized fluid is contained between two finite (in the axial direction) cylinders and is differentially rotating with a velocity profile Ω∝R^−q, R being the radius in cylindric coordinates (R, φ, Z). We deal with three cases: q=1, q=1.5 and q=2. All these flows are stable (only marginally in the case q=2) with respect to the standard hydrodynamic Rayleigh criterion for Taylor-Couette flow, in such a way that an eventual instability is caused by the presence of the magnetic field. The choice q=2 corresponds to a constant specific angular momentum along the radius between the two cylindrical surfaces. Such an angular velocity profile seems to fit the differential rotation for thick discs well (Papaloizou & Pringle 1984) while a Keplerian law (q=1.5) should be more suitable for thin discs. It is generally accepted that the index q in the rotation law depends in turn on the position and that in some zones of the disc, e.g. in the inner part, a profile with q=2 is not unrealistic at all (Abramowicz, Brandenburg & Lasota 1996). In addition, the fluid is supposed to be incompressible, to get rid of the numerical problems related to the presence of sound waves. The system is embedded in an external vertical, uniform magnetic field that should drive the instability.

Flow and field perturbations are developed after the azimuthal Fourier modes exp(imφ). Although the flow and field quantities are then complex numbers, only their real parts have an actual physical meaning. Our main result is that the non-axisymmetric modes also seem to prevail at high values of the Hartmann number, where in a previous work (Kitchatinov & Mazur 1997) the mode m=0 was always found to be the most preferred one. This opens the possibility of believing that, in a fully non-linear approach, the instability could generate a self-excitation of the magnetic field, as was already found in local studies (Brandenburg et al. 1995). Furthermore, a dominance of the mode m=2 at high Hartmann numbers could be relevant for models of the galactic dynamo.

2 Equations and Model

We solve the incompressible, dissipative, linearized MHD equations completely numerically in cylindrical coordinates for a global model. The basic rotation law is  

(1)

with R_Ω as the inner radius of the disc, and the rotating fluid may be threaded by a strictly vertical magnetic field,  

(2)

consistent (as a very special case) with Ferraro's law (cf. Stone & Norman 1994). Extra radially dependent azimuthal fields are also possible in this approximation (Curry & Pudritz 1995; Ogilvie & Pringle 1996; Terquem & Papaloizou 1996; Papaloizou & Terquem 1997). By linearizing around this equilibrium state and using dimensionless quantities, our equations read  

(3)

 

(4)

 

(5)

where p denotes the pressure, u the velocity perturbation and B the magnetic perturbation.

In the equations above, A and L are the linearized advection term and Lorentz force given by  

(6)

and  

(7)

while  

(8)

and  

(9)

are the linearized electromotive forces induced by the basic differential rotation and the flow perturbations. In the equations above, lower-casers are the dimensionless radius rRH, and aspect ratio r_ΩR_ΩH, respectively. Finally,  

(10)

is the Hartmann number (η, ρ and ν being the magnetic diffusivity, density and viscosity),  

(11)

is the shear-flow parameter (‘dynamo number’) and Pm is the magnetic Prandtl number. The quantity H is the half-thickness of the disc, which is also the unit length for the non-dimensionalization. The unit time is the magnetic diffusion time τH²η, and velocities are normalized with HΩ₀. The permeability of a vacuum is written as μ₀ but substituted by 4π for calculations in cgs. It is a time-stepping code using a global energy quenching procedure of the form to find the lowest excitation regime, where E is the total energy of the system (see Elstner, Meinel & Rüdiger 1990). A staggered mesh of 81 by 81 grid points was used to solve the system.

In the equations (4)–(5), Ha, C_Ω and Pm are free parameters. We solved the equations in order to obtain a representation of the curves of marginal stability in the Ha-C_Ω plane. The other free parameter in our model represents the radius of the inner rotating cylinder; in particular r_Ω fixes the aspect ratio of the accretion disc. Standard values of this parameter range between r_Ω∼30 for T Tauri stars and r_Ω∼100 for neutron stars (Lioure & Le Contel 1994). For computational reasons we set r_Ω=30 and only considered a disc reaching from r=30 to r=40.

We chose to impose rigid as well as stress-free boundary conditions for the velocity field, both applying atall boundary layers, and set pseudo-vacuum boundary conditions on the magnetic field, i.e. the components of the magnetic flux B tangential to the surface are vanishing. This choice (‘pseudo-vacuum’) is also made for simplicity in imposing the conditions for the flow, since by using a staggered mesh the grid-points for the velocity are not defined directly at the boundaries. It also ensures the necessary vanishing of the Poynting flux vector in the normal direction (cf. Ruzmaikin, Shukurov & Sokoloff 1988). It is generally accepted that the Balbus-Hawley instability is influenced by the presence of the boundaries (Curry & Pudritz 1995). It is thus of high importance to find out how the various boundary conditions influence (say) the symmetry of the solutions with respect to the rotation axis or with respect to the equator. Local simulations for closed boxes cannot answer such questions, with their high relevance to observations.

3 Results

3.1 The neutral-stability lines

The neutral-stability lines for m<;3 in the C_Ω-Ha plane for q=1 and q=2 and various Pm are given in Fig. 1 for rigid boundary conditions. As there is little qualitative change with q, we only show results for q=1.5 for the stress-free case in Fig. 2. Above the curves the rotation law is unstable while below the eigenvalues the stratification is stable. Instabilities in a disc with stress-free boundary conditions are excited more easily than in a rigid-boundary disc. A direct comparison with the plots of Kitchatinov & Mazur (1997) for m=0 shows that they are similar to ours, in spite of the difference in the definitions. We checked that the curves have no intersection with the C_Ω-axis, because we found the energy decreasing exponentially for Ha=0. This is because, in the absence of a magnetic field, the basic flow is linearly stable by the Rayleigh criterion.

Figure 1.

Neutral-stability lines for various azimuthal wavenumbers (m=0,1,2) and rigid boundary conditions. The parameters are q=1 (left) and q=2 (right). r_Ω=30, Pm=0.1 and Pm=1. The isolated symbols for q=2 denote simulations for the small Prandtl number Pm=0.01.

Figure 2.

Neutral-stability lines for various azimuthal wavenumbers (m=0, 1, 2) and stress-free boundary conditions. A Keplerian rotation law is used (q=1.5). Again, r_Ω=30, Pm=0.1 and Pm=1.

We find for higher values of the Hartmann number (for stronger magnetic fields) either the mode m=1 or m=2 with slightly lower excitation eigenvalues. As a main result of the present study, very similar excitation conditions for both axisymmetric and non-axisymmetric modes are found. It is worth noting that in the case of the sphere (another case in which the fluid is bounded in all directions) Kitchatinov & Rüdiger (1997) found that the mode m=1 was more easily excited for a strong magnetic field. Also the mode m=2 was found to be preferred to the axisymmetric one, though only for a rather strong magnetic field.

Fig. 1 yields linear relations for the Hartmann number as a function of the ‘dynamo number’ such as  

(12)

for the maximal possible Hartmann number. ε is of order unity (ε≲0.8) for Pm=1. It is possible to see from the figure that the scaling of the quantity ε with Pm is approximately  

(13)

If the curves in Fig. 1 are taken for as a function of Ha, the neutral-stability lines are very close together for all values of both Pm and m, even for the (few) examples obtained for Pm=0.01. Note that in this formulation the critical number for the onset of instability proves to be  

(14)

In terms of the Alfvén velocity V_A and the speed of soundc_ac is  

(15)

similar to the value of Papaloizou & Szuszkiewicz (1992) for stabilization of the magnetic shear flow.

For the largest possible values of the magnetic field, one finds from (12)  

(16)

The sound velocity and also the radial density profile are known for the standard accretion disc model (cf. Frank, King & Raine 1985), hence

 

17

M₁M/(10⁹ cm), M₁M/M_? /(10¹⁶ g s⁻¹). /(10¹⁶ g s⁻¹). All three quantities are of the order of unity for discs around white dwarfs. The definition of the viscosity parameter a_SS is given below in (19). With its standard value of 0.01, the maximal possible magnetic field for the Balbus-Hawley instability is  

(18)

in accordance with the observed surface field strength for white dwarfs of 10^5…9 G. After (16), the ratio of the gas pressure to the maximally possible magnetic pressure is of order unity.

3.2 The angular momentum transport

Although our study is linear, it is worth analysing the efficacy of the instability for the angular momentum transport. The effective viscosity that should bring the angular momentum transport is generally parametrized through the quantity α_SS(Shakura & Sunyaev 1973), defined in terms of the Reynolds and Maxwell stress tensors by  

with the normalization  

(19)

where α_SS is called the viscosity alpha. Hence  

(20)

i.e. in dimensionless units  

(21)

(Brandenburg et al. 1996). In our linear approach, the amplitudes of all the quantities are unknown as they are free of a common positive or negative factor. Hence, the sign of quadratic terms can be computed, as well as ratios of fluctuating quantities.

In Fig. 3, the calculated values of α_ss are shown for a point on the marginally stable curves for the modes m=0 and m=2, respectively. It is normalized to unity. The averages are taken in the azimuthal and vertical directions, so that α_ss is a function of the radius. In both cases α_ss is positive, meaning that the instability is effective in transporting angular momentum outwards. Owing to the fact that the eigenfunctions are concentrated close to the inner boundary of the integration domain, the correlation between the R and φ components of u and B is non-vanishing only in proximity to r_Ω.

Figure 3.

The (normalized) viscosity alpha for q=1.5, Ha=32, Pm=0.1, m=0 (solid line) and m=2 (dashed line), which is positive, obtained with stress-free boundary conditions.

It is also possible to give an estimate of the relative strength of the Reynolds stress in comparison with the Maxwell stress by measuring the quantity  

(22)

plotted in Fig. 4. The sign is positive in both cases, meaning that both the Reynolds and Maxwell stresses are working to transport the angular momentum in the same direction. The maximum γ is of order 1.0 to 2.0 so that the contributions of the velocity field and magnetic stress to this transport are similar.

Figure 4.

The ratio between kinetic and magnetic angular momentum transport for q=1.5, Ha=32, Pm=0.1, m=0 (solid line) and m=2 (dashed line), obtained with stress-free boundary conditions.

The same can be done with the total energy. The ratio of kinetic and magnetic energy (in our units) is  

(23)

taken over the complete cylinder. The result for the same model as that in Figs 3 and 4 is plotted in Fig. 5. The magnetic energy slightly dominates the kinetic energy.

Figure 5.

The ratio between kinetic and magnetic energy for q=1.5, Ha=32, Pm=0.1, m=0 (solid line) and m=2 (dashed line), obtained with stress-free boundary conditions.

3.3 The eigenfunctions

In Fig. 6, the induced toroidal magnetic field is shown. It is, obviously, an equatorially symmetric (‘S’) field. Correspondingly, the associated zonal kinetic flow has an antisymmetric structure with respect to the equator (Fig. 7). The induced fields and flows here are concentrated at the inner boundary.

Figure 6.

The (equatorially symmetric) eigenfunction for the toroidal (left) and poloidal (right) magnetic field. We used Ha=32, q=2, Pm=0.1, m=0 and rigid boundary conditions for this example.

Figure 7.

The (equatorially antisymmetric) eigenfunction of the kinetic zonal flow (left) and the meridional drift (right). We used Ha=32, q=2, Pm=0.1, m=0 and rigid boundary conditions for this example.

The zonal flow with stress-free boundary conditions is very similar to the flow with rigid boundary conditions shown in Fig. 7. Those contour lines, which are strongly bent near the upper and lower boundaries of the disc in Fig. 7, are just open in the z-direction in the stress-free case. The similarity of flow patterns helps us to understand the similarity of the respective neutral-stability lines in Figs 1 and 2.

The most important question relevant to such eigenfunctions is the fate of electrically charged dust grains in the flow field system. Do they concentrate in the equatorial plane or not under the influence of the various electromagnetic forces? The ionization states of protoplanetary discs have been studied by Dolginov & Stepinski (1994) and Stepinski & Valageas (1996).

4 A Disc Dynamo Model?

Our numerical simulations of thin accretion discs revealed the existence of stationary laminar flow and field patterns for a certain amplitude range of external magnetic fields. There are even indications that for stronger magnetic fields (increasing Hartmann number) non-axisymmetric modes are more easily excited. Owing to the Cowling theorem, only these modes could excite, in a non-linear regime, a self-maintained axisymmetric magnetic field, so that the results seem to be promising in view of a possible non-linear extension of the code for investigating the possibility of a dynamo effect. After our experiences with large-scale dynamo models, however, the existence of a vertical density stratification should be necessary for dynamo self-excitation (cf. Brandenburg et al. 1995; Hawley, Gammie & Balbus 1996; Stone et al. 1996). Before moving to a non-linear study, however, there are several possible improvements that could be made. The velocity profile could be generalized and different choices for the structure of the external background magnetic field (e.g. a dipole) may be made. A forthcoming extension of the study concerns the choice of the parameters in a model that fits a more ‘realistic’ situation. Runs are in progress with different values of the r_Ω parameters to simulate thinner accretion disc structures and values of the Prandtl number of Pm≲10⁻², closer to typically estimated values for accretion discs.

Let us finally stress that the dominance of the mode m=2 for the higher Hartmann numbers shown in Figs 1 and 2 could be relevant in the case of a galactic dynamo, although the intrinsic limitations of our model prevent us from any direct conclusion.

Acknowledgments

Our thanks is expressed to L. L. Kitchatinov (Irkutsk) and M. Schultz (Potsdam) for support during the work for the present study. We also acknowledge the interesting remarks of the anonymous referee.