High-order schemes for non-ideal $3+1$ GRMHD: a study of the kinematic dynamo process in accretion tori

0.4 MB·

by  L. Del Zanna

**VolumeTitle**
ASPConferenceSeries,Vol.**VolumeNumber**
**Author**
(cid:13)c**CopyrightYear**AstronomicalSocietyofthePacific
High-orderschemesfornon-ideal3+1GRMHD:astudyofthe
kinematicdynamoprocessinaccretiontori
4 L.DelZanna123,M.Bugli14,N.Bucciantini23
1
0 1DipartimentodiFisicaeAstronomia,Universita` diFirenze,Italy
2
2INAF-OsservatorioAstrofisicodiArcetri,Firenze,Italy

n
a 3INFN-SezionediFirenze,Italy
J
4 4Max-Planck-Institutfu¨rAstrophysik,Garching,Germany
1

Abstract. We present the first astrophysical application of the ECHO code in its
]
E recent version supplemented by a generalized Ohm law, namely a kinematic study of
H dynamoeffectsinthickaccretiondisks.High-orderimplicit-explicitRunge-Kuttatime-
steppingroutinesareimplementedandvalidatedwithin3+1GeneralRelativisticMag-
h. netoHydroDynamics(GRMHD).Theschemeisappliedtoadifferentiallyrotatingtorus
p orbitingaKerrblackhole,wherethemean-fielddynamoprocessleadstostrongampli-
- fication of seed magnetic fields. We show that the interplay between the toroidal and
o
poloidal components occurs qualitatively in the same fashion as in the Sun, butterfly
r
t diagramsarereproduced, andatypicaltime-scaleforthefieldevolutionisfound, de-
s
pending on the dynamo and resistivity numbers, which could explain periodicities as
a
[ observedinseveralaccretingsystems.

1
v
3
2 1. Introduction
2
3
Large-scale, ordered magnetic fields are believed to play a crucial role to activate the
.
1 mechanisms responsible for the observed high-energy emission, from active galactic
0
nucleitogamma-raybursts. However,itisstillnotclearwhichisthepreciseoriginof
4
such fields. While the process of collapse to the compact objects which are believed
1
: to power these sources is certainly able to amplify any existing frozen-in field, the
v
question of its origin just shifts to the progenitors. Turbulent motions and small-scale
i
X instabilities,suchastheMagneto-RotationalInstability(MRI)inaccretiondisksorthe
r Taylerinstabilityinproto-neutronstars,areindeedgoodcandidatestoenhancealsothe
a
levelofturbulentmagneticfields,buttheresultingconfigurationswillbehighlytangled
andthusnotappropriatetogeneratetherequiredlarge-scalefields.
Aprocessnaturallyabletoamplifyorderedmagneticfieldsisthatofdynamo. Here
weshallconsidertheso-calledmean-field dynamoprocesses(Moffatt1978),i.e. when
small-scale turbulent (correlated) motions, invariably arising in astrophysical plasmas
withveryhighfluidandmagneticReynoldsnumbers,giverisetoaneffectivepondero-
motive force, along B, in the Ohm law. This is capable in turn to increase any initial
seedfieldwhencoupledtotheFaradayevolutionequation. Theresistive-dynamoOhm
lawcanbewritten,forclassicalMHD,as(herec → 1,4π → 1)
E(cid:48) ≡ E+(cid:51)×B = ηJ +ξB, (1)
1
2 L.DelZannaetal.
whereηisthecoefficientofresistivity(duetocollisionsandtothemean-fieldeffects),
andξ isthemean-fielddynamocoefficient(intheliteratureα = −ξ ismostcommonly
employed). Thelattertermispreciselyresponsibleforthedynamoprocess,leadingto
the generation of exponentially growing modes in the kinematical phase, that is when
the feedback of the increasing field on the other quantities can still be neglected. Dif-
ferentialrotationalsoplaysaroleinthedynamoprocess, andwhenbotheffectsareat
workwerefertoanαΩdynamo.
WhileinclassicalMHDtheaboverelationsimplydefinestheelectricfield,which
is a derived quantity of (cid:51) and B alone (since J = ∇ × B), in the relativistic case the
displacement current in the Maxwell equations cannot be neglected, and the evolution
equationfortheelectricfieldmustbesolvedaswell. Thisdifficultywasavoidedinpre-
vious relativistic studies of kinetic dynamo in accretion disks (Khanna & Camenzind
1996; Brandenburg 1996). To our knowledge the first generalization to full General
Relativity and its formulations for 3 + 1 numerical relativity can be found in Buc-
ciantini & Del Zanna (2013) (BDZ from now on), where also a systematic validation
of the second-order numerical scheme, within the ECHO (Del Zanna et al. 2007) and
X-ECHO (Bucciantini & Del Zanna 2011) codes, is presented. Here we briefly sum-
marizetheformalismandweimprovethenumericalmethodemployedtohigher-order.
Moreover,wepresent,asafirstnumericalapplication,astudyofthekinematicdynamo
actioninthickaccretiontoriaroundKerrblackholes. Moredetailsandadditionalsim-
ulationswillbepresentedinaforthcomingpaper(Buglietal.2014).
2. Basicequations
The fully covariant formulation for a resistive plasma with dynamo action, first pro-
posed in BDZ, is a relation for the quantities measured in the local frame comoving
withthefluid4-velocityuµ,namely
eµ = ηjµ+ξbµ, (2)
tobecomparedwithEq.(1),whereeµ,bµ,and jµ aretheelectricfield,magneticfield,
and current as measured by the observer moving with uµ. When written in the 3+1
formalism(Gourgoulhon2012),theMaxwellequationsbecome
(cid:16) (cid:17)
γ−1/2∂ γ1/2B +∇×(αE+β×B) = 0, (∇·B = 0), (3)
t
(cid:16) (cid:17)
γ−1/2∂ γ1/2E +∇×(−αB+β×E) = −(αJ −qβ), (∇·E = q), (4)
t
wheretheelectromagneticfieldsEandBarethosemeasuredbytheEulerianobserver,
αisthelapsefunction,βthespatialshiftvector,γ the3-metricwithdeterminantγ,q
ij
isthechargedensityand J theusualconductioncurrent. Inthe3+1language,Ohm’s
law(2)translatesinto
Γ[E+v×B−(E·v)v] = η(J −qv)+ξΓ[B−v×E−(B·v)v], (5)
herevistheEulerianvelocityandΓ = (1−v2)−1/2theusualLorentzfactor,fromwhich
onecanfinallycompute J,sothatthenewequationfor Eis
γ−1/2∂(γ1/2E)+∇×(−αB+β×E)+(αv−β)∇·E =
t
−αΓη−1{[E+v×B−(E·v)v]−ξ[B−v×E−(B·v)v]}, (6)
reducingto E = −v×Bforη = ξ = 0andtoEq.(1)inthenon-relativisticcase.
High-orderschemesfornon-ideal3+1GRMHD 3
3. Numericalsettingsandconvergencetests
IntypicalastrophysicalsituationstheresistivitycoefficientisverysmallandthusEq.(6)
is a stiff equation: terms ∝ η−1 can evolve on timescales τ (cid:28) τ , where the latter is
η h
the hyperbolic (fluid) timescale. We have then to solve, even for the simple kinematic
case,asystemofhyperbolicpartialdifferentialequationscombinedwithstiffrelaxation
equations,thatcanberewrittenintheform
X = γ1/2E : ∂ X = Q (X,Y)+R (X,Y), (7)
t X X
Y = γ1/2B : ∂Y = Q (X,Y), (8)
t Y
where Q are the non-stiff terms, whereas R are the stiff terms ∝ η−1 requiring
X,Y X
some form of implicit treatment in order to preserve numerical stability. To solve this
system, we here extend the numerical methods in BDZ (up to second order in time
and space) to higher orders. For spatial integration we use the full set of high-order
shock-capturingmethodsimplementedintheECHOcode,whereasfortimeintegration
we use the IMplicit-EXplicit (IMEX) Runge-Kutta methods developed by Pareschi &
Russo(2005)andfirstintroducedtoresistive(special)relativisticMHDbyPalenzuela
etal.(2009).
InFig.1weshowacoupleof1DnumericaltestsinMinkowskispacetime,already
describedinBDZ.Thefirstoneistheresistivediffusionofasmoothcurrentsheet,fol-
lowedfromt = 1tot = 10. Hereη = 0.01,ξ = 0,and p (cid:29) B2 toreducecompressible
effects(thistestismadebycouplingtheMaxwellequationstothefullsetofGRMHD,
wedonotconsiderthekinematiccase). Inthislimitthesystemremainsnon-relativistic
andtheanalyticalsolutionis
(cid:32) (cid:33)
x
By(x,t) = B erf √ , (9)
0
2 ηt
likefortheclassicaldiffusionequation. Inthetoprow,leftpanel,weshowtheanalytical
solution at the initial and final times, and the numerical one at t = 10, in the case
B = 0.01. In the right panel we show the L errors on By (comparing with a run
0 1
at extremely high resolution) as a function of the grid points number N. When we
usethethird-orderSSP3(4,3,3)IMEXmethod, anoverallsecondorthirdorderspace-
time accuracy is correctly achieved, depending on the spatial reconstruction method
employed. The best performing scheme is REC-MPE5 combined to DER-E6 (Del
Zannaetal.2007).
The second test concerns exponentially growing dynamo modes in a static, mag-
netizedbackground. Theanalyticalsolution,foragivenwavevectork,is
(cid:112)
By(x,t) = B exp(γt)cos(kx), γ = (2η)−1[ 1+4ηk(ξ−ηk)−1], (10)
0
whereγ isthedynamogrowthfactor. InFig. 1weshowinthebottomrow, leftpanel,
thetheoreticalandmeasuredamplificationof By(t)forvariouscombinationsofη(dia-
monds for η = 0.05, triangles for η = 0.1, squares for η = 0.25), ξ = 0.5, and k. The
correspondence is always matched. The case of η = 0.1 and k = 5 should not have
growing modes, due to the stabilizing effect of resistivity, but at late times (t (cid:39) 60)
small-scale fluctuations due to truncation errors triggers the fastest possible growing
modes anyway. In the right panel we show the L errors on By for various choices of
1
theIMEXschemes(andforREC-MPE5,DER-E6). Thenominalsecondorthirdorder
isachieved.
4 L.DelZannaetal.
Figure 1. Top panels: convergence study for the 1D current sheet test. Bottom
panels: convergencestudyforthe1Dsteadydynamotest.
4. Kinematicdynamoinaccretiontori
We choose as the background equilibrium the differentially rotating thick disk in Kerr
metricandBoyer-LindquistcoordinatesdescribedinDelZannaetal.(2007). Weselect
the maximally rotating case with a = 0.99, resulting in an orbital period of P = 76.5
c
(incodeunits)atthecenter. Thedomainisr ∈ [2.5,25](logarithmicallystretched)and
θ ∈ [π/4 − 0.2,3π/4 + 0.2], with a numerical resolution as low as 120 × 120 due to
thelong-termrunsneededtocapturethefulldynamics. TheIMEX-RKisSSP3(4,3,3),
thusthird-orderintime,andweuseREC-MPE5andDER-E6forspatialreconstruction
andderivation(thusfifthorderinspaceforsmoothprofiles).
In the kinematic case the dynamics is basically determined by the two dynamo
numbers, corresponding to the respective importance of the α and Ω effects with re-
spect to diffusion. We take Cξ = ξR/η ≥ 1 and CΩ = ∆ΩR2/η (cid:29) 1, where R is the
radial extension of the thick disk (and vertical too), ∆Ω the difference in the angular
velocitybetweenthediskcenterandtheinnerradius. Thechoiceofthesetwonumbers
determinethevaluesofξ andηatthediskcenter. Thecoefficientsarethenmodulated
within the disk as the square root of the mass density and proportional to it, respec-
tively,whilewetakeξ = 0andη = 10−5 inthelow-densityatmosphere(co-rotatingto
minimize shear flow numerical dissipation). The dynamo coefficient is further multi-
pliedbycosθbecausewewanttoobtainacombinationoftheαΩdynamoeffectwhich
issymmetricalwithrespecttotheequator.
The initial seed field (B ∼ 10−5) can be taken either as purely toroidal or purely
poloidal,inthelattercaseitisderivedfromapotentialA proportionaltothesquareof
φ
thelocalpressure. InFig.2weshowasimulationwithξ = η = 10−3 (maximumvalues
High-orderschemesfornon-ideal3+1GRMHD 5
Figure2. Dynamoevolutionofthetoroidalfieldforξ =η=10−3.
atthediskcenter),correspondingtoCξ = 5andCΩ = 400,withaninitialtoroidalfield.
This is shown as a function of time (in units of P ) and we can clearly see the period-
c
ical formation and migration of structures from the center to higher latitudes, always
confined within the disk. In Fig. 3 we show the maximum value of the toroidal and
poloidal orthonormal components of the field in the disk and their ratio, as a function
oftime. Weobservetheexponentialincreaseofbothcomponentsandthesaturationof
theratioaroundavalue B /B (cid:39) 0.13,whiletheoscillatingbehaviorissimplydueto
P T
adifferenceinphasebetweenthedynamomodesforthetwocomponents.
The periodical generation and migration of structures resembles very much what
happens in the Sun, where the αΩ is believed to be responsible for the 11-years sun-
spotscycle. Thisisusuallytrackedintheso-calledbutterflydiagrams,withperiodical
Magnetic Field Exponential Growth Ratio between B and B
105 0.20 P T
104 BBTP
0.15
103 B)T
max(B)(t) 110012 max(B)/max(P0.10
0.05
100
10-1 0.00
0 5 10 15 20 25 0 5 10 15 20 25
Time [PC] time [PC]
Figure3. GrowthofthetoroidalandpoloidalcomponentsandB /B ratio.
P T
6 L.DelZannaetal.
Model 1 Model 10
4 4
2 2
s [r]+ 0 s [r]+ 0
-2 -2
-4 -4
0 5 10 15 20 25 0 5 10 15 20 25
Time [P] Time [P]
C C

-0.5 0.0 0.5 -0.5 0.0 0.5
Figure4. Butterflydiagramsforξ =±10−3andη=10−3.
field reversaland migrationtowardsthesolar equator. InFig. 4we reportthe position
of the peaks along an appropriate new coordinate s as a function of time. In the left
paneltherunpreviouslydescribedisshown,withmigrationfromtheequatortohigher
latitudes (ξ > 0 in the northern hemisphere). In the right panel the sign of ξ has been
reversed(andthefieldinitializedtoapurelypoloidalone),tobetterreproducethesolar
situation(werecallthatα = −ξ insolardynamoapplications).
To conclude, in this paper we have shown the first implementation of high-order
IMEX-RKmethodsfortheresistive-dynamonon-idealGRMHDversionoftheECHO
code. PreliminaryrunsofthekinematicαΩdynamoeffectinthickaccretiontoriaround
maximally rotating Kerr black holes have ben performed as an astrophysical test case.
Similarlytothesolarcycle,wehavefoundperiodicalgenerationandmigrationofmag-
neticstructures,withexponentialgrowthofthemaximumvalueofthefield,expectedto
bequenchedbynon-linearfeedbackwhenthefullsetofGRMHDwillbesolved. The
(half) cycle of the dynamo process is seen to be around 6−7 orbital periods, but this
valuecanchangeaccordingtotheparameterschosen. Inanycasetheturbulenceprop-
ertiescansetanadditionaltimescale, throughthemean-fielddynamoeffect, unrelated
tothelarger-scalequantities. Thispropertymayhelptoexplainperiodicalmodulations
intheaccretingprocessandinthusthesourceluminosityitself(e.g.Gilfanov2010).
References
Brandenburg,A.1996,ApJ,465,L115
Bucciantini,N.,&DelZanna,L.2011,A&A,528,A101.1010.3532
—2013,MNRAS,428,71.1205.2951
Bugli,M.,DelZanna,L.,&Bucciantini,N.2014,MNRAS,inpress
DelZanna,L.,Zanotti,O.,Bucciantini,N.,&Londrillo,P.2007,A&A,473,11.0704.3206
Gilfanov, M. 2010, in Lecture Notes in Physics, Berlin Springer Verlag, edited by T. Belloni,
vol.794ofLectureNotesinPhysics,BerlinSpringerVerlag,17.0909.2567
Gourgoulhon, E. 2012, 3+1 Formalism in General Relativity, vol. 846 of Lecture Notes in
Physics,BerlinSpringerVerlag(Springer)
Khanna,R.,&Camenzind,M.1996,A&A,307,665
Moffatt, H. K. 1978, Magnetic field generation in electrically conducting fluids (Cambridge
UniversityPress)
Palenzuela,C.,Lehner,L.,Reula,O.,&Rezzolla,L.2009,MNRAS,394,1727.0810.1838
Pareschi,L.,&Russo,G.2005,JournalofScientificComputing,25,129