Growth of the MRI in Accretion Discs -- the Influence of Radiation Transport

0.73 MB·English

by  M. Flaig

Mon.Not.R.Astron.Soc.000,1–11(????) Printed26January2009 (MNLATEXstylefilev2.2)
Growth of the MRI in Accretion Discs – the Influence of
Radiation Transport
M. Flaig,1 R. Kissmann,1 W. Kley,1
1 Institut fu¨r Astronomie & Astrophysik, Universit¨at Tu¨bingen, Auf der Morgenstelle 10, 72076 Tu¨bingen, Germany
9
0
0
2 Accepted??.Received??;inoriginalform??

n
a
J ABSTRACT
In this paper we investigate the influence of radiative transport on the growth of the
6
magnetorotational instability (MRI) in accretion discs. The analysis is performed by
2
use of analytical and numerical means. We provide a general dispersion relation to-

] getherwiththecorrespondingeigenfunctionsdescribingthegrowthratesofsmalldis-
R turbances on a homogeneous background shear flow. The dispersion relation includes
S compressibility and radiative effects in the flux-limited diffusion approximation. By
. introducing an effective speed of sound, all the effects of radiation transport can be
h
subsumedintoonesingleparameter.Itcanbeshownthatthegrowthratesofthever-
p
ticalmodes–whicharethefastestgrowingones–arereducedbyradiativetransport.
-
o For the case of non-vertical modes, the growth rates may instead be enhanced. We
r quantify the effects of compressibility and radiative diffusion on the growth rates for
t
s thegas-pressuredominatedcase.Theanalyticaldiscussionissupplementedbynumer-
a
ical simulations, which are also used for a first investigation of the non-linear stage of
[
the MRI in gas-pressure dominated accretion discs with radiation transport included.

1
Key words: accretion, accretion discs – instabilities – magnetic fields – MHD –
v
9 radiative transfer
2
0
4
.
1 1 INTRODUCTION Nowadays it is the magnetorotational instability (Bal-
0 bus & Hawley 1991, hereafter BH) that appears to be the
9 most promising candidate. The instability is local, linear,
0 Accretiondiscsarebelievedtobepresentinavarietyofas- powerful (the growth rates being of the order of the orbital
: trophysicalobjects,namelyactivegalacticnuclei,binaryX-
v period),anditexistseveniftheinitialmagneticfieldisvery
ray sources, cataclysmic variables and protostellar systems.
i weak. Both three-dimensional local shearing-box (Hawley
X Theclassicalprobleminthetheoryofaccretiondiscsiswhy
et al. 1995; Brandenburg et al. 1995) and global disc sim-
they are accreting, or, stated in other words, what is the
r ulations (Armitage 1998; Hawley 2000; Fromang & Nelson
a mechanism which transports angular momentum outwards
2006) show that the resulting turbulence indeed transports
so that matter can spiral inwards? It cannot be ordinary
significant amounts of angular momentum outwards, lead-
molecular viscosity because it is several orders of magni-
ing to values of the α-parameter that are in the range of
tude too small to account for the observed accretion rates
α∼10−3···10−2, depending especially on the initial mag-
(Pringle 1981; Balbus 2003). In order to resolve this dis-
neticfieldconfiguration.Whencomparedtoobservations,it
crepancy, Shakura & Syunyaev (1973) introduced the phe-
seems that numerical simulations generally tend to yield α
nomenologicalα-discmodelwherealargekinematicviscos-
values that are too small, possibly due to numerical effects
ityassociatedwithturbulentmotionsinthediscisthought
(King et al. 2007). The situation now looks rather compli-
todriveoutwardangularmomentumtransport.Whatmech-
cated, since it has turned out that α does not only depend
anismwouldmakethediscsturbulent,however,stillremains
onphysicalparameterssuchasviscosityandresistivity,but
a matter of debate. It is well known that in general shear
alsoondetailsofthenumericalsetup,namelythenumerical
flows are subject to nonlinear hydrodynamical instabilities,
dissipation and resolution (Pessah et al. 2007; Fromang &
butitseemsthatacompellingcaseforhydrodynamicaltur-
Papaloizou 2007). Therefore it is very important to further
bulencehasnotbeenmadeyet.Infact,theassumptionthat
investigateandtoquantifytheimpactofthevariousphysi-
Kepleriandiscsarehydrodynamicallyunstablecanbechal-
calandnumericalparametersontheMRI-inducedturbulent
lenged both on analytical and numerical grounds (see, e.g.,
transport.
Hawley et al. 1999; Balbus & Hawley 2006). Since most ac-
cretion discs consist of ionised gas, it is thus only logical to Despite the huge amount of numerical and analytical
focus on magnetohydrodynamical instabilities instead. work that has been devoted to exploring the MRI in accre-
(cid:13)c ????RAS
2 M. Flaig et al.
tion discs, comparatively little attention has been paid to As our starting point we take the equations of ideal
the issue of radiation transport. Especially when it comes magnetohydrodynamics, namely the continuity equation
to simulations which are global in nature or at least in-
∂ρ
clude vertical stratification, there are relatively few works +∇·(ρv)=0, (1a)
∂t
that include radiation transport (like Turner et. al. 2003,
Turner2004fortheradiationpressure-dominatedcase,Kro- the equation of motion
liket.al.2007,Blaeset.al.2007foradiscwithcomparable ∂v 1 „ B2 « B·∇B
radiation and gas pressure, and Hirose et. al. 2006 for the +v·∇v+ ∇ p+ − +∇Φ=0, (1b)
∂t ρ 2µ µ ρ
0 0
gas-pressuredominatedcase).Insteadusuallyanisothermal
equation of state is used (as, for example, in the works of and the induction equation
Miller & Stone 2000; Papaloizou & Nelson 2003; Fromang ∂B
&Nelson2006).Yetfromanisothermalmodelitisnotpos- =∇×(v×B). (1c)
∂t
sible to derive the true temperature profile, nor can one be
In Eq. (1b), Φ denotes the static gravitational potential of
sure that the resulting vertical structure is correct. If, on
the central object. It vanishes when linearising the above
the other hand, one chooses an adiabatic equation of state,
equations,sowedonotneedtospecifyit.Theothersymbols
but does not include cooling, then the temperature in the
have their usual meaning.
disc will rise quite rapidly due to turbulent heating, which
In order to perform a local stability analysis, we con-
is also not very realistic. Only by including radiative trans-
sider a patch of an accretion disc that is small enough so
fer can we hope to model the energetics of MRI-generated
thatthebackgrounddensityρ,pressurepandmagneticfield
turbulence correctly and achieve a quasi-steady turbulent
B=B φˆ+B zˆcanbetakenconstant.Wedonotincludea
state. In this work, we assume that the heat generated by φ z
radialcomponentinthemagneticfield;otherwisetheback-
the turbulence is transported by radiative diffusion.
ground solution would not be time independent because an
Yetanotherpointis,thatincorporatingradiativediffu-
initiallyradialfieldwouldbeshearedintoatime-dependent
sionmightwellyieldunexpectedconsequences.Forexample
azimuthal field. For a Keplerian disc, the background flow
it has turned out recently that the migration of planets in
is locally v = −3Ωrφˆ, with Ω the angular orbital fre-
aprotoplanetarydiscissignificantlyinfluencedbyradiative 2
quency. We impose plane-wave axisymmetric perturbations
effects(see,e.g.,Kley&Crida2008).Itisnotunreasonable
δρ,δv,δB,δpproportionaltoexp[ı(k r+k z)+σt],yielding
to expect that radiative diffusion might have other impor- r z
the following linearised equations:
tant influences, for example when considering the growth
of planetesimals out of dust in a turbulent disc (this has σδρ=−ρık·δv, (2a)
been pointed out recently in Brandenburg 2008). For these „ «
3 ık δB·B
reasons, including radiation transport should be considered σδv− Ωδv φˆ+ · C2 δρ+
2 r ρ eff µ
imperativeontheroadtowardsbuildingmorerealisticmod- 0
els of accretion discs. − ıkzBz δB+2Ωzˆ×δv=0, (2b)
The plan of our paper is as follows: In Sec. 2 we de- µ0ρ
rive the dispersion relation for axisymmetric MRI modes „ 3Ω «
σδB=−ık·δvB+ık B δv− δv φˆ ; (2c)
in a radiative, gas-pressure dominated accretion disc in the z z 2σ r
flux-limited diffusion approximation and calculate the cor-
wherewehaveintroducedtheeffective sound speed C via
responding eigenfunctions and MRI growth rates. In Sec. 3 eff
thedefinitionC2 ≡δp/δρ.NotethatEq.(2c)ensuresthat
we describe numerical simulations of single MRI modes us- eff
the divergence-free constraint k·δB=0 be satisfied.
ing two different numerical schemes and demonstrate their
AftereliminatingδρandδBfromthelinearisedmomen-
excellent agreement with the analytical prediction. We also
tumequation(2b)bytheuseof (2a)and(2c),theremaining
investigatethedependenceofthesaturationleveloftheMRI
systemofthreeequationscanbereadilysolved.Bydefining
on the strength of radiative diffusion in a series of 3D local
a wavenumber K via
shearing box simulations. We conclude in Sec. 4, by dis-
cussingsomeconsequencesofouranalysisandprovidingan Kδv ≡k·δv (3)
r
outlook for further research. In the appendix we relax the
assumption that gas pressure dominates radiation pressure (whichmeansthatK basicallyconstitutesameasureofthe
andprovideageneraldispersionrelation,includingboththe compressibility of the perturbations), we can write the re-
radiation dominated and the gas-pressure dominated case. sulting dispersion relation in a form that closely resembles
the incompressible dispersion relation of BH:
„ k K «k2
2 LINEAR ANALYSIS σ˜2− r σ2 σ˜2+Ω2σ˜2−4Ω2k2v2
k2 k2 z Az
z
2.1 Dispersion relation
−2Ωσv v k K =0; (4)
Aφ Az z
We begin our analysis by deriving the MRI dispersion √
wherek=|k|,v =B/ µ ρdenotestheAlfv´enspeedand
A 0
relation with radiation transport included in the one- σ˜2 ≡σ2+k2v2 .Theform(4)ofthedispersionrelationex-
temperature flux-limited diffusion approximation (see Lev- z Az
plicitlydemonstratesthedependenceonthecompressibility.
ermore&Pomraning1981).Thismeansweassumethatgas
The expression for K turns out to be
andradiationareatthesametemperatureandthatthera-
diation energy is small compared to the internal energy of K = σ˜2σ2kr−2Ωσkz3vAφvAz , (5)
the gas, i.e. the disc is gas-pressure dominated. σ˜2(σ2+k2C2 )+σ2k2v2
z eff z Aφ
(cid:13)c ????RAS,MNRAS000,1–11
Growth of the MRI in Accretion Discs – the Influence of Radiation Transport 3
where σ is one of the roots of Eq. (4). In contrast to what
onemayexpectfromitsdefinitioninEq.(3),Eq.(5)shows
K does not couple the modes δv and δv . In the incom-
r φ
pressiblelimitC →∞,wehaveK =0andindeedrecover
eff
the dispersion relation of BH.
The value of C is determined by the internal energy
eff
equation. In the one-temperature approximation suitable,
e.g., for protoplanetary discs, it reads:
∂e
+∇·(ev)=−p∇·v−∇·F, (6)
∂t
withtheradiativefluxvectorF intheflux-limiteddiffusion
approximation
∇T aT4 λc
F =−eD , with D≡4 . (7)
T e κρ
Here,κistheopacity,cthespeedoflightandλthefluxlim-
iter, which depends on the flux-limited diffusion model one
decides to use. In the linear growth phase, λ is a constant,
which is due to the fact that we consider a uniform back-
ground. This and the uniformity of the other background
quantitiesallowsusingDasaparameterforourstudy,which
willbedesignatedasthe“coefficientofradiativediffusion”.
Linearising(6)andusingp=(γ−1)e,aswellastherelation
Figure 1. Effective adiabatic index γ as a function of
eff
δT/T =δp/p−δρ/ρ we obtain km2axD/Ω, where kmax is the wavenumber of the fastest grow-
σδp=−ρC2 ık·δv; (8) ingmode.Inthisplot,theazimuthalAlfv´enspeedistakenequal
eff totheisothermalsoundspeed,v2 =p/ρ(theglobalfeaturesof
Aφ
where Ceff is given by the following expression: thisplotarenotsensitivetothisparticularchoice).
p γ+k2D/σ
C2 =γ , with γ ≡ . (9)
eff effρ eff 1+k2D/σ
Blaes & Socrates (2001) considered the case of a
We see that in the linear regime, the influence of radiative
radiaton-dominated disc where matter and radiation inter-
diffusion can be described by introducing an effective adi-
act only by momentum exchange and are, therefore, at dif-
abatic index γ . Without radiative diffusion (D = 0), we
eff ferent temperatures, while we consider the opposite limit
haveγ =γ,andC isjusttheadiabaticgassoundspeed.
eff eff wherematterandradiationarecloselycoupledviaemission
With D increasing, the effective adiabatic index γ is re-
eff andabsorptionprocesses,andarethereforeatthesametem-
ducedtowardstheisothermalvalueγ =1atD=∞.This
eff perature.Intheappendixamoregeneraldispersionrelation
reduction happens mostly in the regime where the nondi-
isprovided,fromwhichboththeBlaes&Socratesdispersion
mensional diffusion coefficient k2D/σ ∼ 1, as can be seen
relation and the one given in Eq. (4) can be derived.
fromFig.1,wheretheeffectiveadiabaticindexisplottedas
a function of D. Thus, in the linear regime, the only effect
ofradiativediffusionistomakethephysicalsituationmore
isothermal. 2.2 Eigenfunctions
Inordertoaidcomparisonwiththedispersionrelations
The corresponding eigenfunctions in terms of the radial ve-
derived by other authors, we define (following Blaes & Bal-
locityperturbationδv aregivenbythefollowingformulae:
bus 1994; Blaes & Socrates 2001) r
ω≡ıσ, (10a)
K
δρ=−ıρ ·δv , (12a)
ω˜2 =−σ˜2, (10b) σ r
Dms ≡ω2−kz2(Ce2ff +vA2φ), (10c) δvφ = 2Ωσk2 „k2σσ˜22 −krK«·δvr, (12b)
k2 z
DBH ≡ k2ω˜4−Ω2ω˜2−4Ω2kz2vA2z; (10d) δv = K−kr ·δv , (12c)
z z k r
z
and write Eq. (4) in an alternative way as
K
δp=−ıρC2 ·δv , (12d)
k2 eff σ r
D D −k2v2 v2 (k2ω˜2+3k2Ω2)− rω˜2ω4 =0. (11)
ms BH z Aφ Az z kz2 δBr =ıkzσBz ·δvr, (12e)
Intheappropriatelimitsourdispersionrelationisconsistent
k B „k2σ˜2 k K Ω2« KB ff
with that of other authors, like for example the two fluid δB =ı z z − r −3 − φ ·δv ,
φ 2Ω k2σ2 k2 σ2 σ r
dispersion relation of Blaes & Balbus (1994) in the limit of z z
(12f)
D→0andzerocoupling,ortheresistivedispersionrelation
of Sano & Miyama (1999) in the limit D → 0 and zero δB =−ıkrBz ·δv . (12g)
resistvity. z σ r
(cid:13)c ????RAS,MNRAS000,1–11
4 M. Flaig et al.
Using once more the relation δT/T = δp/p−δρ/ρ, we can
derive the following expression for the temperature pertur-
bation:
δT γ−1 δρ
= . (13)
T 1+k2D/σ ρ
Assaid,radiativediffusiontendstomakethesituationmore
isothermal, and thus for D → ∞ we have δT → 0, as is to
be expected.
2.3 Change of growth rates due to compressibility
and radiative diffusion
Inordertoquantifytheeffectsofafinitecompressibilityand
radiative diffusion on the growth rates, we first look at the
verticalmodes(k =0)whichhavethelargestgrowthrates.
r
Bydifferentiatingthedispersionrelation(4)withrespectto
C−2,wefindforthechangeofthegrowthrateσwithrespect
eff
to a change in the effective sound speed:
1 „ dσ « 2k2v2 /σ˜2
=− Az v2 . (14)
σ dC−2 1+2σ˜2/Ω2 Aφ
eff Ceff=∞
From this result we conclude that for vertical modes the
effect of a nonzero compressibility leads to a decrease of
the growth rates in the presence of a nonzero azimuthal
field. For the fastest growing mode, which, for C = ∞,
eff
has a growth rate σ = 3/4Ω with wavenumber k =
√ max max
( 15/4)Ω/v , this means that the change ∆σ in the
Az max
growth rate as compared to the incompressible case C =
eff
∞ is approximately: Figure2. SketchoftheMRIeigenmodes.Plottedisthedensity
perturbationδρ(grayscale)andthevelocityperturbationδv(ar-
∆σmax ≈−1vA2φ. (15) rows) in the r-z plane. The plots on the left show the case of a
σ 5C2 verticalmodeinthepresenceofabackgroundazimuthalmagnetic
max eff
field;where(a)correspondstotheincompressiblecase,while(b)
The relative dampening of the growth rates is, thus, in this
isforfinitecompressibility.Theplotsontherightshowthecase
case of the order of O(vA2φ/Ce2ff). ofanonverticalmodewithkr =kz wherenoazimuthalmagnetic
Curiously, if we consider non-vertical modes (kr (cid:54)= 0) field is present. Here, (c) corresponds to the case of maximum
and a vanishing azimuthal field, we discover the opposite. compressibilityand(d)totheincompressiblecase.
In this case we find by an analogous analysis:
1 „ dσ « k2 σ˜2σ2/k2Ω2
σ dC−2 = 2kr21+2k2σ˜2z/k2Ω2. (16) Forexample,recentresultsshowthatinzero-netfluxshear-
eff Ceff=∞ z z ing box calculations this depends critically on the value of
The corresponding shift of the maximum growth rate be- the magnetic Prandtl number (Fromang et al. 2007).
comes: To summarise, radiative diffusion changes nothing fun-
damentalregardingthegrowthoftheMRI.Thisisbecause
∆σ 27k2 v2
max ≈ r A . (17) the radiative effects enter only via the effective adiabatic
σ 240k2C2
max eff index γ and, thus, change the compressibility, but noth-
eff
This effect should not be considered too important, be- ing else. This statement, that the growth of the MRI under
cause the growth of the MRI will be dominated by the the influence of radiative diffusion can be characterised by
fastest growing modes, which are the vertical ones. In both an effective sound speed, remains true also if we drop the
casesconsidered,radiativediffusioncontributesafractionof one-temperature approximation (cf. the appendix).
(γp/ρ)/C2 −1=γ/γ −1tothetotalshiftinthegrowth In gas-pressure dominated discs, the growth rates will
eff eff
rates. never be dramatically changed by radiative diffusion, since
The critical wavelength, which is obtained by setting thevariationinC2 coversonlyafactorofγ.Thisisdiffer-
eff
σ =0 in the dispersion relation, is the same as in the non- entfromthesituationencounteredinaradiationdominated
radiative case. This means that radiative diffusion will gen- disc,wheretheanalogofC2 mayvarygreatly.Inthepres-
eff
erally not affect the regime in which the MRI may possibly ence of a strong enough azimuthal field, the growth rates
operate.Oneshouldbearinmindhowever,thatnotthelin- maythusbeseverelyreduced(cf.theappendixandthepa-
earanalysis,butonlynumericalsimulations,canprovideus per of Turner et al. 2002).
withtheanswersif,inagivensituation,theturbulencethat In order to understand the change of the growth rates
isinitiatedbytheMRIwillbesustainedforalongperiodof duetoradiativediffusioninaqualitativemanner,letuscon-
time and how effectively it transports angular momentum. sidertwocases:First,thecaseofaverticalmode(k=k e )
z z
(cid:13)c ????RAS,MNRAS000,1–11
Growth of the MRI in Accretion Discs – the Influence of Radiation Transport 5
in the presence of an azimuthal field. In the incompress- e.g., Balsara & Spicer 1999; Londrillo & Del Zanna 2000)
ible limit (C → ∞), the motion of the fluid is confined by which we assure the solenoidality of the magnetic field.
eff
to the plane perpendicular to the perturbation wavevector, This method uses the hyperbolic fluxes to compute electric
so that δv = 0 (Fig. 2 (a)). When the compressibility is fields on a staggered grid. These are then used to evolve
z
nonzero, the Lorentz force due to the azimuthal magnetic the magnetic induction, the components of which are also
fieldcausesafluidflowintheverticaldirection(Fig.2(b)). given on a (different) staggered grid. The stability of the
The higher the compressibility (and the smaller therefore code and its capability to resolve steep gradients without
the gas pressure), the more vertical will the velocity vector introducing artificial oscillations have been proven, e.g., in
become. This in turn makes the buildup of the magnetic Kissmann (2006).
field less effective, resulting in a smaller growth rate. Ra- Being interested in the evolution of the gas under the
diativediffusionincreasesthecompressibility,andtherefore influenceofradiationtransportwealsohavetodealwiththe
decreases the growth rate. corresponding source terms to the energy equation. These
Next consider the case of a nonvertical mode (k (cid:54)=0), were implemented using two alternative approaches. In the
r
with no background azimuthal field. If we now start with first approach we included the diffusive radiation transport
the limit of maximum compressibility (C = 0), we have explicitlyinthescheme,whereasweusedanimplicitsolver
eff
δv = 0 (Fig. 2 (c)), since the z-component of the Lorentz for the second method. This dual approach is motivated by
z
force vanishes. When now increasing C , the gas pressure the fact that we are interested in being able to cross-check
eff
will act to push the velocity vector into the plane perpen- the two different methods.
dicular to the perturbation wavevector, this effect becom- Fortheimplicitimplementationofthesourcetermswe
ing stronger and stronger as the compressibility is further splittheevolutionoftheradiationfromtheevolutionofthe
decreased (Fig. 2 (d)). Therefore, this time increasing the hyperbolic part of the system of equations:
compressibilitymakesthebuildupofthemagneticfieldmore „ «
∂e
effective,whichmeansthatthegrowthrateincreasesdueto e(∗) =e(n)+ (19)
∂t
radiative diffusion. hyp
„ «
When considering the general case of a nonvertical e(n+1) =e(∗)+ ∂e (20)
mode in the presence of a nonzero azimuthal field, both ef- ∂t
rad
fects are present and the result will be either an increase
whereweevolvetheradiationtransportstepusingaparallel
oradecreaseofthegrowthrate,dependingonthestrength
ω-Jacobimatrixsolver.Theω-Jacobimethodisverysimilar
of the azimuthal field and the direction of the perturbation
tothemethodofsuccessiveoverrelaxation(SOR).Ithasthe
wavevector.
advantage that it is simple to implement and very easy to
parallelise. We have also a multigrid method available, in
case that the ω-Jacobi solver should prove insufficient for
3 NUMERICAL SIMULATIONS
high-resolution simulations.
In this section we investigate the growth of the MRI using Having an implicit solver for the radiation transport
the appropriate numerical tools. Here we especially aim at is of special importance with regard to future applications.
a verification of these numerical tools by a comparison to Thatisforhighresolutionsimulationswithstrongradiation
the analytical results derived above. After discussing the transportthetime-stepofanexplicitschemewouldbecome
numericalschemewewillpresenttheresultsofthenumerical prohibitivelysmall,whereasanimplicitsolverwouldonlybe
simulations of the MRI. limited by the Courant condition (see Courant et al. 1928)
for the hyperbolic part of the system.
Our second implementation is nonetheless of explicit
3.1 Numerical approach form, because we wanted to have a completely different
method at hand to be able to compare it to the implicit
Theinvestigationoftheanalyticalresultsusinganumerical
scheme for more general simulations for which no analyti-
method is clearly not very demanding in this case owing to
cal solution is available. For this purpose we extended the
the fact that the analytical model is restricted to the lin-
flux-function for the energy equation by:
ear growth phase. Here we decided, however, for a scheme,
whichiswelladaptedtohighMach-numberturbulencesim- e „∂T ∂T ∂T «
FRad =− D e + e + e (21)
ulations, because we are also interested in more general ac- T ∂x x ∂y y ∂z z
cretion disc simulations including radiation transport. One
Sincethefluxfunctionsareonlyneededatthecellfacesthe
of the most important constraints for a numerical scheme
occurringgradientsareeasilycomputedascanbeseenfrom
used for compressible turbulence simulations is the correct
the example of the x-component:
handling of sharp discontinuities. Here we use a second or-
„ «
der finite volume scheme based on the work by Kurganov e T −T
FRad =− D i+1,j,k i,j,k (22)
et al. (2001) for the hyperbolic part of the system of equa- x,i+1/2,j,k T ∆x
tions. This is a central conservative scheme for the solution
This scheme proved to be stable for:
of equations of the form:
∂u +∇·F(u)=0. (18) ∆t< ∆x2 ρT/e (23)
∂t 8D γ−1
It does not require a Riemann solver and a characteris- In the following section we will show that both implemen-
tic decomposition. This scheme was extended by a con- tations for the radiation transport are well-behaved, yield-
strained transport description for the magnetic field (see, ingcorrectresultwhenappliedtotheanalyticalproblemat
(cid:13)c ????RAS,MNRAS000,1–11
6 M. Flaig et al.
D=∞
Ω D=1.6·10−5
/
σ D=8·10−6
D=4·10−6
D=2·10−6
D=10−6
D=0
kzvAz/Ω
Figure 4. Comparison of the numerical growth rates (black
crosses)totheanalyticalmodelresults(fulllines).Hereweshow
the data from the simulations with the explicit numerical im-
plementation of the radiation transport terms. All simulations
haveBφ=25Bz andaverticalmagneticfieldcorrespondingtoa
plasmabetaβ=400.ThecurveforD=∞hasbeeninsertedto
Figure 3. MRI growth rates obtained with the implicit solver.
showthemaximumpossiblechangeinthegrowthrates.
Plotted are the analytical solutions for various cases, together
with corresponding data points obtained by numerical simula-
tions.Thecurvesinthemainpartofthegraphicareforvertical
modes(kr =0),wherethegrowthratesarereducedbyradiative the simulations with the implicit scheme, the full radiation
diffusion.Fromtoptobottom,theyrepresentthefollowingcases,
transport was done as prescribed by Eq. (7), with λ = 1/3
Bφ = D = 0; Bφ = 20Bz,D = 10−5; Bφ = 20Bz,D = 10−4; and constant opacity κ. The data points nicely match the
Bφ = 20Bz,D = 1. The vertical magnetic field in terms of the
plasma beta, β =2µ0p/B2, is chosen such that β =400, which analytical prediction, thereby confirming both the validity
meansthatvAz/pp/ρ(cid:39)0.1.Theinsetshowssimulationresults of the numerical scheme and the analytical calculation. In
for nonvertical modes with no azimuthal magnetic field, where particular it is evident that both implementations of the
radiativediffusionincreasesthegrowthrates.Wechosekr =kz, radiation transport are in good agreement to the analytical
p
Bφ = 0 and Bz such that vAz = p/ρ. From bottom to top, predictions.Therefore,wecanusebothmethodstosimulate
we have D = 0, D = 10−5, D = 10−4, D = 1. Although the morecomplexsituationsinthefuturegivingusthegoodop-
differenceinthegrowthratesisonlysmallforthelattercase,it portunitytohavetwodifferentmethodstobeappliedtothe
isaccuratelyreproducedbyournumericalcode.
same problem. Thus, it will be easier to decide if observed
effects might be due to some numerical problem.
hand.Thus,theycanbeusedforfuturesimulationsofmore
complex radiation transport phenomena.
3.3 Saturation level
3.2 Growth rates In order to study the impact of radiation transport on the
saturationleveloftheMRI,weswitchtoa3Drepresentation
For our numerical simulations, we use the shearing box ap-
of the shearing box with a computational domain [-0.5,0.5]
proximation(see,e.g.,Hawleyetal.1995).Inordertocom-
x[0,4]x[-0.5,0.5].Concerningthevaluesofρ,p,Ω weusea
pare the analytical growth rates of the MRI to the results
similar setup as in the 2D simulations described in the pre-
from our numerical code, we perform simulations of single
vious section. We choose an initially vertical magnetic field
MRI eigenmodes. Since only variations in the radial and in
corresponding to a plasma beta of β = 400 in all simula-
thez-directionshavetobetakenintoaccount,itissufficient
tions.Insteadofprescribingeigenfunctions,weinitialisethe
to use a 2D model. We set up the problem using the eigen-
problem with random velocity and pressure perturbations
functions as given by (12). The size of the computational
of order 10−6.
domain is [-0.5,0.5] x [-0.5,0.5] (x- and z-direction respec-
The strength of the angular momentum transport is
tively) with a typical resolutionof 64x256 cells. Initiallywe
described by the α-parameter (Shakura & Syunyaev 1973).
set Ω = 10−3, density ρ = 1 and pressure p = 10−6. The
Inoursimulations,itismeasuredaccordingtothefollowing
simulations were performed for different wavenumbers and
prescription:
forvariousstrengthsofmagneticfieldandradiativediffusion
coefficient (for the explicit values see in the figures). α=(cid:104)T +T (cid:105)/p , (24)
R M 0
The results are plotted in Figs. 3 and 4 for the im-
plicit and the explicit implementation of the source terms, where p is the initial gas pressure, T = ρδv δv denotes
0 R r φ
respectively, together with the corresponding analytical so- the Reynolds stress and T = −δB δB is the Maxwell
M r φ
lution. In the simulations with the explicit scheme a con- stress.Theangularbrackets(cid:104)···(cid:105)indicateaspatialaverage.
stant diffusion coefficient was used [cf. Eq. (21)], while in We also use the dimensionless stresses α and α defined
R M
(cid:13)c ????RAS,MNRAS000,1–11
Growth of the MRI in Accretion Discs – the Influence of Radiation Transport 7
Table 1.Time-andspace-averagedquantitiesfrom3Dshearing Table 2.3Dshearing-boxsimulationsperformedwithZEUS.In
boxsimulationsperformedwithourcode.Inallofthesimulations someofthesimulations,anon-zeroconstantB wasused,leading
φ
B = 0. In the first three lines, only the resolution is varied, to a higher saturation level. Averages were taken from 10 to 40
φ
whilethelastthreelinesshowthedependenceonthestrengthof orbits.
radiativediffusion.Averagesweretakenover10to40orbits.
Resolution D Emag/p0 α αM αR Resolution D BBφz Emag/p0 α TM TR
32x64x32 0 0 0.398 0.248 0.210 0.038
32x64x32 1 0.400 0.278 0.248 0.044 32x64x32 1 0 0.245 0.146 0.123 0.023
128x256x128 1 0.505 0.299 0.279 0.036 32x64x32 0 20 2.557 1.125 0.983 0.141
64x128x64 1 0.337 0.226 0.195 0.031 32x64x32 1 20 1.542 0.577 0.498 0.079
64x128x64 0 0 0.552 0.335 0.288 0.047
64x128x64 0 1.007 0.719 0.599 0.119
64x128x64 1 0 0.392 0.218 0.187 0.031
64x128x64 10−5 0.358 0.244 0.210 0.034
64x128x64 0 20 2.474 1.137 1.006 0.131
64x128x64 1 0.337 0.226 0.195 0.031
64x128x64 1 20 2.315 1.075 0.938 0.137
as as in the runs with our code. Despite these issues, all runs
unambigouslyshowthatthesaturationleveldecreaseswith
α =(cid:104)T (cid:105)/p , α =(cid:104)T (cid:105)/p . (25)
R R 0 M M 0 increasing radiative diffusion.
A plausible explanation for the phenomenon of the re-
Note that with this definitions in general we will have
ductionofthesaturationlevelduetoradiativediffusioncan
α (cid:54)= α +α . Due to the fact that in our simulations the
R M
be given as follows: In the case of an initial magnetic field
magnetic field posesses a net flux, the resulting α is quite
large, α(cid:38)0.1. withnonzeronetflux,theturbulenceinthesaturatedstate
isstillpartiallypumpedbythetwo-channelmode(theverti-
In order to check if there exists a trend for the turbu-
calmodewiththelongestwavelengththatfitsintothebox).
lent activity with different resolution, we performed three
The channel mode reappears every few orbits; a manifesta-
runs where the resolution was successively doubled. Fig. 5
tion of this are the oscillations that can be seen in the time
shows the time evolution of the magnetic energy and the
evolutionofthemagneticenergyandtheα-parameter(see,
α-parameter for this runs. In Table 1 long-term averages
forexampleSano&Inutsuka2001).InSec.2.3wehaveseen
of this quantities as well as of the Maxwell and Reynold
that radiative transport tends to decrease the growth rates
stresses can be found. The results displayed in Fig. 5 and
of vertical modes by effectively increasing the compressibil-
Tab.1indicatenonettrendforthechangeofthesaturation
ityofthefluidandthusmakingiteasierfortheLorentzforce
level with resolution (the first three lines in the table).
topushthevelocityvectorsoutofther-φplane.Thereforeit
Aswehaveseen,theanalyticalcalculationpredictsthat
istobeexpectedthatradiativetransport,bydecreasingthe
radiativediffusionmodifiestheMRIgrowthrates.Theques-
growth rate of the recurrent channel mode, acts to reduce
tion if the saturation level of the MRI-induced turbulence
the turbulent activity.
isalsochangedbyradiativeeffectscanonlybeansweredby
numericalsimulations.Wehaveperformedsimulationswith
varyingradiativediffusioncoefficient,wherewepickedares-
3.4 Temperature distribution
olutionof64x128x64.Theresultsshowamoderatedecrease
of the turbulent activity with inreasing radiative diffusion For the full 3D simulation the influence of the radiation
coefficient D, see last 3 entries in Tab. 1. A qualitative ex- transportcanalsobevisualisedbythetemperaturedistribu-
planation for this is given below. tion.ThisisshowninFig.7forasimilarsetupasdescribed
Wehavealsoperformedadditionalsimulationswiththe in the previous section. Here we used an extent of the com-
ZEUScodeinordertocheckthisresultusingadifferentnu- putational domain of [-0.5,0.5] x [0,6] x [-0.5,0.5] with the
mericalscheme.Weusedthesamesetupconcerningρ,p,Ωas same spatial resolution of 64x128x64. The z-component of
in the simulations done with our code. The strength of the the magnetic induction was initialised with β = 800, and
constant verticalmagneticfieldwastakenascorresponding we chose B =20B . We followed the evolution for several
φ z
to a plasma beta of β = 800. We performed simulations orbits. In Fig. 8 we show the distribution in the density -
with and without azimuthal magnetic field B , using dif- thermal energy plane for one simulation without and one
φ
ferent resolutions (see Table 2). In Fig. 6, the magnetic en- with radiative transport. Apparently the plasma gets more
ergy is plotted for the high resolution simulations. As has isothermal when radiation transport is present – for a fully
already been stated in the introduction, it is known that isothermalsimulationthethermalenergywoulddependlin-
the outcome of an MRI simulation does depend on numeri- early on the density, thus, yielding a line in this plot. This
cal issues such as the resolution or the numerical magnetic factalsobecomesapparent,whencomparingtheresultsfor
Prandtl number of the code (Fromang & Papaloizou 2007). the temperature distribution function.
Indeed, if we compare for example the simulations in lines This is done in Fig. 7 for the same simulations. Here,
3 and 4 of Tab. 1 with the corresponding simulations done the temperatures are normalised to the initial temperature
with the ZEUS code, lines 5 and 6 of Tab. 2, the ratio of (that is, initially we had a delta-peak at T = 1). The shift
the saturation levels in the case of the ZEUS code turns of the distribution functions with regard to the initial tem-
out to be about a factor of two smaller than in the simu- perature results from the heating by (numerical) viscosity
lations done with our code. We have checked that this fact and resistivity. Obviously the distribution for the disc with
remainstrueevenifwechoosethesameinitialplasmabeta radiationtransportisconsiderablynarrower.Note,however,
(cid:13)c ????RAS,MNRAS000,1–11
8 M. Flaig et al.
32x64x32
64x128x64
0 128x256x128
p
/
ag α
m
E
32x64x32
64x128x64
128x256x128
Figure5. Resolutiondependenceoftheturbulentactivityonthenumericalresolution.Theleftandrightpanelsshowthetimeevolution
oftheperturbedmagneticenergyandthealphaparameter,respectively.Incaseofthealphaparameterthecurveshavebeensmoothed
over1orbit.Whenconsideringlong-termaveragesofthesequantities,thereappearstobenotrendwithvaryingresolution.
T T
f f
T T
Figure7.Temperaturedistributionfortwo3Dsimulationswithoutradiationtransport(left)andwithradiativetransport,D=4·10−5
(right).Inbothplots,thedistributiontotheleftcorrespondstot=3.5orbits,whiletherightcorrespondstot=4.1orbits.Inthecase
withradiativetransport,thewidthofthedistributionfunctionissmallerandtheheatingoccurslessrapidly.
that also for this case the temperature after several orbits 4 DISCUSSION
differs markedly from the initial temperature. This conlu-
sion is also strengthend by the fact that both distributions Weanalysedtheinfluenceofradiationtransportonthemag-
arecenteredatnearlythesametemperature.Thisisdueto netorotational instability. We started by investigating the
thefactthatthegasisstillheatedbyviscosityandresistiv- linear growth of small disturbances on an underlying con-
ity. In contrast to the case without radiation transport the stant shear flow. The system under investigation is a mag-
local temperature enhancements are smoothed out by the netised shearing box, where we explicitly retained the com-
radiation transport. pressibilityofthegas.Additionallywetookradiationtrans-
port via flux-limited diffusion into account. For this system
wederivedageneraldispersionrelation,whichintheappro-
priate limits can be reduced to several dispersion relations
known from the literature.
(cid:13)c ????RAS,MNRAS000,1–11
Growth of the MRI in Accretion Discs – the Influence of Radiation Transport 9
h h
et et
ρ ρ
Figure 8. Distribution in the density- thermal energy plane. Values are shown in normalised units. Here we show the distribution
function for a simulation without radiative transport (left) and one where we used D = 4·10−5 (right). Both snapshots are taken at
timet=4.1orbits.
predictionnotonlyconfirmedtheaccuracyoftheemployed
numericalschemesbutalsothecorrectnessoftheanalytical
solution.
This analysis was then further extended by numerical
simulations, which also took the non-linear phase of the in-
stability into account. These simulations were done for a
more general case, where several modes of the instability
where excited in a fully three dimensional setup. There we
concentrated on the resulting saturation levels of the tur-
bulence and on the temperature distribution, which is also
influenced by the radiation transport.
Fromouranalysisitisapparentthattheinfluenceofthe
radiationtransportisarichphenomenon.Itnotonlyweak-
ens the magnetorotational instability, but can also enhance
it in a certain parameter regime. This shows that radiation
transport – even in the one-temperature flux-limited diffu-
sionlimit–doesnotsimplyactasadissipativeprocesslike
viscosity and resistivity. Another difference between these
processes is that the latter can completely suppress the in-
stability, whereas radiation transport can only weaken it.
This is due to the fact that a diffusive radiation transport
Figure 6. 3D shearing box simulations performed with ZEUS only smoothes the temperature, whereas the other dissipa-
at resolution 64x128x64. Shown is the perturbed magnetic en- tive processes directly act onto either velocity or magnetic
ergy normalised to the initial gas pressure. The solid, dashed, fieldfluctuations.Thus,wehave,forstrongradiationtrans-
dot-dashed and dotted curves correspond to the combinations port,approximatelyanisothermalstateoftheplasma.This,
Bφ = 0/D = 0, Bφ = 0/D = 1, Bφ = 20Bz/D = 0 and however,stilldiffersfromtheuseofanisothermalequation
Bφ=20Bz/D=1,respectively.Inthelinearphase,theobtained of state with the difference arising from the fact that when
growthratesareonlyafewpercentsmallerthanthegrowthrate
using an energy equation with γ differing from one, we still
ofthefastestgrowingmodethatfitsinthebox,whichmeansthat
haveheatingduetodissipativeprocesses(eitherphysicalor
the initial growth is dominated by the fastest growing mode, as
numerical).Thus,theeffectofradiationtransportisjustto
istobeexpected.
smooth out the local peaks of heat appearing due to local
dissipation.
Here, we placed special emphasis on the one- This picture will probably become different, when in-
temperature flux-limited diffusion approximation suited for vestigating the MRI in an open box instead of the periodic
accretion discs not dominated by the radiation. The corre- shearingbox.Thefuturetaskwill,thus,betheinvestigation
sponding dispersion relation can also be derived from the ofastratifiedshearingbox,wherewecanallowforoutgoing
general relation given in the appendix. From this we com- radiation in the halo of the disc. This way even a nearly
puted the growth rates for the linear case, which we then isothermal core of the disc might occur, which might then
also compared to the results from numerical simulations. becomparedtosimulationsusinganisothermalequationof
Theexcellentagreementbetweensimulationsandanalytical state.Here,however,wehaverestrictedourselvestothein-
(cid:13)c ????RAS,MNRAS000,1–11
10 M. Flaig et al.
vestigation of the local shearing box, since only this can be From these linearised equations we find after some algebra
approached by analytical means. δp = (1+αP+αF)`γ+ 43Ee´− 43Ee(1−3αP)(1+αF)p,
δρ 1+α +α +4α E(1+α ) ρ
P F Pe F
Acknowledgements (A7a)
This research has been supported in part by the Deutsche δE = 1+4αP`γ−1+ 43Ee´ E. (A7b)
ForschungsgemeinschaftDFGthroughgrantDFGForscher- δρ 1+α +α +4α E(1+α ) ρ
P F Pe F
gruppe 759 “The Formation of Planets: The Critical First
Inconclusion,letuslookatthelimitsofgas-pressuredomi-
Growth Phase”. Computational resources wereprovidedby
nated and radiation-pressure dominated discs, respectively.
theHighPerformanceComputingClusteroftheUniversity
In the limit α →∞, E/e→0, suitable for a gas-pressure
ofTu¨bingenandintheprojecthbo25bytheForschungszen- P
dominated disc, we get
trumJu¨lich.Wethanktheanonymousrefereeforproviding
valuable suggestions that helped to improve the paper. δp = γ+4EeαFp, δE =0; (A8)
δρ 1+4Eα ρ δρ
e F
andinthiswayrecoverourEqs.(9).Forthecaseofaradi-
APPENDIX A: GENERAL DISPERSION
ation dominated disc where α →0,
P
RELATION
γp 4λE/3ρ
Inthisappendix,weprovidethegeneraldispersionrelation Ce2ff = ρ + 1+α . (A9)
F
byincludingtheradiationenergyequationanddroppingthe
one-temperatureapproximation.Thismeansthatinthemo- WhenusingthisvalueforCeff,andsettingλ=1/3,Eq.(11)
mentumequation(1b),wehavetoreplace∇pby∇p+λ∇E, becomesidenticaltotheBlaes&Socrates(2001)dispersion
whereE denotestheradiationenergydensity.AsinEq.(7) relation.Thegeneraldispersionrelationprovidedinthisap-
λ denotes the flux limiter, which is constant in the linear pendix, thus, encompasses both the case of a gas-pressure
growthphaseduetotheconstantbackground.Theeffectof dominateddisc,suchasaprotoplanetarydisc,andtheoppo-
this is that in the linearised momentum equation (2b), we site extreme of a radiation-pressure dominated system like
should now define the effective sound speed C as an accretion disc around a black hole.
eff
δp δE
C2 ≡ +λ . (A1)
eff δρ δρ
REFERENCES
To determine it, we need the internal energy equation and
the radiation energy equation: Armitage P. J., 1998, ApJ, 501, L189+
Balbus S., Hawley J., 1991, ApJ, 376, 214
∂e
+∇·(ev)=−p∇·v−κ ρ(4πB−cE), (A2a) Balbus S. A., 2003, ARA&A, 41, 555
∂t P
Balbus S. A., Hawley J. F., 2006, ApJ, 652, 1020
∂E E
+∇·(Ev)=− ∇·v+κ ρ(4πB−cE)−∇·F. BalsaraD.S.,SpicerD.S.,1999,JournalofComputational
∂t 3 P
Physics, 149, 270
(A2b)
Blaes O., Hirose S., Krolik J. H., 2007, ApJ, 664, 1057
Here, B denotes the Planck function and κ the Planck Blaes O., Socrates A., 2001, ApJ, 553, 987
P
meanopacity.Thetermscontainingκ describethecoupling Blaes O. M., Balbus S. A., 1994, ApJ, 421, 163
P
ofthegastotheradiationviaemissionandabsorption.The Brandenburg A., 2008, Physica Scripta, p. 014016
radiative flux vector now reads Brandenburg A., Nordlund A., Stein R. F., Torkelsson U.,
1995, ApJ, 446, 741
λc
F =− ∇E, (A3) CourantR.,FriedrichsK.,LewyH.,1928,Math.Ann.,100,
κ ρ
F
32
where κF is the flux-mean total opacity. We define Fromang S., Nelson R. P., 2006, A&A, 457, 343
κ ρc k2 λc Fromang S., Papaloizou J., 2007, A&A, 476, 1113
αP ≡ Pσ and αF ≡ σ κ ρ, (A4) FromangS.,PapaloizouJ.,LesurG.,HeinemannT.,2007,
F A&A, 476, 1123
The dimensionless quantity αP gives the rate at which ra- Hawley J. F., 2000, ApJ, 528, 462
diation and matter are thermally coupled with respect to HawleyJ.F.,BalbusS.A.,WintersW.F.,1999,ApJ,518,
the growth rate of the instability while αF constitutes a 394
non-dimensionaldiffusioncoefficient.Weassumethatinthe HawleyJ.F.,GammieC.F.,BalbusS.A.,1995,ApJ,440,
equilibrium the gas temperature Tg and radiation temper- 742
ature Tr defined by E = aTr4 are equal: Tg = Tr ≡ T and Hirose S., Krolik J. H., Stone J. M., 2006, ApJ, 640, 901
4πB = aT4. Now the linearised equations corresponding to King A. R., Pringle J. E., Livio M., 2007, MNRAS, 376,
c
Eqs. (A2) become: 1740
„ « KissmannR.,2006,PhDthesis,Ruhr-Universita¨tBochum
p δT δE
δp=γρδρ−(γ−1)αPE 4 Tg − E , (A5) Kley W., Crida A., 2008, A&A, 487, L9
Krolik J. H., Hirose S., Blaes O., 2007, ApJ, 664, 1045
„ «
4E δT δE
δE = δρ+α E 4 g − −α δE; (A6) Kurganov A., Noelle S., Petrova G., 2001, SIAM J. Sci.
3ρ P T E F Comput., 23, 707
(cid:13)c ????RAS,MNRAS000,1–11