Accretion disc viscosity: how big is alpha?

0.16 MB·English

by  A. R. King

Mon.Not.R.Astron.Soc.000,1–8(2006) Printed5February2008 (MNLATEXstylefilev2.2)
Accretion disc viscosity: how big is alpha?
A. R. King1, J.E. Pringle1,2,3 and M. Livio3
1Theoretical Astrophysics Group, Universityof Leicester,LeicesterLE1 7RH
2Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA
3Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
7
0
0
5February2008
2

n
a ABSTRACT
J

9 Weconsiderobservationalandtheoreticalestimatesoftheaccretiondiscviscosity
2 parameterα. We find that in thin, fully-ionizeddiscs, the best observationalevidence
suggests a typical range α ∼ 0.1−0.4, whereas the relevant numerical simulations
1
tend to derive estimates for α which are an order of magnitude smaller. We discuss
v
possible reasons for this apparent discrepancy.
3
0
Key words: accretion, accretion discs
8
1
0
7
0 1 INTRODUCTION is aquantitywhose properties are tobe determined experi-
h/ mentally.1
Accretion discs are believed to be present in a wide variety
p Thisalpha–prescriptionallowsformalclosureofthesys-
of astronomical systems, and have been a major research
- temofequationsdescribingathindisc,eventhoughthereis
o topic for several decades (see e.g. Pringle, 1981 and Frank,
no presumption that α is anything other than an unknown
r King & Raine, 2002). For much of this time theorists have
t dimensionless scaling variable, although it is often assumed
s hadproblemsunderstandingthefundamentaldrivingmech-
to be a constant. Gratifyingly, many of the properties of
a
anism transporting angular momentum outwards and thus
: steadythindiscsturnouttohaveratherweak dependences
v allowing mattertospiralinwardsinadisc.Thismechanism
on α. Ignorance of the physical properties or strength of
i is usually called ‘viscosity’, and appears in virtually all of
X the angular momentum transport process, represented in
disc theory. Ideas about discs nevertheless gained credibil-
dimensionless fashion by α, is thus much less of an obsta-
r ity because of two factors. The first was that the radial
a cle to practical application of this simple picture than one
distribution of effective temperature across a steady disc
might suppose. The perceived success of these applications
(T(R) ∝ R−3/4) is independent of the viscosity, being just
amounts to noting that fairly similar values of α appear to
astatementofenergyconservation,andisinreasonableac-
give reasonable agreement with observations of many sys-
cord with both continuum spectra and eclipse mapping of
tems.
cataclysmicvariables(CVs).Theseareclosebinarysystems
A physically plausible theory of the underlying causes
where a white dwarf accretes from a low–mass companion
of disc‘viscosity’ has emerged overthelast fifteen years. In
via a disc. The second was the devising of a physically–
theirseminal paperShakura& Sunyaev(1973) argued that
motivated, dimensionless scaling of the kinematic viscosity
magnetic fields are the likely way in which a shearing disc
ν as
flow transports angular momentum from rapidly–rotating
fluidtomoreslowly–rotatingfluidfurtherout.Thisconcept
was given impetus with the realization (Balbus & Hawley,
ν =αc H (1)
s 1991) that what is now called the magnetorotational insta-
bility(MRI)canprovidethenecessaryfeedbacktomaintain
a magnetic dynamo in accretion discs. The MRI forms the
(Shakura&Sunyaev,1973).Herec isthelocalmeansound
s
basis of all current theoretical simulations of this process.
speed in the disc, and H ∼(c /v )R is thescaleheight per-
s φ
These have not yet reached the point where direct com-
pendicular to the disc plane at radius R, where v is the
φ
parison with observation is possible, in terms for example
azimuthal velocity. In a thin disc (i.e. one with cooling effi-
of being able to predict the spectrum of radiation emitted
cientenoughthatH ≪R)wherethevelocityv isveryclose
φ
to the Kepler value v = (GM/R)1/2 (where M is the ac-
K
cretingcentralmass),thesequantitiesarewell–defined,soα 1 It is important to realise that α is a mean quantity, averaged
isadimensionlessquantityspecifyingthelocalrateatwhich perpendiculartothediscplane.Itiswell-definedonlyifthedisc
angular momentum (strictly speaking, that component or- isthin, i.e.H ≪R. In athick disc, withH ∼R, α has no clear
thogonaltothediscplane)istransported.Theparameterα meaning
2 A. R. King, J. E. Pringle & M. Livio
by accretion discs. Instead the main reason for optimism weighted towards larger radii where the viscous timescales
has been the belief that these simulations do demonstrate are longest.
thefeasibilityofaself-maintainingprocesswhichtransports
angular momentum in the required manner.
2.1 Dwarf nova outbursts
The point we wish to address in this paper is that it
is also often assumed that numerical simulations produce The largest body of evidence here comes from the light
formal valuesof αwhich agree with thoseinferred from ob- curves of dwarf novae, which are a subclass of cataclysmic
servation. The purpose of this paper is to examine how far variables (CVs) which undergo outburstsat irregular inter-
thisistrue.InSection2wediscussthoseastronomical phe- vals (Warner, 2003). There is now general agreement that
nomena which give thestrongest observational evidence for theseoutburstsresult from thepresence of ionization zones
the value of α. These are time-dependent discs involved in withinthedisc,allowing thistoswitch betweenacool, low-
the outbursts of dwarf novae and of the X-ray transients. ionization, low–viscosity state, and a hot, highly ionized,
The discs in these systems are fully ionized, because of the high–viscosity state (see Lasota, 2001 for a recent review).
natureoftheoutburstmechanism,andsocorrespondclosely In thehot state thedisc evolves on the viscous timescale
tothethemajority ofthenumericalsimulationswhichcon-
R2
sider full MHD. For cooler discs which are not sufficiently tvisc∼ (2)
ν
ionizedsothatthemagneticfieldisnottiedstronglyenough
tothediscgas,theMRIislikelytobelessvigorousoreven for a time, before a cooling front propagates through the
non-operative, so reducing the expected value of α (Gam- disc and returns it to the cool state. The initial slow decay
mie,1996).InSection3wediscusstheestimatesofthevalue of the outburst thus allows estimates of α in this hot state.
ofαderivedfromnumericalsimulations.Weconcentrateon Estimates of α can therefore be obtained from theoretical
thosesimulationswhichdonotimposeanexternalseedfield modelling of the outburst lightcurves. Since the disc sizes
threadingthewholedisctodrivetheMRI,asneitherinthe are known from the system properties, and since the disc
dwarf novae nor in the X-ray transients is there a plausi- temperatures are known from thespectra thusdetermining
ble source for such a global field. We note that the most H/R,observationoftheevolutiontimescaleoftheoutbursts
recent computations obtain values of α smaller than those gives a reasonably well-determined estimate of the viscous
required byobservations byat least an orderof magnitude, timescale and hence of α.
andoftenmore.InSection4wediscussthelimitationsunder Smak (1999) considers the observed relation (Bailey,
whichnumericalsimulationsofaccretiondiscshavetooper- 1975) between the decay rate and the orbital period. Since
ate, drivenboth bythespeed ofcurrent computersand the the latter largely fixes the orbital separation and thus the
nature of the numerical algorithms, especially the bound- discsizeR,thedecayratemeasures αdirectly through(2).
ary conditions. We point out that most of the limitations Other estimates (Smak, 1998, 1999) use the delay between
are likely toact in thedirection of reducingα. In Section 5 the peak of the outburst in the optical and the UV, which
wediscussthepossibilitythatthefieldsgeneratedbyaccre- resultsfromthediscclosingasmallcentralholearoundthe
tiondiscsaremoreglobalthancanbeeasilyaccommodated white dwarf on a viscous timescale tvisc. These estimates
incurrentnumericalsimulations,togetherwiththepossible come from the collective properties of a large sample of
consequences of, and complications arising from, the pres- dwarf novae. Other estimates use detailed observations of
ence such global fields. theoutburstlightcurvesofindividualsystems(Schreiberet
al.,2003, 2004 [SSCygandVWHyi];Cannizzo2001a [VW
Hyi,UGemandSSCyg];Cannizzo2001b [WZSge];Buat–
M´enardetal.,2001[ZCam]).Allofthesepapersagreethat
2 OBSERVED CONSTRAINTS ON ALPHA α must lie in a fairly narrow range α≃0.1−0.3.
Since α is a dimensionless measure of the viscosity2 its
propertiesneedtobedeterminedfromobservations.Webe- 2.2 Outbursts of X-ray transients
gin therefore by considering observationally determined es-
timates of α. Since as we remarked above, steady thin disc A second class of accreting binaries that have outbursts is
theory has only a fairly weak dependence on α, by far the thatofthesoftX–raytransients(SXTs)inwhichtheaccre-
most reliable and direct way of estimating this quantity is torisablackholeorneutronstarratherthanawhitedwarf
toconsidertime–dependentdiscbehaviour.Thesizeof αis as in CVs (Lewin & van der Klis, 2006). Even for similar
directly proportional to the rate at which angular momen- orbital parameters, SXT outbursts are considerably longer
tumistransportedwithinthedisc,andsoisdirectlyrelated than those of dwarf novae (months rather than days), and
to the timescale on which a disc can evolve. We note that have a different lightcurve shape (exponential for short or-
even so it is not possible to determine detailed properties bital periods). This at first presented a challenge to theory.
of α, for example how α might depend on radius, or H/R. There were initial attempts to explain this by devising dif-
Thus the observed values of α correspond to some appro- ferentad hocformsofviscosity(forexamplesettingαtobe
priate mean valueoverthediscbeing modeled, presumably a function of H/R; Cannizzo, Chen & Livio, 1995). How-
ever observations show that the discs in SXTs (and indeed
in all low–mass X–ray binaries) are optically muchbrighter
2 Notethatαisnotameasureofanisotropicviscosityasappears, than expected on the basis of the accretion rate revealed
forexample,intheNavier-Stokesequation,butis,rather,strictly bytheirX–rayluminosities. Thisextralightcanbedirectly
onlyameasureoftheverticallyaveragedratiobetweenthe(R,φ)- attributed to irradiation of the outer parts of the disc by
components ofthestressandtherateofstraintensors. some of the central X–rays (van Paradijs & McClintock,
Accretion disc viscosity: how big is alpha? 3
1994). King & Ritter (1998) pointed out that this would aplanet(Lodato&Clarke,2004).Thus,eitherthereissome
force most of the disc to remain in the hot, high–viscosity subtletyatworkhere,orthethermaldiscoutburstmodelfor
stateuntilasignificantfractionofthediscmasshadaccreted FUOrioutburstsdoesnotwork.Inthisregardwenotethat
on to the central black hole/neutron star. This explained otherpossibilities forcausingFUOrioutburstshaveindeed
thelonger duration of SXToutburstscompared with dwarf been discussed, such as a collision of a protostar and disc
novae, and indeed their exponential shape at short orbital with another star (cf. Bonnell & Bastien, 1992; Reipurth &
periods(wherethewholeofthedisccanbeefficientlyirradi- Aspin,2004).
ated). Under the assumption of efficient irradiation Dubus
et al. (2001) made detailed models of complete SXT light
curvesand found α≃0.2−0.4. 2.4 Summary
Weconcludethatinthemostclearcutcasesthereappearsto
bestrongobservationalevidencethatvaluesofα=0.1−0.4
2.3 Other systems which yield estimates of α
are required to providea good description of thebehaviour
2.3.1 Variability in AGN of fully-ionized, thin accretion discs.
A somewhat more indirect method of estimating α is sug-
gested by Starling et al. (2004) who look at the optical
variability of active galactic nuclei (AGN) on timescales of 3 THEORETICAL ESTIMATES FROM
monthstoyears.Starlingetal.(2004)measureatwo–folding NUMERICAL SIMULATIONS
timescale,definedasthetimescaleoverwhichtheopticallu-
minosity changes by a factor of two. They assume that the SincetheworkofBalbus&Hawley(1991) therehasbeena
optical emission is generated in a standard, fully ionized, large numberof publications investigating theproperties of
thin disc, and that the two–folding timescale is given by theMRI and its relevance toviscosity in accretion discs.
disc’s local thermal timescale (Pringle, 1981) Theoretical simulations of disc viscosity come in two
flavours – those which assume a superimposed seed net
1 R poloidal field, and those which do not. Hawley, Gammie &
tth ∼ . (3)
α(cid:18)vK(cid:19) Balbus(1995)showedthatsimulationswithasuperimposed
net B tend to yield estimates of α larger by an order of
They find that 0.01 ≤ α ≤ 0.03 for 0.1 ≤ L/L ≤ 1. Star- z
E magnitude than those which do not. They also found that
ling et al. notethat thesevaluesof α arereally lower limits
the value of α produced depended almost linearly on the
because data sampling means that they might miss shorter
magnitudeoftheexternallyimposedB .Thosesimulations
timescales. z
which do not have an externally imposed B mostly start
z
with either a small seed toroidal field, or alternate regions
of positive and negative vertical field B .
2.3.2 Protostellar accretion discs z
It seems to us an unlikely proposition that each disc
Youngpre-main sequencestars are often surrounded byac- is subject to a superimposed, immovable, net poloidal field
cretion discs (e.g. Hartmann, 1998). Estimates of the life- component,ofexactlytherightmagnitudetogiverisetothe
times of these discs can be obtained by comparing disc fre- α in the right range. For a typical dwarf nova in outburst
quencies among stars of different ages. Estimates for α in we find the equipartion field in the outer disc regions to be
protostellar(TTauri)discs,basedonevolutionarylifetimes, (Shakura& Sunyaev,1973; Frank,King & Raine, 2002)
aregivenbyHartmannetal.(1998).Theygiveestimatesof
pαle≈sd0i.s0c1usasteddiisncSreacdtiiionRs∼2.110an−d120.20iAnvUo.lvAelalcocfretthieonexdaismcs- Beq ≃5.2×104(cid:16)0α.1(cid:17)−9/10(cid:18)101M8˙g/s(cid:19)17/40(cid:18)MM⊙(cid:19)7/16
which are sufficiently hot that they are fully ionized. How- R −21/16
ever,atsuchlargeradiitheprotostellardiscsarecoolenough × G. (4)
1010cm
that they are unlikely to be fully ionized. If the ionization (cid:16) (cid:17)
fraction is sufficiently low the numerical MHD simulations HereM˙ istheaccretionrate,M themassofthecentralwhite
are not strictly applicable and the MRI which is thought dwarf, and R the disc radius, with typical values assumed.
to drive the viscosity mechanism is significantly suppressed In these numerical simulations, assumed values of βz (here
(Gammie, 1996). βz =8πP/Bz2, where P is the pressure) in the plane of the
disc are in the range 25 – 400. This implies typical vertical
fields in the range Bz ≃ 0.05 – 0.2 Beq. For a dwarf nova
2.3.3 FU Orionis outbursts in outburst this corresponds to fields of several hundred to
several thousand Gauss. There seems no obvious source for
Models for the outburstsof the pre-main sequence FU Ori-
fieldsofsuchamagnitudeonsuchaglobalscale.Moreover,a
onis stars in terms of thermally driven disc outbursts of
verticalglobalfieldthreadingadisccanbeexpectedtogive
the kind seen in CVs and SXTs are given by Clarke et al.
rise to a wind or jet from the disc (Lovelace, Romanova &
(1990), Bell & Lin (1994), and Lodato & Clarke (2004). In
Contopoulos, 1992; Pelletier & Pudritz, 1992) and a rough
order to fit the timescales of the outbursts the models re-
estimate of the mass loss of such a jet/wind is given by
quireα≃0.001−0.003. Thesevaluesaresignificantlylower
Pringle (1993) as
thoserequiredfordiscsundergoingphysicallyanalogousout-
burstsin thebinary systems, although theoutburstsin FU M˙ R
jet ∼ β−1α−1. (5)
Orisystemsseemtoneedtobemediatedbythepresenceof M˙ H z
4 A. R. King, J. E. Pringle & M. Livio
For typical values of α∼0.1 and H/R∼1/30 it is evident 5,−0.14<z/R∗ <0.14.Thevolumeaveragedstressparam-
that the typically assumed values of the globally imposed eter αis found to be≃0.005 throughout the simulation.
fieldwouldbeabletogeneratewindmasslossescomparable A discrepant note is set by Hawley & Krolik (2001)
to thedisc accretion rates. who find a value of α∼0.1 (their Fig 13) for radii R/R ≤
s
In view of all this it seems sensible to assume that the 30. They perform a global simulation of a disc around a
numerical simulations most likely to correspond tophysical pseudo–Newtonian black hole. They start with a torus at
reality are those with no net poloidal field. In the following R = 30Rs (with Rs being the Schwarzschild radius). The
we consider some representative simulations. These all use grid is in (R,φ,z) coordinates, with −10<z/Rs <10. The
full MHD, and thuscorrespond to fully ionized discs. transversemagneticfieldcomponentsaresettozerooutside
thecomputationaldomain.Theinitialfieldispoloidalalong
EarlysimulationswerecarriedoutbyStoneetal.(1996)
equal density contours, which in effect implies that there is
whoconsiderashearingbox,withverticalstructureconfined
anetB throughtheinnerhalftorusR/R ≤30(andanet
within −2H <z<2H.Theyadopt periodic boundarycon- z s
B ofoppositesignthroughtheouterhalf).Thesimulations
ditionsintheverticaldirection,essentiallyfornumericalrea- z
ran for t = 1500 (in units with GM = R = 1), which is
sons.Inthesesimulationstheseedfieldseitherhavezeronet s
about 7 orbits at R/R =10. Thus it seems likely that the
vertical field or are purely poloidal, and α is defined as the s
initial global poloidal seed field is still present throughout
time and volume averaged stress, normalized to the initial
thecomputation, in which case it is not surprising that the
mid-planepressure. The measured α is stated as <0.01 for
resultingvalueofαcorrespondsmorecloselytothosefound
most of thesimulations. Atabout thesame time, Branden-
in shearing box runswhich havea superimposed seed B .
burg et al. (1995) also considered a shearing box, but now z
Itisapparentthereforethat,exceptperhapsforthoseof
withmagneticfieldkeptverticalatthez-boundaries,though
Hawley & Krolik (2001), theoretical simulations relevant to
with zero net vertical flux.They found 0.001<α<0.005.
fully-ionizeddiscswithnoimposedverticalmagneticfieldall
Morerecentworkgivesrisetosimilarvalues.Forexam- producevalues of α<0.02, and often considerably smaller.
ple,Sanoet al.(2004) considerashearingboxwith vertical ∼
periodic boundary conditions and no vertical gravity and
examine the dependence of the saturation level of the MRI
on gas pressure. For simulations with no net vertical flux 4 THEORETICAL LIMITATIONS
they find 5×10−5 < α < 0.01, with an average value of
α=0.001. The general result in need of explanation is that for fully
ionized discs, fitting the observations appears to require
Miller&Stone(2000)considerashearingboxwithver-
α∼0.1−0.4, whereas simulations consistently yield values
ticalgravityandextendthecomputationaldomain to±5H
whichareanorderofmagnitude,ormore,belowthisvalue.
in the vertical direction. They find that the regions away
This also implies that the simulations have much smaller
from thediscplane are strongly magnetic, with β as low as
magneticfieldsthanareactuallypresent,sothatdiscstruc-
around 0.02, and term this region a ‘corona’. This coronal
turesanddissipationpatterns,aswellastimescales,arenot
region is however quite unlike the solar corona, or stellar
beingsimulatedcorrectly.Thiscouldbethereasonwhysim-
coronas, in that the field is quiescent (being strong enough
pleatmospheremodelsareunabletofittheobservedspectra
to stabilise the MRI) and is well ordered, being predom-
ofCVdiscsinoutburst,especiallythelackofBalmerjumps,
inantly toroidal. They find typically that α ∼ 0.02. This
and the ultraviolet continuum (e.g. Wade, 1984). We must
work is extended by Hirose, Krolik & Stone (2005) who
therefore ask whether the simulations are missing some es-
studytheverticalstructureofgaspressuredominatedaccre-
sential ingredient. We consider various possibilities in turn.
tion discs with local dissipation of turbulenceand radiative
transport. They have a shearing box, with vertical gravity,
and−8H <z <8H.Theyfindasimilardiscstructurewith
magnetically strong regions (‘coronas’) at large |z|, and ob- 4.1 Restrictions of scale
tain α≃0.016.
Shearingboxsimulationsmissoutonlowvaluesofazimuthal
Other workers have needed to move away from the wavenumber m. This is because the azimuthal box size is
shearing box approximation and to consider more radially typically around 2πH, and so these simulations can only
extended computational domains. Winters et al. (2003) in- handlem=0,R/H,2R/H,3R/H,.... Of course R/H is not
vestigategapformation byplanetsinturbulentprotostellar actually defined for the simulations, as in effect the limit
discs.TheyusefullMHDandthusstudyfullyionizeddiscs. R → ∞ is required for the box to be Cartesian. But the
They consider a cylindrical annulus 0.25 < R/R∗ < 3.75 magneticstructureswhichgenerallyemergehavestructures
and a vertical grid −2H < z/R∗ < 2H, and before adding up to and including the box size (see Figure 16 of Hirose,
a planet theyestablish abackground flow. They donot use Krolik and Stone2005).
aninitialverticalfieldasthisisknowntoproduceaseriesof Thenetresultofthisisstriking.Theglobal3Dsimula-
evacuated gaps in the disc (Steinacker & Papaloizou 2002; tionsbyHirose, Krolik andDeVilliers (2004) dealwith the
Nelson&Papaloizou2003).Theyuseaninitialseedtoroidal evolution of a (quite thick) torus. This is close to a black
field, and their vertical boundary conditions are periodic. hole and so not in any sort of equilibrium. Looking at the
They quote a global average value of α=0.02. fieldstructureweseeinFigure6thatjustabovethedisc(in
Nelson (2005) studiestheorbital evolution oflow mass theregion calledthecorona)thefieldisstronglyazimuthal.
protoplanets in turbulent,magnetized discs. There is no z– Andindeedin thebody ofthedisc, all m valuesareclearly
dependence,and he uses vertically periodic boundaries and present,withmostofthepoweratlowvaluesofm.Thusthe
atoroidalseedfield.Thegridiscylindricalwith1<R/R∗ < main effect of a small box is to eliminate the possibility of
Accretion disc viscosity: how big is alpha? 5
large–scale field structures, and thus transmission of power tice, periodic boundary conditions act much like rigid walls
in thespatial spectrum to low valuesof m. in that there can be no net loss of mass or magnetic flux
Similar considerations apply also on radial scales. In through them.” An attempt at circumventing this problem
Figure 5of Hiroseet al (2004) we seethat large–scale mag- was made by Miller & Stone (2000) in order to deal with
neticlinkagescan occurin theradial direction. Thusshear- thestronglymagneticdiscregionsclosetotheboundary.In
ing boxes prohibit the generation of large–scale magnetic order to reduce the limitations of the Courant condition in
structureseitherbyinversecascade(PringleandTout,1996) regions where the Alfv´en speed, v becomes unacceptably
A
orbyfootpointtwisting(Lynden-Bell,2003). Itmaybesig- large they introduced the concept of an Alfv´en speed lim-
nificant that the one computation that looked at 3D global iter.This limitation iseffected inpractice byincreasing the
simulations (Hawley and Krolik, 2001) did appear to get a momentumdensityofthefluidbyafactor[1+v2/c2 ].This
A lim
largervalueofα.Itisdifficulttobeconclusivehere,sinceas implies,ofcourse,thatnotalltheconservationpropertiesof
we remarked above the simulation probably still contained theMHDequationscanberetained.Hiroseetal.(2006)face
a net poloidal field in the relevant region. similar problems which they address by imposing a density
WenotefurtherthatthefulldisccomputationsofWin- floorandbyimposingavelocity capofaround30timesthe
ters et al. (2003) and of Nelson (2005) are not restricted gassoundspeedonthediscmidplane.Theyhaveanoutflow
to low m, but do restrict the vertical structure (periodic boundarycondition,andinlinewiththecommentsofStone
boundary conditions vertically, and Nelson has no vertical etal.,theynotethatitneedscarefulhandlingtoensuresta-
gravity).Fromthiswecanconcludethatsimplyallowingall bility.Theyalsoaddadiffusivityintheghostcells,andnote
azimuthal valuesof m to bepresent does not byitself solve thatthesignofthePoyntingfluxacrosstheboundaryisnot
theproblem. restricted.
Fromang and Nelson (2006) present global 3D models.
Their radial extent is a factor of 8, and their azimuthal
4.2 Boundary conditions
extent an angle π/4 (thus only m = 0,8,16,24,32,... are
Shearingboxsimulationshaveperiodicazimuthalandradial present). Their vertical extent is 0.3 – 0.4 times the inner
boundaryconditions.Theradialoneisphase–shiftedtotake radius, with at most 25 grid cells per vertical scale height.
account of the shear, which is acceptable within the limita- They have H/R = 0.07−0.1. They use both outflow and
tions discussed above, although it should perhaps be noted periodic vertical boundary conditions. They comment that
that the radial force is discontinuous at the radial bound- thelatterislessphysical,buthastheadvantageof preserv-
aries. The calculations of Armitage (1998), however, which ing thetotal flux of themagnetic field and thevanishingof
model a radially extended disc with no vertical structure, itsdivergence,whichisevidentlydifficulttoensurewiththe
butwithaverticallyimposedfield,suggesttheshearingbox outflowcondition. Forthelattertheyusetheapproachesof
assumptionmightalsohaveasignificanteffectbyrestricting Miller and Stone (2000). During the simulation the upper
thescale of thefield in theradial direction. layers of the disc develop very strong fields, forcing them
However, the vertical boundary condition poses quite tousean‘Alfv´enspeedlimiter’.Thisseemstoindicatethat
a different problem for attempts to represent the relevant fluxistrappedbytheboundaryconditions.Indeedtheyfind
physics realistically. Usually for the shearing box this too (theirSection4)thatthefinalstatesofthemagneticcorona
is taken to be periodic (implying a stack of accretion discs, are the same for both sets of boundary conditions. In line
rather like an old juke box). This prevents magnetic flux with otherwork they find an average effectiveα = 0.004.
from escaping. Another approach (e.g. Brandenburg et al.,
1995)assumesthatthefieldiskeptverticalattheboundary,
4.3 Convergence of the simulations
and so again one cannot have flux loops escaping. Thus in
general the vertical boundary conditions serve to suppress Many papers on the application of MRI to accretion discs
magnetic buoyancy and Parker–typeinstabilities. often do not include an explicit magnetic diffusivity and
Stone et al (1996) describe early attempts to use free so allow numerical diffusivity (at the size of the grid cells)
boundaries, and their attendant difficulties: they write “In to provide the small scale limit for the turbulence process.
principle,freeboundariesthatdonotinhibitoutgoingwaves Thus the saturation level of the turbulence (and therefore
or mass motions would be the most appropriate for model- the value of α) depends on the grid size. Most interesting
ing an astrophysical accretion disk. However, we have en- here is the paper by Sano et al (2004) who find that the
countered numerical difficulties when strong (β <1) highly saturation level of the MRI turbulence depends on the gas
tangled fields are advected across free boundaries. When pressureinthe(shearing)box.However,sinceallsimulations
a strong flux loop begins to cross the boundary, the tip find similar (too small) values for α, it seems unlikely that
is ‘snipped’ off, releasing the two ends. Magnetic tension convergenceis a major issue.
forces which previously confined the loop are now unbal- There is a possible problem here however. In hydrody-
anced, and the ends of the loop ‘snap’ straight, imparting namic turbulence, we are used to thinking that the details
a large Lorentz force to the fluid near the boundary. These at thesmallest scales donotgreatly affect thebehaviourat
forcescanproducefluidmotionsthatdisrupttheentiredisk. largescales,wheretransportpropertiessuchastheReynolds
Since strong, highly tangled fields are an unavoidable con- stressaredetermined.Butitisnotclearthatthisistruefor
sequence of the nonlinear evolution of the instability, we MHD turbulence. Schekochihin et al. (2004, 2005) discuss
have found that free outflowing boundaries cannot be used thismainlyfortheinterstellarmedium.Theyarguethatthe
to studythe long-term evolution of disks. Instead,for most magnetic field structure created by the turbulence depends
of the simulations described in this paper we adopt peri- criticallyonthePrandtlnumber(i.e.theratiobetweenmag-
odic boundary conditions in the vertical direction. In prac- netic diffusivity and viscosity) and that this in turn affects
6 A. R. King, J. E. Pringle & M. Livio
the saturation level of the dynamo. There is considerable intoregionswheretheMHDapproximationispoor,oreven
discussion of the fact that MHD turbulence gives rise to breaksdown.Thesebuoyantloopsgiveriseallkindsofcom-
inverse cascades, and therefore to magnetic field structures plicatedmagneticphenomena,includingprominences,flares
which are much larger than thetypical drivingscales. and the solar wind, many of the details of which are still
poorlyunderstood.Butitisevidentthatthefieldsandflux
loops extending outside the solar surface are quite global
4.4 The breakdown of the MHD approximation
in extent. And the radial extent of flux loops can be much
greater in stars which are rapidly rotating (for example ex-
The MHD formulation used in these simulations assumes
tendingat far as5timesthestellar radiusinABDor;Hus-
fluid velocities v ≪ c and Alfv´en velocities v ≪ c. It ex-
A
sain et al., 2002). It is well known that small-scale twisting
plicitly excludes the displacement current and so removes
of magnetic footpoints on the solar surface can give rise to
the possibility of electromagnetic waves. A good discusson
large-scale changes in the global structure of the field (Aly,
of the extension of the MHD approximation to the regime
1984; Sturrock, 1991), and numerical simulation of the be-
wherethefluidvelocities approachthespeedoflightispre-
haviour of magnetic fields in the low-β regions of the so-
sentedbyGammie, McKinney&T´oth(2003).Inparticular
lar atmosphere, driven by motions in the high-β regions, is
they note the distinction between behaviour in the MHD
an active area of research (e.g. Galsgaard & Parnell, 2005;
limit when fluid density ρ → 0 and in the vacuum case
T¨or¨ok & Kliem, 2005; Mackay & van Ballegooijen, 2006a,
ρ = 0. Thus the usual MHD formulation is unable to deal
b).Inaccretiondiscs,asstressedbyUstyogovaetal.(2000)
with regions where the densities are low enough that the
and by Lynden-Bell (2003), all these phenomena are likely
approximation ofinfiniteconductivitybreaksdown andthe
tobepresentbutdrivenmuchmorevigorouslybythestrong
distinctionbetweenρ→0andρ=0becomescritical.There
discshearandbythefactthatthediscisrotatingsorapidly
areseverenumericalproblemsimplementingitinregions of
that it is centrifugally supported (and so disc velocities are
largedensitycontrastsuchastheinterfacebetweenthesolar
around 70 percent of the local escape speed).
interiorandthesolarcorona,wherethefieldsareessentially
force–free. Tout&Pringle(1996)arguethatalthoughthedynamo
There are also significant numerical problems deal- process in the disc is likely to give rise to magnetic field
ing with magnetically dominated outflows, or Poynting– structures with predominant poloidal length-scales of order
flux dominated jets. The problem is not simply that the ∼ H, it is reasonable to assume that the interaction and
timestepsgetverysmallinnumericalsimulationswherethe reconnectionofsuchstructuresoutsidetheplaneofthedisc
densitygetssmall(andsov islarge)–thepointisthatthe can give rise to an inverse cascade, producing field struc-
A
MHDapproximationmaybreakdown.Wenotethatnumeri- turesoforder∼Randgreater.Theysuggestthatthelarge-
calstudiesofthesolarcoronaencountertheseproblems(e.g. scale poloidal fields generated in thismanner can be strong
Galsgaard&Parnell,2005;T¨or¨ok&Kliem,2005;Mackay& enoughtopoweroutflowsandjets.Ofcoursethepresenceof
vanBallegooijen,2006a,b)andconsiderablecomplicationis a magnetically driven outflow can remove angular momen-
involved in dealing with them. tumfromthedisc,andsodriveaccretion,evenwithoutafor-
malkinematicviscosityorα.Butinneitherthedwarfnovae
nor the SXTs is there evidence for such outflows. However,
even without driving an outflow, large-scale poloidal fields
5 DISCUSSION
outside the plane of the disc can provide a non-negligible
We have shown that there is a large discrepancy between contributiontoangularmomentumtransport,andhenceα,
the values of the viscosity parameter α which is required of the kind that does not become apparent from, and can-
to model observations of fully ionized, time-dependent ac- noteasilybeaddressedin,thekindofnumericalsimulations
cretion discs (Section 2: α ≈ 0.1−0.4) and those which discussedabove.Thistransportcomesaboutforthreemain
are generally obtained from numerical MHD simulations reasons.
without including a superimposed magnetic field (Section First, such a process can more easily magnetically link
3: α ≤ 0.02). In Section 4 we have described some of the disparatediscradiithanprocesseswhichrequireradialpene-
limitations inherent in the numerical simulations and have trationofdiscmaterial(cf.Armitage,1998;Fromang&Nel-
noted that most of these would indeed tend to lead to un- son,2006).Oncelinked,differentialshearwindsupthefield,
derestimating the value of α. We here discuss some other increasing themagnetic energy at theexpenseof rotational
theoreticalmeansofangularmomentumtransportinvolving energy, and therefore transfers angular momentum (cf. the
globalfieldstructureswhichlieoutsidethescopeofcurrent disc-magnetosphere interaction; Livio & Pringle, 1992).
numerical simulations. Second, even if such large-scale poloidal flux loops do
WehavenotedinSection4thatoneofthemajorrestric- not give rise to a steady outflow or wind, they may well
tionsinthenumericalsimulationsisontheglobalstructure give rise to intermittent outward acceleration of disc mate-
of the magnetic field. This suppression is driven in part by rial. As noted by Blandford & Payne (1982), once a mag-
the exigency of limited computer resources and in part by netic field line in a centrifugally supported disc makes an
the limitations imposed by the boundary conditions. Per- angle of greater than 60◦ to the vertical, material can be
hapstheclosest analogy wehavetowhatmightbehappen- centrifugallyacceleratedalongit.WhileBlandford&Payne
inginanaccretiondisccomesfromlookingatthebehaviour considered a constant field structurewith a steady outflow,
of the magnetic field at the surface of the sun. Here fluid the same physics applies equally well to intermittent field
motions below the surface, where the MHD approximation structures. Thus one can envisage a continuous process by
holdsreasonablywell,drivethegenerationofbuoyantloops which small patches of disc material are from time to time
of magnetic flux which rise up through the surface layers accelerated outwards for brief periods as and when the lo-
Accretion disc viscosity: how big is alpha? 7
calfieldconfigurationbecomesfavourable.Thediscmaterial Buat-M´enard,V.,Hameury,J.-M.,Lasota,J.-P.,2001,A&A369,
acquiresangularmomentumintheprocess,butnotenough 925
energy to be expelled, and presumably falls back onto the Cannizzo,J.K.,2001, ApJ,556,847
discatsomelargerradius.Suchaprocessleadstooutwards Cannizzo,J.K.,2001, ApJ,561,L175
transportofangularmomentumbyadirectfluxofmaterial Cannizzo,J.K.,Chen,W.,Livio,M.,1995,ApJ,454,880
Clarke,C.J.,Lin,D.N.C.,Pringle,J.E.,1990, MNRAS,242,439
moving in regions out of the disc plane. Wenote in passing
Dubus,G.,Hameury,J.-M.,Lasota, J.-P.,2001,A&A,373,251
that such a process provides a more plausible radial trans-
Frank, J., King, A.R., Raine, D., 2002, Accretion Power in As-
port process, required perhaps to explain the properties of
trophysics(Cambridge:CambridgeUniversityPress)
crystalline silicate grains in the pre-solar nebula, than try-
Fromang,S.,Nelson,R.P.,2006,A&A,457,343
ingtomixmaterialupstreaminanaccretiondisc(e.g.Gail,
Gail,H.-P.,2001,A&A,378,192
2001).
Galsgaard,K.,Parnell,C.E.,2005,A&A,439,335
Third, the global field envisaged by Tout & Pringle Gammie,C.F.,1996, ApJ,457,355
(1996) generated by an inverse cascade from the tangled Gammie,C.F.,McKinney,J.C.,To´th,G.,2003,ApJ,589,444
disc field is of necessity much weaker than the mean dy- Hartmann,L.,1998,AccretionProcessesinStarFormation,Cam-
namo field in the disc. But, as the numerical experiments bridgeUniversityPress
seem to indicate, the presence of a weak poloidal field can Hartmann,L.,Calvet,N.,Gullbring,E.,D’Alessio,P.,1998,ApJ,
serve to increase the strength of the dynamo process, and 495,385
therefore the local value of α. Thus it may be that the disc Hawley,J.F.,Gammie,C.F.,Balbus,S.A.,1995,ApJ,440,742
Hawley,J.F.,Krolik,J.H.,2001,ApJ,548,348
itself is capable of generating and sustaining, at least in a
Hirose,S.,Krolik,J.H.,deVilliers,J.P.,2004,ApJ606,1083
time-averaged sense, the kind of weak global poloidal field
Hirose,S.,Krolik,J.H.,Stone,J.M.,2006,ApJ,640,901
required to enhancethe valuenumerically estimated of α.
Hussain, G.A.J., van Ballegooijen, A.A., Jardine, M., Collier
CameronA.,2002, ApJ,575,1078
KingA.R.,Ritter,H.,1998,MNRAS,293,L42
Lasota,J.-P.,2001, NewAR,45,449
6 CONCLUSION
Lewin,W.,vanderKlis,M.,2006,CompactStellarX-raySources,
Over the last decade, thanks mainly to numerical simula- CambridgeUniversityPress.
tions, we now havea much betterunderstandingof what is Livio,M.,Pringle,J.E.,1992,MNRAS,259,23
Lodato,G.,Clarke,C.J.,2004,MNRAS,353,841
the likely driving mechanism for accretion discs. We have
Lovelace,R.V.E.,Romanova,M.M.,Contopoulos,J.,1993,ApJ,
notedherethatthereis,however,roughlyafactorof10dis-
403,158
crepancy between observational and theoretical estimates
Lynden-Bell,D.,2003,MNRAS,341,1360
of the accretion disc viscosity parameter α. We have sug-
Mackay,D.H.,vanBallegooijen,A.A.,2006a, ApJ,641,577
gested possible lines for resolving this problem. While rec-
Mackay,D.H.,vanBallegooijen,A.A.,2006b, ApJ,642,1193
ognizing that this is at best close to the limits of what is
Menou,K.,Hameury,J.M.,Lasota,J.P.,Narayan,R.,2000,MN-
currently computationally feasible, we suggest that it is es- RAS,314,498
sential to undertake fully three-dimensional, global simula- Miller,K.A.,Stone, J.M.,2000,ApJ,534,398
tions, preferably in a large enough computational domain Nelson,R.P.,2005, A&A,443,1067
that the boundary conditions have little effect on dynam- Nelson,R.P,Papaloizou,J.C.B.,2003,MNRAS,339,983
ics of the thin disc. We have also noted reasons why even Pelletier,G.,Pudritz,R.E.,1992, ApJ394,117
thismaynotbeadequate,andhavedrawntheanalogywith Pringle,J.E.,1981,ARA&A,19,137
currentattemptstounderstandthedrivingofchromospheric Pringle,J.E.,1993,in’AstrophysicalJets’,eds.D.Burgarella,M.
Livio&C.P.O’Dea,CambridgeUniversityPress,1
and coronal solar activity by subphotospheric motions. Ev-
Pringle,J.E.,Tout,C.A.,1992, MNRAS,259,604
identlytheremaybesomewaytogobeforewehaveatruly
Pringle,J.E.,Tout,C.A.,1996, MNRAS,281,291
predictivetheory of accretion discs.
Reipurth,B.,Aspin,C.,2004, ApJ,608,L65
Sano,T.,Inutsuka,S.,Turner,N.J.,Stone,J.M.,2004,ApJ,605,
321
Schekochihin, A.A., Cowley, S.C., Taylor, S.F., Maron, J.L.,
7 ACKNOWLEDGMENTS
McWilliams,J.C.,2004, ApJ,612,276
ARKacknowledgesaRoyalSociety–WolfsonResearchMerit Schekochihin, A.A., Haugen, N.E.L., Brandenburg, A., Cowley,
Award.JEPthanksSTScIforhospitality andforcontinued S.C.,Maron,J.L.,McWilliams,J.C.,2005,ApJ,625,L115
Schreiber,M.R.,Hameury,J.-M.,Lasota,J.-P.,2003,A&A,410,
support under the Visitors’ Program. We thank the referee
239
for helping to clarify thecontents of thepaper.
Schreiber,M.R.,Hameury,J.-M.,Lasota,J.-P.,2004,A&A,427,
621
Smak,J.(1998)ActaAstr,48,677
Smak,J.(1999)ActaAstr,49,391
REFERENCES
Starling, R.L.C., Siemiginowska, A., Uttley, P., Soria, R., 2004,
Aly,J.J.,1984,ApJ,283,349 MNRAS,347,67
Bailey,J.A.,1975,J.Brit.Astron.Soc.,86,30 Steinacker, A.,Papaloizou,J.C.B.,2002, ApJ,571,413
Bell,K.R.,Lin,D.N.C.,1994, ApJ,427,987 Stone, J.M., Hawley, J.H., Gammie, C.F., Balbus, S.A., 1996,
Blandford,R.D.,Payne,D.G.,1982,MNRAS,199,883 ApJ,463,656
Bonnell,I.,Bastien,P.,1992,ApJ,401,L31 Sturrock,P.A.,1991,ApJ,380,655
Brandenburg,A.,Nordlund,A.,SteinR.F.,Torkelsson,U.,1995, To¨ro¨k,T.,Kliem,B.,2005,ApJ,630,L97
ApJ,446,741 Tout,C.A.,Pringle,J.E.,1996, MNRAS,281,219
8 A. R. King, J. E. Pringle & M. Livio
Ustyogova,G.V.,Lovelace,R.V.E.,Romanova,M.M.,LI,H.,Col-
gate, S.A.,2000,ApJ,541,L21.
Wade,R.A.,1984,MNRAS,208,381
Warner, B., 2003, Cataclysmic Variable Stars, Cambridge Uni-
versityPress
Winters,W.F.,Balbus,S.A.,Hawley,J.F.,2003,ApJ,589,543