Table Of Content

Astronomy&Astrophysicsmanuscriptno.x-echo˙rev2 c ESO2011 (cid:13) January27,2011 GRMHD in axisymmetric dynamical spacetimes: the X-ECHO code N.Bucciantini1 andL.DelZanna2 1 NORDITA,AlbaNovaUniversityCenter,Roslagstullsbacken23,10691Stockholm,Sweden e-mail:[email protected] 1 2 DipartimentodiFisicaeAstronomia,Universita`diFirenze,LargoE.Fermi2,50125Firenze,Italy 1 e-mail:[email protected] 0 2 Received...;accepted... n ABSTRACT a J Wepresent anew numerical code, X-ECHO,for general relativisticmagnetohydrodynamics (GRMHD) in dynamical spacetimes. 6 Thisisaimedatstudyingastrophysicalsituationswherestronggravityandmagneticfieldsarebothsupposedtoplayanimportant 2 role,suchasfortheevolutionofmagnetizedneutronstarsorforthegravitationalcollapseofthemagnetizedrotatingcoresofmassive stars,whichistheastrophysical scenario believedtoeventually leadto(long) GRBevents. Thecode isbased ontheextension of M] theEulerianconservativehigh-order(ECHO)scheme[DelZannaetal.,A&A473,11(2007)]forGRMHD,herecoupledtoanovel solverfortheEinsteinequationsintheextendedconformallyflatcondition(XCFC).Wesolvetheequationsinthe3+1formalism, I assuming axisymmetry and adopting spherical coordinates for the conformal background metric. The GRMHD conservation laws h. aresolvedbymeansofshock-capturingmethodswithinafinite-differencediscretization,whereas,onthesamenumericalgrid,the p Einsteinellipticequationsaretreatedbyresortingtosphericalharmonicsdecompositionandsolved,foreachharmonic,byinverting - banddiagonalmatrices.Asasideproduct,webuildandmakeavailabletothecommunityacodetoproduceGRMHDaxisymmetric o equilibriaforpolytropic relativisticstarsinthepresenceof differentialrotationandapurelytoroidal magneticfield.Thisusesthe r sameXCFCmetricsolverofthemaincodeandhasbeennamedXNS.BothXNSandthefullX-ECHOcodesarevalidatedthrough t s severaltestsofastrophysicalinterest. a [ Key words.Gravitation–Relativisticprocesses–Magnetohydrodynamics (MHD)–Stars:neutron–Gamma-rayburst:general– Methods:numerical 2 v 2 1. Introduction ablyduetoefficientdynamoactionduringtheneutrinocooling 3 phaseinthehot,deleptonizingproto-NS(Duncan&Thompson 5 3 The most spectacular phenomenain high-energyAstrophysics, 1992).Ontheobservationalside,magnetarsaretheacceptedex- . like those associated to active galactic nuclei (AGNs), galactic planation for anomalousX-ray pulsars and soft gamma-rayre- 0 X-ray binary systems, or gamma-ray bursts (GRBs), typically peaters(Kouveliotouetal.1998). 1 involve the presence of rotating compact objects and magnetic On the computational side, the last decade has witnessed 0 fields. In some cases, such as the merging of binary systems a very rapid evolution in the construction of shock-capturing 1 : (formedbyeitherneutronstars,NSs,orblackholes,BHs)orthe codes for general relativistic MHD (GRMHD) in both static v collapseofrotatingcoresofmassivestarstowardsaNSorBH, and dynamical spacetimes, with a wealth of astrophysical ap- i X theinterplaybetweenmatter,electromagneticfields,andgravity plications to the situations outlined above (Font 2008). In the issostrongthattheMHDequationsgoverningthefluidmotions presentpaperwedescribeanovelcodeforGRMHDindynam- r a mustbesolvedself-consistentlywiththeEinsteinequationsfor icalspacetimes, namedX-ECHO, aimedat studyingthe evolu- the spacetime metric. Even for less violent phenomena, as the tionofmagnetizedrelativisticstarsandthegravitationalcollapse oscillations of neutron stars (Font et al. 2002), a self consis- of the magnetized rotating cores of massive stars. X-ECHO is tent solution of the fluid equationstogether with the spacetime built on top of the Eulerian conservative high-order code (Del evolution is essential to properly estimates the frequencies of Zannaetal.2007)forGRMHDinagivenandstationaryback- the eigenmodes.In the case of a binary NS merger,which is a groundmetric(Cowlingapproximation),whichinturnupgraded possible mechanism to account for short GRBs, a strong mag- thepreviousversionforaMinkowskianspacetime(DelZanna& netic field could be produced due to the induced shear (Price Bucciantini2002;DelZannaetal.2003).ECHOreliesonrobust &Rosswog2006).LongGRBsareinsteadassociatedtosuper- shock-capturingmethodswithinafinite-differencediscretization novaeventsandto the core collapseof massive stars(Woosley scheme (two-wave Riemann solvers and limited high order re- & Bloom 2006), leading to the subsequent formation of a ro- construction routines), with a staggered constrained transport tatingandstronglymagnetizedcompactobject.Themainstream method to preserve the divergence-freecondition for the mag- collapsarmodel(Woosley1993)impliestherapidformationof neticfield(tomachineaccuracyforsecondorderofspatialaccu- a maximally rotating Kerr BH at the center, accreting material racy),asproposedbyLondrillo&DelZanna(2000,2004).The from a torusand likely to loose energythroughthe Blandford- ECHO code has been already successfully applied to a variety Znajeckmechanism(Barkov& Komissarov2008).However,a ofastrophysicalsituationsinvolvingmagnetizedplasmasaround promising alternative for the GRB central engine involves the compactobjects,likethedynamicsandnon-thermalemissionof presence of a millisecond magnetar with B & 1015 G (Usov pulsar wind nebulae (Bucciantini et al. 2003; Del Zanna et al. 1992).Theoriginofsuchenormousfieldsformagnetarsisprob- 2004; Bucciantini et al. 2004, 2005a,b;Del Zanna et al. 2006; 1 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode Volpietal.2008),emissionofrelativisticMHDwindsfromro- recentlybytheextendedconformallyflatcondition(XCFC)for- tatingNSs(Bucciantinietal.2006),magnetarwindsproducing mulation(Cordero-Carriońetal.2009):alltheellipticequations longGRBjetsescapingthestellarprogenitor(Bucciantinietal. (noweightratherthanfive)arehierarchicallydecoupledandlo- 2008,2009),andpost-mergeraccretingdisksaroundKerrBHs caluniquenessisensured(seealsoSaijo2004),thusthiswillbe (Zanottietal.2010).Inspitethecodebeingfully3D,duetothe ourchoiceforX-ECHO. natureofthesources,invariablyaplasmasurroundingacentral In the present work we propose and test a new numerical compactobject,alltheaboveapplicationshavebeenperformed solverbasedonXCFCforanaxisymmetricspacetimeinconfor- in2Daxisymmetricspacetimesusingspherical-typecoordinates mally flat spherical-like coordinates. We employ the same nu- (either in Minkowski, Schwarzschild, or Kerr metric). The X- merical grid used for the evolution of the fluid and magnetic ECHOversionpresentedandtestedheresharesthesamephilos- quantitiesthroughtheECHOscheme,andthePoisson-likeequa- ophy,and, in view of the future applicationsmentionedabove, tionsaresolvedthroughahybridmethodbasedonsphericalhar- onlytheaxisymmetriccasewillbeconsidered. monicsdecompositionanddirectinversionofbanddiagonalma- TheEinsteinandGRMHDequationsinX-ECHOarewritten trices,resultingfromasecondorderfinitedifferencediscretiza- byfullyexploitingtheso-called3+1formalism(likeinECHO), tionoftheradialequations,foreachharmonic.Asaside prod- in which the original equations are split in their temporal and uct,webuildandmakeavailabletothecommunityanumerical spatialcomponents.The3+1formalismisnowadaysadoptedin code based on XCFC to produce self-consistent GRMHD ax- basicallyallnumericalschemesforgeneralrelativity(Alcubierre isymmetricequilibriaforpolytropicrelativisticstarsinthepres- 2008;Baumgarte&Shapiro2010),wherethesystemofEinstein ence of differential rotation and toroidal magnetic fields, here equationsistreatedlikeaCauchyproblemwithsomeinitialdata named XNS. Both XNS and the full X-ECHO codes are vali- to be evolved in time through hyperbolic equations. However, dated through several tests of astrophysical interest, including like for the solenoidal condition for the magnetic field, non- accuracy checks in the initial data for various NS equilibrium evolutionaryconstraintsmustbepreservedinthenumericalevo- configurations,accuracyinfindingthefrequenciesoftheirnor- lutionandcomputationalmethodsformoderncodesaredivided mal modesof oscillations, an evolutionarytest of migration of intotwomainclasses: 1)free-evolutionschemes,mainlybased NSunstableequilibriatostablebranches,atestonthestability onhyperbolicequationsalone,wherethisproblemisalleviated ofadifferentiallyrotating,magnetizedNSwithatoroidalfield, byappropriatereformulationsofthe equations(BSSN:Shibata 1D and 2D collapse of an unstable NS toward a BH, and a toy &Nakamura1995;Baumgarte&Shapiro1999),eventuallywith collapseofadifferentiallyrotatingNSwithpoloidalfields,asa theadditionofpropagatingmodesanddampingterms(Z4:Bona firststeptowardsmorerealisticmagneto-rotationalcorecollapse et al. 2003; Bernuzzi & Hilditch 2010); 2) fully constrained simulations. schemes, where the constraints are enforced at each timestep The paper is structured as follows. In Sect. 2 we introduce throughthesolutionofellipticequations(Bonazzolaetal.2004), and review the Einstein equations in the 3+1 formalism, first a morerobustbut computationallydemandingoption,since el- in their general form and then specialized in the CFC approx- lipticsolversarenotoriouslydifficulttoparallelize.Mostofthe imation, and we review the GRMHD equations. In Sect. 3 we state-of-the-art3DcodesforGRMHDindynamicalspacetimes discuss the new numerical XCFC solver assuming axisymme- are based on free-evolution schemes in Cartesian coordinates try and spherical coordinates for the conformal flat 3-metric, (Duez et al. 2005; Shibata & Sekiguchi 2005; Anderson et al. whereas the description of our novel XNS code for NS initial 2006;Giacomazzo&Rezzolla2007;Monteroetal.2008;Farris datacanbefoundinSect.4.Numericalvalidationandtestingof etal.2008),andhavebeenused,inthecaseofmagnetizedplas- variouscasesofNSequilibria,oscillationsandcollapsearepre- mas,forgravitationalcollapse(Duezetal.2006a;Shibataetal. sentedinSect.5,whileSect.6isdevotedtotheconclusions.In 2006a,b; Stephens et al. 2007, 2008), evolution of NSs (Duez the followingwe assume a signature ,+,+,+ for the space- {− } et al. 2006b; Kiuchi et al. 2008; Liebling et al. 2010), binary time metric and we use Greek letters µ,ν,λ,... (running from NSmergers(Andersonetal.2008;Liuetal.2008;Giacomazzo 0to3)for4Dspacetimetensorcomponents,whileLatinletters etal.2009,2010),andaccretingtoriaroundKerrBHs(Montero i, j,k,...(runningfrom1to3)willbeemployedfor3Dspatial etal.2010). tensorcomponents.Moreover,wesetc = G = M = 1andwe Provided the emission of gravitational waves is not of pri- absorbthe √4π factorsin the definitionof the ele⊙ctromagnetic mary interest, a good option in the class of fully constrained quantities. scheme is represented by the conformally flat condition (CFC) schemes(e.g.Wilson&Mathews2003;Isenberg2008),anap- proximation often employed for the study of gravitational col- 2. Basicequationsinthe3+1formalism lapse or NS stability and evolution (Dimmelmeier et al. 2002; In the present section we present the 3 + 1 formalism for Saijo2004;Dimmelmeieretal.2006;Cerda´-Durańetal.2008; Einsteinequations.Furtherdetailscanbefoundinrecentbooks Abdikamalov et al. 2009). CFC is typically associated to ax- and reviews on 3+1 numerical relativity (Gourgoulhon2007; isymmetric configurations and spherical coordinates, and it is Alcubierre 2008; Baumgarte & Shapiro 2010). We briefly dis- exact in the spherically symmetric case. Deviations from full cusstheconstrainedevolutionschemes,focussingontheelliptic GR solutions in the axisymmetric case for this kind of appli- CFC andXCFC solvers,andfinally,in Sect.2.3we reviewthe cations have already been shown to be negligible (Shibata & GRMHDequationsin3+1conservativeform,asimplemented Sekiguchi2004;Ottetal.2007),thoughthenonlinearequations intheoriginalECHOscheme(DelZannaetal.2007). mayshowseriousuniquenessproblemsforhighlycompactNSs ornascentBHs(seeCordero-Carriońetal.2009,andreferences therein). Moreover, CFC requires the solution of all the ellip- 2.1.TheEinsteinequationsforconformalflatness tic equationsfor the metric terms at the same time, usually by meansofiterationofPoissonsolvers,togetherwiththeinversion Thefieldequationsofgeneralrelativityexpressedinthe4Dfully ofconservativetoprimitivefluid/MHDvariables,anotheritera- covariantform tivenumericalprocess.Allthesedifficultieshavebeenresolved G =8πT , (1) µν µν 2 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode where G is the Einstein tensor containing the derivatives of wherethemainideaistomaximizethenumberofellipticequa- µν the metric tensor g and T is the matter and/or electromag- tions,usuallymorestablethanhyperbolicequations.Moreover, µν µν netic energy-momentumtensor, are not appropriatefor numer- in the steady state case the set of equations should easily re- ical computationssince time and space are treated on an equal ducetothoseusedforstationaryspacetimesandfortheconstruc- footing,whereasone wouldlike to cast themin theformof an tionofinitialdata.Anexampleistheso-calledfullyconstrained initialvalue(orCauchy)problemandevolvethemintime.The formalism (FCF) for asymptotically flat spacetimes in full GR mostwidelyusedapproachtoreachthisgoalisbasedontheso- (Bonazzola et al. 2004),which contains the widely used set of called3+1formalism,inwhichthegenericspacetime( ,g ) CFCellipticequationsasanapproximation.Inthefollowingwe µν issplitintospace-likehyper-surfacesΣ.IfnµindicatesthMetime- willdescribethegeneralassumptionsofconformalflatnessand t likeunitnormaltoΣ (alsoknownasthevelocityoftheEulerian wewillreviewthesetofCFCequations. t observer,n nµ = 1), the inducedthree-metric on each hyper- Let us start by applying a Lichnerowiczconformaldecom- µ − surfaceandtherelatedextrinsiccurvaturecanbedefinedrespec- position tivelyas γ :=ψ4f , (9) ij ij γ :=g +n n , (2) µν µν µ ν assumingaflatbackgroundmetric f (timeindependentandnot ij K := γλ n , (3) necessarilyinCartesian coordinates),wheretheconformalfac- µν − µ∇λ ν torsatisfies ψ = (γ/f)1/12, with f := detf . A secondassump- ij where µ isthecovariantderivativewithrespecttogµν (sothat tionistheconditionofmaximumslicingoffoliations ∇ g = 0). In general, any 4-vector or tensor can be decom- λ µν p∇osed into normal and spatial componentsby contracting with K =0. (10) nµorwiththeprojectorγµ,respectively.Inparticular,bothγ − ν µν Undertheseassumptions,tobepreservedduringtimeevolution, andK arepurelyspatial(andsymmetric)tensors. µν thetraceofEq.(5)anditstracelesspartbecome,respectively If xµ := (t,xi) are the spacetime coordinatesadaptedto the foliationof introducedabove,thelineelementisusuallywrit- ∂ lnγ1/2 ∂ lnψ1/6 = Dβi, (11) teninthesoM-calledADMform t ≡ t i 2αK = Dβ +D β 2(D βk)γ , (12) ds2 :=g dxµdxν = α2dt2+γ (dxi+βidt)(dxj+βjdt), (4) ij i j j i− 3 k ij µν − ij where∂tγ=0if,andonlyif,Diβi =0andtheextrinsiccurvature where the lapse function α and the shift vector βi (a purely canbeexpressedintermsofderivativesoftheshiftvectoralone. The next step is to use covariant derivatives associated to the spatial vector) are free gauge functions. In this adapted coor- dinate system the unit normal vector has components nµ = flat 3-metric fij, that will be indicated here by the usual nabla (1/α, βi/α) and n = ( α,0). In the 3 + 1 formalism, the operator i(and kfij =0).WhenEq.(9)holds,itispossibleto µ i ∇ ∇ − − demonstratethattheRicciscalarforγ (thatfor f iszero)is EinsteinequationsofEq.(1)aresplitintoasetofevolutionary ij ij equationsforγ andtheextrinsiccurvatureK ij ij R= 8ψ 5∆ψ, (13) − − ∂tγij = 2αKij+Diβj+Djβi, (5) in which ∆ := i is the usual Laplacian of flat space. The − i ∇∇ aboverelationcombinedfirstwiththeHamiltonianconstraintin Eq.(7)andthenwiththetraceofEq.(6)providesthefollowing ∂tKij =βkDkKij+KikDjβk+KjkDiβk DiDjα+ twoscalarPoisson-likeequationsfortheconformalfactorψand − α[R +KK 2K Kk]+4πα[γ (S E) 2S ], (6) thelapsefunctionα ij ij− ik j ij − − ij ∆ψ= [2πE+ 1K Kij]ψ5, (14) plusasetofconstraintsthatmustbesatisfiedatalltimes − 8 ij R+K2 KijKij =16πE, (7) ∆(αψ)=[2π(E+2S)+ 78KijKij]αψ5, (15) − wherewestillneedtowriteK intermsofflatspacederivatives ij D (Kij Kγij)=8πSi, (8) oftheshiftvector. j − TothetracelessextrinsiccurvatureA := K 1Kγ K named Hamiltonian and momentum constraints, respectively. ij ij− 3 ij ≡ ij HereD :=γµ isthecovariantderivativewithrespecttothe3- isthenappliedaconformaltime-evolutionrescaling i i∇µ metricγij (sothatDkγij = 0),Rij istheRiccitensor,againwith Kij =ψ−4A˜ij, 2αA˜ij :=(Lβ)ij, (16) respecttoγ ,R:=RithecorrespondingRicciscalar,K := Kiis ij i i wheretheconformalKillingoperatorassociatedtotheflatmetric thetraceoftheextrinsiccurvature.Asfarasthefluidsourcesare concernedE := n n Tµν,Si := n γiTµν,andSij := γi γjTµν andappliedtothevectorβi isdefinedas µ ν − µ ν µ ν (oftraceS := Si)are,respectively,theenergydensity,momen- (Lβ)ij := iβj+ jβi 2( βk)fij. (17) tum density, andi the stress-energy tensor as measured by the ∇ ∇ − 3 ∇k Eulerianobservers. This scaling is also employed in the so-called conformal thin Theconstraintsintroducedabovearenotoriouslydifficultto sandwich (CTS) approach to initial data. Moreover, it is quite bemaintainedinthenumericalevolutionofEqs.(5)and(6),and natural,duetoEq.(12),andsincewithinaconformallyflatde- twopossibleapproachescanbefollowed.Themostwidelyused compositionofthemetricwehave one relies on hyperbolicformulationsof the initial value prob- (L β)ij =ψ 4(Lβ)ij, (18) lem,theconstraintsareimposedonlyfortheinitialdataandnu- γ − mericalerrorsarejustmonitoredordampedduringtimeevolu- where L indicates here the conformalKilling operator associ- γ tion(freeevolutionschemes).Ontheotherhand,theconstraints atedtoγ .Ontheotherhand,inconformalflatnesswealsohave ij can be enforced at each timestep during the numerical simulation, leading to the so-called constrained evolution schemes, D Kij =ψ 10 (ψ10Kij), (19) j − j ∇ 3 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode tobeusedwiththeaboverescalinginthemomentumconstraint ∆ βi =16παψ 6fijSˆ +2Aîj (αψ 6), (30) L − j j − tofindanequationforβi. ∇ whereforconveniencewehaveintroducedrescaledfluidsource Thankstoalltherelationsderivedsofar,thefinalsetofCFC termsoftheform ellipticequationsmaybewrittenintermsofthesourcesandof A˜ij(containingαandfirstderivativesofβi)as Sˆ :=ψ6S , Eˆ :=ψ6E, Sˆ :=ψ6S, (31) j j ∆ψ=−[2πE+ 18fikfjlA˜ijA˜kl]ψ5, (20) andweremindthat ∆(αψ)=[2π(E+2S)+ 87fikfjlA˜ijA˜kl]αψ5, (21) Aîj =∇iWj+∇jWi− 32(∇kWk)fij. (32) ∆ βi =16παψ4Si+2ψ6A˜ij (αψ 6), (22) SomecommentsandcomparisonsbetweentheCFCandXCFC L j − ∇ setsofequationsarenowdue. where ∆ βi := (Lβ)ij =∆βi+ 1 i( βj), (23) – The unknown functions are now 8 rather than 5 (Wi, ψ, α, L ∇j 3∇ ∇j βi), and this is reflected by the augmentednumberof ellip- istheso-calledconformalvectorLaplacianoperator,associated ticequations(thereisanewvectorPoissonequationforthe totheflat3-metric fij andappliedtoβi. auxiliaryvariableWi). – WhileinCFCalltheequationswerestronglycoupled,here theequationscanbesolvedhierarchicallyonebyone,inthe 2.2.FromCFCtoXCFC givenorder,since eachrighthandside justcontainsknown A slightly different approach to the Einstein equations for functionsorthevariableitself(inthetwoscalarPoisson-like asymptotically flat spacetimes has been recently presented equationsforψandαψ). (Cordero-Carriońet al. 2009). This involves a rewriting of the – Aswe willsee inthenextsub-section,schemesforgeneral elliptical part of the FCF system for full GR through a differ- relativistic hydrodynamics or MHD (like ECHO), given a ent decompositionof the extrinsic curvature.Here we will just metric in 3+1 form, actually evolve in time the conserva- describeitsconformalflatnessapproximation,leadingtotheso- tivevariablesγ1/2S andγ1/2E,ratherthanSi and E.Since j calledextendedconformalflatnesscondition(XCFC)systemof ψ6 = γ1/2/f1/2 and f1/2 is known and time-independent, elliptic equations,improvingon the CFC onesdescribedin the the sources Sˆ and Eˆ are basically known after each com- j previous section. The set of XCFC equations is our choice for putational time-step without the need of an updated value themetricevolutioninX-ECHO. of ψ. This will be only needed to work out Sˆ = ψ6γ Sij, ij The new approach still relies on the usual conformal de- after the new value of ψ has been provided by Eq. (28) composition in Eq. (9) and on the maximum slicing condition and the inversion of conservativeto primitive variables has of Eq. (10),butthe choicefor the decompositionof the (trace- been achieved. Primitive variables are then updated self- less)extrinsiccurvatureisdifferent.Weuseherethemomentum- consistently together with the new values for the metric, constraintrescalingandtheso-calledYorkconformaltransverse whereasthiswasnotpossibleinCFC.Inthatcase,onecould traceless(CTT)decomposition,firstintroducedforinitialdata, either use Eq. (11) to derive a guess of the updated ψ (a thatis methodeasilypronetobothconvergenceproblemsanddis- Kij =ψ−10Aîj, Aîj :=(LW)ij+AîTjT, (24) cretizationerrors),oroneisforcedtoiteratesimultaneously overthemetricsolver(thewholeCFCset)andtheinversion wheretheconformalKillingoperatorassociatedtotheunknown routinefortheprimitivevariables(typicallyitselfanumeri- vector Wi gives the longitudinal part of Aîj , whereas Aîj TT TT caliterativeRaphson-Newtonmethod). is a transverse (∇jAîTjT = 0), traceless (fijAîTjT = 0) tensor. – Thelast,andcertainlynotleast,issueisrelatedtothemath- ConsistencybetweentheCTSandCTTdecompositions(notice ematicalnatureofthescalarPoisson-likeequations.Inboth thatAîj = ψ6A˜ij)shouldrequireanon-vanishingAîj .However caseswehaveastructureoftheform TT ithasbeendemonstratedthatthisquantityis evensmallerthan thenon-conformalpartofthespatialmetricwithintheCFCap- ∆u=hup, (33) proach, and hence can be safely neglected on the level of the where u is the generic variable (ψ or αψ), h is the generic CFCapproximation.Thusweset,asanadditionalhypothesis sourceterm,andpprovidestheexponentofthenon-linearity Aîj =0 Aîj =(LW)ij, (25) (p = 0fora canonicalPoissonequation).Itcanbe demon- TT ⇒ stratedthattheconditionph 0impliesthatthesolutionuis ≥ locallyunique.WhilethisisalwaystrueinXCFC,sincewe so that Aîj is definedin termsof the auxiliaryvectorWi alone. havetwocontributionswith p = 1and p = 7,bothwith Thelatteristobederivedfromthemomentumconstraintusing − − h 0, in Eq. (28), and one contribution with p = +1 and Eq.(19),thatissimply ≤ h 0inEq.(29),localuniquenesscannotbeguaranteedfor ≥ Aîj =∆ Wi =8πψ10Si, (26) theCFCsystem,sinceEq.(21)containsatermthatcertainly j L ∇ violatestherequirement(thesecondone,duetothepresence tobeaddedtotheotherCFCequations. ofafactorα 1inA˜ij). − The final augmented set of CFC elliptic equations, also knownasXCFCequations,isthenthefollowing 2.3.TheGRMHDequationsandtheECHOscheme ∆LWi =8πfijSˆj, (27) The equations for an ideal, magnetized, perfectly conducting plasmaare ∆ψ=−2πEˆψ−1− 18fikfjlAˆklAîjψ−7, (28) µ(ρuµ)=0, (34) ∇ ∆(αψ)=[2π(Eˆ +2Sˆ)ψ 2+ 7f f AˆklAîjψ 8]αψ, (29) Tµν =0, (35) − 8 ik jl − ∇µ 4 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode F µν =0, (36) retainingtheoriginalconservationform.Weendupwithasub- µ ∗ ∇ setoffluid-likebalancelawsindivergenceform respectively the continuity equation, the conservation law for momentum-energy,andthe sourcelessMaxwellequation.Here ∂U+∂Fi =S, (48) t i ρ is the mass density as measured in the frame comovingwith thefluid4-velocityuµ,thetotalmomentum-energytensoris plusamagneticsub-setwiththeinductionequationincurlform andtheassociateddivergence-freecondition,tobepreservedat Tµν =ρhuµuν+pgµν+FµλFνλ− 41(FλκFλκ)gµν, (37) alltimesduringevolution withh = 1+ε+ p/ρthespecificenthalpy,ǫ thespecificinter- ∂ i+[ijk]∂ =0, ∂ i =0. (49) t j k i nalenergy, p = p(ρ,ε)thethermalpressure(providedbysome B E B form of equation of state, EoS). Moreover, Fµν is the Faraday Thesetofconservativefluidvariablesandthesetofassociated (antisymmetric)electromagnetictensor,withtheassociateddual fluxesare,respectively F µν = 1ǫµνλκF , where ǫµνλκ = ( g) 1/2[µνλκ] is the space- ∗ 2 λκ − − D (αvi βi)D tgTpilmha=esemsLdyaees,vttwegim-µhCνiciiavhsnibtcdaleoc[psµoseνmedλueκbds]yoa-iOtsceohntnhmssoet’rrsaa(lliǫtanµewtνrλfnκofaort=iranagv−ap(Ln−eiersgvfh)eii1-cn/Ct2gl[iyvµeilνcteλaocκnts]rdy)iuc,mcwfitbieiontllhdg. U :=γ1/2SEj, Fi :=γ1/2ααSSiji−−−ββiiSEj, (50) whereasthesetofsourcetermscontainthederivativeofthemet- intheframecomovingwiththefluid ricandthusthecurvatureeffects Fµνuν =0. (38) 0 S:=γ1/2 1αSik∂ γ +S ∂ βi E∂ α . (51) TwhhiesIrnebotahsredice4ar-lcltyourdrreeepnriltvaJceνetsihsetahGdeeRMrMivaeHxdwDqeuellaqnuetaqittuyioaatnissoininnc∇3lµa+Fssµi1νcafl=oMrm−H,JDaνs,.  2 αKj iijkSij−iSjj∂j−α j  HereK istheextrinsiccurvatureintroducedintheprevioussec- ij employedinECHO,wemustdecomposeall4-vectorsandten- tion,whoseevolutionisdirectlyprovidedbytheEinsteinequa- sorsintotheirspatialandtemporalcomponentsoneachsliceΣ t tionsinthe3+1formalism,togetherwithγ ,ormaybegiven ij of the time evolution.Thiscan be easily achievedby using the intermsofthederivativesofthemetricterms.Foradynamical unitnormalvectornµintroducedintheprevioussection spacetime under the maximal slicing condition, from Eq. (12) wecanwrite ρuµ := D(vµ+nµ), (39) αK Sij = 1βjSik∂ γ +Sj∂ βi 1[γ 1/2∂(γ1/2βi)]S, (52) Tµν :=Sµν+nµSν+Sµnν+Enµnν, (40) ij 2 j ik i j − 3 − i Fµν :=nµEν Eµnν+ǫµνλκB n , (41) and the same expressionwithoutthe last term may be used for λ κ − anystationaryspacetime(Cowlingapproximation).Asfarasthe F∗µν :=nµBν Bµnν ǫµνλκEλnκ, (42) inductionequationisconcerned,wehavedefined − − whereeverynewquantityisnowpurelyspatial,asmeasuredby i:=γ1/2Bi, (53) theEulerianobserver.Inparticular,D:=ρΓistherestmassden- B fsaitcyt,ovri,iwshthoesefludiedfivneitlioocnitiys,dΓu:e=to(1t−hevicvoi)n−d1/i2tioisnthueuuµsu=alL1o.reTnhtez Ei:=αEi+ǫijkβjBk=−γ1/2[ijk](αvj−βj)Bk, (54) µ stress-energy 3-tensor Sij, its trace S, the momentum−density wherethevectorαvi βi issometimescalledtransportvelocity. − Si, and the energydensity E are the same quantitiesappearing The induction equation may be also written in the equivalent divergenceformas inthe3+1Einsteinequations(thisiswhynotationshavebeen slightlymodifiedwithrespecttotheoriginalECHOpaper),and ∂ j+∂ γ1/2[(αvi βi)Bj Bi(αvj βj)] =0, (55) arerespectivelygivenby tB i{ − − − } andtheantisymmetricnatureofthemagneticfluxesreflectsthat S =ρhΓ2vv +pγ EE BB +1(E Ek+B Bk)γ , (43) ij i j ij− i j− i j 2 k k ij of the originalelectromagnetictensor. Whateveris the adopted choicefortheformoftheinductionequation,wehaveafinalset S =ρh(Γ2 1)+3p+ 1(EEi+BBi), (44) − 2 i i of eight hyperbolic equations. Usually the corresponding vari- S =ρhΓ2v +ǫ EjBk, (45) ablesarenamedprimitivevariables,forexampletheset i i ijk E =ρhΓ2 p+ 1(EEi+BBi). (46) P:=[ρ,vi,p,Bi]T, (56) − 2 i i The electric and magnetic fields as measured by the Eulerian forwhichU = U(P),Fi = Fi(P),andS = S(P),wherefor observerare defined as Ei := n γiFµν and Bi := n γiF µν, simplicity we consider here the augmented system with j in andduetotheconditioninEq.(−38µ),tνheelectricfieldi−saµdeνriv∗ed U. B quantitypreciselyasinclassicalMHD As discussed in the previous section, the inversion of the non-linear system U = U(P) to recover the set of primitive E = ǫ vjBk, (47) variablesis achieved througha numericaliterative scheme and i ijk − requirestheknowledgeofthe volumeelementγ1/2, forconsis- where ǫ = γ1/2[ijk] (ǫijk = γ 1/2[ijk]) is the Levi-Civita tency updated at the same time level as U. Moreover, as dis- ijk − pseudo-tensor for the 3-metric γ , and [ijk] is the alternating cussedinDelZannaetal. (2007),theconservativetoprimitive ij symboltakingvalues+1, 1,or0. variablesinversionisthemostdelicatepartofarelativisticMHD − Thankstotheabovedecompositions,theGRMHDequations code,ashighLorentzfactorflowsorstrongmagneticfields(i.e. canbeentirelyrewrittenintermsofpurelyspatialvectors,while in a NSmagnetosphere)mayeasily lead toerrorsin thevalues 5 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode of the conservative variables. The whole procedure can be re- while the relation to covariant components still involves the ducedtothesolutionoftwocouplednonlinearequations,forany function ψ, since X = ψ4f Xj. As far as covariant derivatives i ij givenEoS.InX-ECHOweleavecompletefreedominthechoice areconcerned,thechangeofbasisallowsonetousethe op- ∇î p= p(ρ,ε),howeverinthepresentpaper,giventhenatureofthe eratorofsphericalcoordinates.Inparticular,theLaplacianofa numericaltestsproposed,welimitthischoicetoeitheranideal genericscalarfunctionu(r,θ)is γ-lawEoS ∆u=∂2u+2r 1∂ u+r 2(∂2u+cotθ∂ u), (65) p(ρ,ε)=(γ 1)ρε, (57) r − r − θ θ − ortoapolytropiclawEoS whereas the orthonormal components of the conformal vector Laplacianare,respectively 1 p(ρ)= Kργ, ,ε(ρ)= Kργ−1, (58) (∆ X)rˆ =∆Xrˆ 2r 2(Xrˆ+∂ Xθˆ+cotθXθˆ)+ 1∂ ( X), (66) γ 1 L − − θ 3 r ∇· − wthheerroeoγt-afinnddKingarpergoicveednucreonssotlavnetss.fWorhtehnetvhaerfiiarbstleopxti=onvisvui,saecd-, (∆LX)θˆ =∆Xθˆ+2r−2∂θXrˆ−(rsinθ)−2Xθˆ+31r−1∂θ(∇·X), (67) cordingtothethirdmethoddescribedintheECHOpaper,iwhere (∆LX)φˆ =∆Xφˆ (rsinθ)−2Xφˆ, (68) − thenestedsecondequation(acubic)fory = ρhΓ2 canbeeither wherethedivergenceofXis solvedanalyticallyoriterativelywithaninnerloop(andwetyp- icallyadoptthislatterchoice).NoticethatwhenSiBi = 0,orin X=∂rXrˆ+r−1(2Xrˆ+∂θXθˆ+cotθXθˆ), (69) generalforapurelyfluidsimulation,theequationforyislinear ∇· and can be readily solved. On the other hand, for a polytropic precisely the formulae of vector calculus in spherical coordi- lawtheenergyequationbecomesredundant,sinceallthermody- nates. namical quantities are now functions of the density alone, and The Poisson-like elliptic equations used in XCFC to com- theoverallinversionproceduresimplyreducestoaroot-finding putetheeightmetrictermsconsistoftwoofscalarequationsin iterativesolverforx. the formof Eq. (33), for the variablesu = ψ and u = αψ, and two vector equations for the generic unknown vector Xi = Wi andXi =βi.Duetonon-linearity,thescalarequationsarebetter 3. TheXCFCsolverforaxisymmetricspacetimes solved iteratively for the quantity q := u 1, which in both n n − cases gives the deviation from asymptotic flatness ψ 1 and The X-ECHO code for GRMHD in dynamical spacetimes is → α 1,wherethenewvalueatstepniscomputedusingthepre- builtupontheECHOscheme,coupledtoanovelsolverforthe → vious value at step n 1 in the source term, until convergence XCFCequationsdescribedinSect.2.2.Allthenecessarydefini- − is reachedwithin some prescribedtolerance.Summarizing,the tionsandequationshavebeenprovidedintheprevioussection, metricequationsareexpressedinoneofthetwogenericforms herewespecializetotheparticularimplementationofX-ECHO wearemostlyinterestedin,thatisundertheassumptionofax- ∆q = H h(1+q )p, (70) n n 1 n 1 isymmetricGRMHDconfigurations,adoptingspherical-likeco- − ≡ − ordinatesxi =(r,θ,φ)fortheconformalflatmetric. and,usingtheorthonormalbasisintroducedabove, Thus,asafirststepweassumeinEq.(9)theusualspherical (∆ X)î = Hî, (71) coordinates L f =diag(1,r2,r2sin2θ), (59) ij where h and H are generic scalar and vector source terms, to sothat beprovidedbytheXCFCequations,sothatboththerighthand sidesaboveareknownfunctionsof(r,θ). ds2 = α2dt2+ψ4[(dr+βrdt)2+r2(dθ+βθdt)2+r2sin2θ(dφ+βφdt)2], Numerical methods available for solving these elliptic par- − (60) tialdifferentialequations(PDEs)canbedividedintothreemain with f1/2 = r2sinθ and γ1/2 = ψ6r2sinθ. Moreover, we are categories(see,formethodsanddiscussionsGrandcle´mentetal. goingtospecializetotheaxisymmetriccase,thusthecondition 2001;Dimmelmeieretal.2005;Grandcle´ment&Novak2009): directinversion,fullrelaxation,spectral.Directinversioncodes ∂φ 0, (61) areabletosolvethecompletesystemofCFC(orXCFC)atonce ≡ usingNewton-Raphsonsolver(oranyotherinversiontechnique) will be assumed throughout the paper. Within the flat metric ontheentirecomputationalgridoverwhichthemetricsolution fij, it is convenient to introduce the orthonormal basis eî := is desired. Theyhave verygoodconvergenceto machine accu- (erˆ,eθˆ,eφˆ),with racywithinfewsteps,buttheysufferseriouslimitations:theini- tial guess must be close enough to the solution to avoid con- erˆ :=∂r, eθˆ :=r−1∂θ, eφˆ :=(rsinθ)−1∂φ, (62) vergenceon localminima (instead of globalminima);memory requirementformatrixallocationistypicallyverylarge(usually for which f = diag(1,1,1), where a similar notation as in îˆj a sparse matrix is needed in the whole 2D domain);in general Bonazzola et al. (2004)has been assumed. Any generic vector direct inversionschemes solve the metric on smaller grids that Xcanbethenexpressedintheusualformas the one over which the fluid variablesare evolved, and require interpolationbetween the two, with problemsthatmay arise at X := Xrêrˆ+Xθêθˆ+Xφêφˆ, (63) boundaries. Full relaxation codes use SOR, multigrid or other relaxation techniques. The schemes are fast, they require little wheretheorthonormal(withrespectto f ) componentsXî are, memory allocation, but usually suffer from poor convergence ij respectively properties:theconvergenceisin generalslow(comparedtodi- rectorspectralschemes)andmightfailontheaxisoratthecen- Xrˆ := Xr, Xθˆ :=rXθ, Xφˆ :=rsinθXφ, (64) ter,duetothesingularnatureofsomemetricelements.Spectral 6 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode schemesdecomposethesetofCFC(orXCFC)equationsusinga forC(r)alone.Thenewsourcetermsaregivenby l combinationofsphericalharmonics(basedonLegendrepolyno- mials) in the angulardirectionsand Chebyshev polynomialsin Hrˆ(r):= Hrˆ(r,θ)Y(θ)dΩ, (80) theradialdirection.Thisensuresacorrectbehaviorontheaxis l I l andatthecenterevenwithalimitednumberofeigenfunctions, 1 buttheyrequirespecificgridsandsometimescomplexcompacti- Hlθˆ(r):= l(l+1)I Hθˆ(r,θ)Yl′(θ)dΩ, (81) ficationsormulti-domaindecompositiontechniqueswithappro- priateboundaryconditionsforeachmultipoleandeachdomain. Hφˆ(r):= 1 Hφˆ(r,θ)Y (θ)dΩ. (82) The metric solver in X-ECHO uses a mixed technique. In the l l(l+1)I l′ angular direction we decompose using spherical harmonics, to These integrals, as well as that in Eq. (75), are computed in preservethe correctasymptotic form on axis. However,the set X-ECHO by using Gaussian quadrature points, so the original ofordinarydifferentialequations(ODEs)obtainedforeachhar- sourcetermsmustbefirstinterpolatedontherequiredlocations. monicisthensolvedusingdirectinversionoverthesameradial Asanticipated,weusefinitedifferencesandasecondorder gridusedintheGRMHDcode,withnoneedofinterpolationor approximationforfirstandsecondspatialderivatives,sothat,for compactification.Atsecondorderaccuracyinafinitedifference eachharmonicl,thetwocoupledpoloidalequationsarereduced discretization,thescalarequationsreducetothesimpleinversion totheinversionofasparsematrixofbandwidth7,whereasthe oftridiagonalmatrices(banddiagonalmatricesforthepoloidal toroidal equation leads to a tridiagonal matrix (standard open componentsofthevectorPoissonequations),whereappropriate source routines are employed). Boundary conditions at r = 0 solvers are fast, require little memory allocation, and typically andr=r aregivenbytheparityandasymptoticpropertiesof convergewithhighaccuracy. max themultipolecorrespondingtotheharmonicl.Inparticular,we ForthescalarPoisson-likeEq.(70)wethendecompose,for assumethatatthecenterA(r)hasparity( 1)l,whereasB(r)and eachlevelnofiteration(hereomitted),theunknownqas C(r)haveparity( 1)l+1,alndatlargedis−tancesfromthelcentral l sourcesthemultip−olesareforcedtodecayasr l(l+1). ∞ − q(r,θ):= [A(r)Y(θ)], (72) l l Xl=0 4. AxisymmetricGRMHDequilibria:theXNScode where One of the most common application of the fully constrained Y(θ) Y0(θ):= 2l+1P(cosθ), (73) l ≡ l q 4π l formalism is the search of self-consistent(axisymmetric)equi- libriumconfigurations(i.e.fluidquantitiesandmetric)forcom- with P the Legendre polynomial of degree l and the axisym- l pactrelativisticstars(wewillsimplyrefertotheseobjectswith metry condition has been imposed (m = 0). The PDE, where thetermNS)atypicalcaseofinitialdataprobleminGR(Cook theLaplacianisprovidedinEq.(65),maybethensplitintothe 2000;Stergioulas2003;Gourgoulhon2010).Severalcodeshave seriesofradialODEsforeachharmonicl been presented in the years addressing this issue (Komatsu d2A 2dA l(l+1) etal.1989a,b;Cooketal.1994;Stergioulas&Friedman1995; l + l A = H, (74) dr2 r dr − r2 l l Nozawa et al. 1998; Bonazzola et al. 1998; Kiuchi & Yoshida 2008),andallofthem,despitethedifferentapproachesandup- wherethenewsourcetermis grades(e.g.differentialrotation,toroidalmagneticfield),adopt theso-calledquasi-isotropiccoordinates.Undertheconditions H(r):= H(r,θ)Y(θ)dΩ, (75) l l I ∂ 0, ∂ 0, (83) t φ ≡ ≡ with dΩ = 2πsinθdθ and the integral running from θ = 0 to wewriteherethecorrespondinglineelementintheform θ=πduetoaxisymmetry. As far as the vector Poisson equation in Eq. (71) is con- ds2 = α2dt2+ψ4(dr2+r2dθ2)+R2(dφ ωdt)2, (84) − − cerned, the unknown vector X is first decomposed into vector whichresemblesthatoftheCFCmetricforaxisymmetricspace- sphericalharmonics,thatis,intheaxisymmetriccase timesofEq.(60)insphericalcoordinatesandreducestoitwhen X(r,θ):= ∞ A(r)Y(θ)e +B(r)Y (θ)e +C(r)Y (θ)e , R=ψ2rsinθ, (85) l l rˆ l l′ θˆ l l′ φˆ Xl=0 h i where only in this case the function ψ in Eq. (84) recover the (76) meaningofconformalfactor.Ingeneral,wecanthinkthemetric where Y := dY/dθ. As it is apparent from the operators in Eqs. (66l′-69), thelset of equations split into a series of ODEs functionR := γ1/2 tobeasortofgeneralizedcylindricalradius, φφ with,foreachharmonicl,acoupledpoloidalpartfortheradial whereasω:= βφistheintrinsicangularvelocityaboutthesym- − functionsA(r)andB(r) metry axis of the zero angularmomentum observers (ZAMOs: l l Bardeen et al. 1972), that are the normal (Eulerian) observers 4d2Al + 8 rdAl A l(l+1) A + r dBl 7B =Hrˆ, for an axisymmetric spacetime. When ω = 0 the spacetime is 3 dr2 3r2 dr − l!− r2 l 3 dr − 3 l! l sphericallysymmetricandthetwometricsbothreducetothatin (77) isotropicSchwarzschildcoordinates.Inthefollowingwebriefly d2Bl + 2dBl 4l(l+1)B + 1 dAl + 8 A = Hθˆ, (78) re-derive the GRMHD equilibrium condition (a Bernoulli-like dr2 r dr − 3r2 l 3r dr 3r2 l l integral)forpurelytoroidalflowsandmagneticfieldsinquasi- isotropiccoordinates,moregeneralderivationscanbefoundin andatoroidalpart (Kiuchi&Yoshida2008)or(Komissarov2006),thislatterwork dd2rC2l + 2r ddCrl − l(lr+21)Cl =Hlφˆ, (79) anpupmlieerdictaolmteasgtninetDizeeldZtoanrinaaroeutnald.r2o0t0at7i)n.gBHs(seealsothefinal 7 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode Theonlynonvanishingcomponentsofthespatial(Eulerian) non-magnetized case, it has been demonstrated that, even for 3-vectorsviandBiaretheazimuthalones,andthecorresponding rapid rotators close to the mass shedding limit (see Sect. 5.3), moduliare solutions obtained with the RNS numerical code (Stergioulas & Friedman 1995) show very little deviations from a CFC v:=(v vφ)1/2 =Rvφ =α 1R(Ω ω), (86) φ − metric, so in principle one would expect the CFC limit to − provideareasonablygoodapproximationofthecorrectsolution B:=(B Bφ)1/2 =RBφ, (87) φ for compact relativistic stars. Given the above mentioned where we have used the 3 + 1 relations uφ = Γ(vφ + ω/α), computational advantages in the solution of the XCFC system ut = Γ/α, and we have defined Ω := uφ/ut = dφ/dt, the an- ofequationswith respectto the originalsetof coupledPDE of gularvelocityofthefluidseenbyanobserveratrestatinfinity fullyconstrainedschemes,ourchoicehasbeentobuildanovel (afunctionofrandθforadifferentiallyrotatingNS).Theequi- numerical code (written in Fortran90),which we name XNS libriumconditionwearelookingforcanbedirectlyderivedfrom andthatcanbefreelydownloadedat the3+1conservativeformoftheGRMHDequations,inwhich the only non vanishing contribution is due to the two poloidal http://sites.google.com/site/niccolobucciantini/xns componentsofthemomentumconservationinEq.(48) This takes advantage of the same XCFC metric solvers γ 1/2∂ [γ1/2α(p+ 1B2)] α(p+ 1B2)1γii∂ γ developed for the X-ECHO scheme described in the present − j 2 − 2 2 j ii section.Acomparisonbetweenequilibriafoundwith XNSand α(ρhΓ2v2 B2)1R 2∂ R2+ρhΓ2Rv∂ ω − − 2 − j j RNSispresentedinSect.5.3,togetherwithadiscussionofthe +(ρhΓ2 p+ 1B2)∂ α=0, (88) results. − 2 j XNSemploysaself-consistentmethodinthesearchforthe where jiseitherrorθandweremindthattheelectricfieldvan- axisymmetricequilibriumsolutionsofrelativisticcompactstars, ishes. Recalling that 12γii∂jγii = γ−1/2∂jγ1/2 and differentiating inthepresenceofdifferentialrotationandapurelytoroidalmag- thedefinitionofv,afterafewalgebraicstepsthefollowingcon- neticfield,thusmetrictermsandfluid-likequantitiesarederived ditioncanbefound: at the same time. Given the values of the six free parameters ∂ p Γ2v2∂ Ω ∂ (α2R2B2) K,n,A,Ωc,Km,m,plusaguessforthecentraldensityρc,thefol- j j j +∂ lnα ∂ lnΓ+ + =0. (89) lowingstepsaretaken: ρh j − j Ω ω 2α2R2ρh − Integrabilityofthisequationdemandsthefollowingconditions: – A starting guess for the CFC metric terms (α,ψ,ω) is pro- videdfromthepreviousstep,withRgivenbyEq.(85).The – A barotropicEoS, p = p(ρh). The simplest choice and the firsttime, the sphericallysymmetricTolman-Oppenheimer- mostcommonassumptionisapolytropiclaw Volkoff (TOV) solution in isotropic coordinates (Tolman 1934),forthe metric ofa non-rotatingandnon-magnetized p= Kρ1+1/n h=1+(n+1)Kρ1/n, (90) NS with a central density value ρ , is computed through a ⇒ c shootingmethodforODEs. wherenisthepolytropicindex(thecorrespondingadiabatic – Ontopofthesemetrictermsandforeachgridpoint,Ωisde- indexisγ=1+1/n). rivedbyinvertingEq.(91),thenv(andΓ)canbedetermined – ThequantityF :=u ut Γ2v2/(Ω ω),relatedtothespecific angularmomentumφℓ :=≡ u /u,is−afunctionoftheangular fromEq.(86).Finally,hc isaknownfunctionofρc andthe − φ t BernoulliintegralinEq.(93)issolvedviaaNewtonmethod velocityΩalone.Acommonlyadopteddifferentialrotation tofindthelocalvaluesofhandρ,allowingustodetermine law(e.g.Stergioulas2003)is alsothemagneticfieldstrengthfromEq.(92). R2(Ω ω) – Combiningtheupdatedfluidquantitieswith theoldmetric, F(Ω)= A2(Ωc−Ω)= α2 R2(Ω− ω)2, (91) thenewconservedquantitiesEˆ,Sˆjarederived. − − – The set of XCFC equations is solved for this set of con- whereΩ isthecentralangularvelocityand Aisameasure served variables, and a new metric is computed. Here only c ofthedifferentialrotationrate.ForuniformrotatorsΩ Ω the azimuthalcomponentsfor the vector Poisson equations c (and A ), and this contribution can be excluded≡from aretreated,inordertofindWφ andthenβφ.Weremindhere → ∞ theBernoulliintegral. thattheconservativetoprimitivevariablesinversionmustbe – A sortofmagneticbarotropiclaw,whereαRBisafunction enforcedbetweenthesolutionstothetwoscalarPoisson-like ofα2R2ρh.Thesimplestchoiceisamagneticpolytropiclaw equations,namelyafterthenewvalueofψisfoundandbe- forethatofα. αRB= K (α2R2ρh)m, (92) m Thesestepsarerepeateduntilconvergencetoadesiredtolerance where m is the magnetic polytropic index, with m 1 ≥ isachieved.However,a wordofcautionisin orderhere.Since (Kiuchi&Yoshida2008). XNSisbasedontheXCFCmetricsolver,thatworksonthecon- Using the above prescriptions we can easily derive the final servative variables densitized with ψ6, convergence is actually GRMHDBernoulliintegral enforcedonthecentralvalueofthequantityDˆ := ψ6D =ψ6ρc. Therefore, given that the final conformal factor ψ for the self- h α A2 mK2 consistent 2D equilibrium may be quite different from that de- ln +ln lnΓ (Ω Ω)2+ m (α2R2ρh)2m 1=0, h α − − 2 c− 2m 1 − rivedfromtheradialTOVsolutionatthefirststep,thefinalvalue c c − (93) ofρatthecenterisexpectedtodifferfromtheparameterρ inthe c whereagainwiththesubscript weindicatevaluesatthecenter. Bernoulliintegral.Ifonewantstofindanequilibriumconverg- c The above derivation is exact for quasi-isotropic coordi- ingexactlytoacentraldensityρ ,anadditionaloveralliterative c nates and obviously applies also to the CFC sub-case. In the loopisneeded. 8 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode Before concluding the section, some remarks are in order. have started to appear.This, togetherwith the intrinsic degrees One might question the choice of a purely toroidal field ver- offreedominthechoiceofgaugeandcoordinatesystems,typi- sus a more realistic configuration including a poloidal compo- calofGR,haveresultedinalackofwelldefined,agreedupon, nent.However,equilibriumwithpoloidalfieldsisonlypossible standardnumericaltests(inthe spiritofwhatshock-tubeprob- for uniform rotators, with magnetic field fully confined within lemsareinflatspacetimeoraccretion/outflowssolutionsarefor the star. It is well known (Goldreich & Julian 1969) that any a stationary curved metric). While in 1D a few problems have poloidalmagneticfieldextendingoutsidearotatingNSwilllead emergedasstandardbenchmarks,thiscannotbesaidaboutmul- eventuallytoaspin-downofthesame,evenforadipolealigned tidimensionalcasesyet,alsobecauseofalackofanalyticalso- with the rotation axis, unless the star is surrounded by a true lutionsforfullymultidimensionaldynamicalproblems.Someof vacuum. A small charge density of order of 1020cm 3 pairs is theteststhatwehaveselectedhere,havebeenalreadypublished − enough to break this assumption, even in the case of millisec- in the literature, using different numerical schemes, both fully ond rotators with B 1017 G, and the timescale to replenish constrained (Dimmelmeier et al. 2002; Stergioulas et al. 2004; ∼ an evacuated magnetosphere is of order of the rotation period. Dimmelmeieretal.2006;Cordero-Carriońetal.2009)andhy- For such strong magnetic fields, affecting the overall equilib- perbolic (Font et al. 2002; Bernuzzi & Hilditch 2010), so that rium of the NS and providing a non negligible contribution to we can validate our code against existing results. Some other the globalstress-energytensorin the Einstein equations(lower have been done in the perturbative regime and compared with fields have little dynamicaleffects and can be easily treated as theresultsoflineartheory,andwhenpossiblealsowithprevious perturbations)theproblembecomesparticularlyseverebecause numerical simulations. However, in an attempt to evaluate the the typical spin-down time for millisecond rotators can be as performancesunder different conditions,we also present some shortas50ms.Weakermagneticfieldcanleadtomuchlonger novelcases. spin-downtimes,andthoseconfigurationsmightbeconsidered quasi-equilibriumcases,atleastforthetimeofatypicalnumer- Unless otherwise stated, in all our numerical tests we will icalrun(afewthousandsM).Howevertheyareoflittleinterest, useaCourantnumberof0.4,anidealgasEoSwithγ = 2(cor- asstatedabove. responding to a polytropic index n = 1 for the initial data in One might also question if purely toroidal configurations XNS), and we will solve the full GRMHD system, including are stable, and, if not, what is the growth rate of instabilities. the equation for the total energy density E, often neglected in A recent study of the stability properties of neutron star with isentropictests fora givenpolytropiclaw.For theapproximate strongtoroidalfieldhasbeenpresentedbyKiuchietal.(2008). Riemannsolver,weusehereforthefirsttimetheHLLCsolver Their results show that non-rotating neutron stars with a mag- by (Mignone& Bodo 2006),never applied before to GRMHD netic polytropic law with m 2 are always unstable, and that ≥ studiesin curved,evolvingspacetimes, to ourknowledge.This (uniformly)rotatingsystems are stable only if their kinetic en- choice is imposed by the sharp transition between the rotating ergyexceedbyatleastafactor5themagneticenergy.However, NSandtheexternalatmosphere,sincethetwo-wavesolverHLL eveniftheirconclusionsaboutstabilityaresupportedbynumer- isfoundtobetoodiffusiveonthecontactdiscontinuity. icalsimulations,itmustbepointedoutthattheircriterionisde- rivedanalyticallyonlyintheNewtonianregime,wheremagnetic fieldandspacetimemetricarenotcoupled.Moreover,theirfull Given that the astrophysical problems we are mostly inter- GRresultsdonotconsiderintermediatecaseswith1 < m < 2, estedin,andtowardwhichtheX-ECHOschemehasbeendevel- thusfurtherinvestigationiscertainlyneeded. oped,involvetheabilityofsimultaneouslyhandlethehighden- sityNSandanylowdensityoutflow/atmosphere/magnetosphere that might surround it, in almost all of our tests we have in- 5. Numericalresults cludedinthedomainanextendedatmosphere,whichisleftfree toevolveandrespondtotheevolutionoftheNS.Weareaware WepresenthereasetofnumericaltestsoftheX-ECHOscheme. that in the literature a common practice is to reset floor values StandardHD/MHDtestsinastaticbackgroundmetric(Cowling outsidetheNS,butwebelievethatthisproceduremayinprin- approximation)havebeenalreadypresentedelsewherebothfor cipleleadtoviolationsofconservationpropertiesofthescheme a flat metric (Del Zanna & Bucciantini 2002; Del Zanna et al. attheNSsurface.Forthesamereasons,asstatedabove,simula- 2003) and for a given stationary curved spacetime (Del Zanna tionshavebeenperformedusingaidealgasEoS,whichallows etal.2007).Asdiscussedintheintroduction,theperformances us to handle systems where matter in different thermal condi- oftheECHOcodeinastrophysicalscenariosinvolvingavariety tions(acoldNSversusahotatmosphere)ispresent. of different conditions like strong shocks, relativistic outflows, strong gravity, and highly magnetized systems, have also been alreadyassessed innumerouspapers.Forthesereasons,there- Reconstructionatcellboundariesisachievedforsimplicity sults presented here focus mostly on evaluating the quality of throughamonotonized-central(MC)algorithm,thoughtheother ournovelmetricsolver,andtheperformancesoftheHD/MHD choicesdescribedintheappendixof(DelZannaetal.2007)are ECHO algorithms when coupled with a dynamical spacetime. alsopossible.TheXCFCmetricsolverisinvokedevery10steps Moreover, the original ECHO scheme was designed for high- oftheHD/MHDRunge-Kuttaevolutionscheme.In2Drunswe order accuracy, whereas the metric solver implemented in X- use50Gaussianquadraturepointsand20sphericalharmonics. ECHO is formally only second order in space. Given that the With these settings, we have managedto make comparablethe performances of high-order reconstruction techniques has al- computationaltimestakenbythemetricsolver,usuallyslower, ready been tested (Del Zanna et al. 2007), for simplicity, we andbyasingleRunge-Kuttacycleofthefluidsolver.Thus,solv- havedecidedtolimitourset-uptosecondorderaccuracy,both ingXCFCateveryfluidstepwillonlydoubletheoveralltimes, inspaceandtime,whichisalwaysagoodcompromisebetween butno significantimprovementhas beennoticed in the results. efficiency,accuracy,androbustness. Finally,gridspacingwillbeconstantbothalongtheradialdirec- As discussed in the introduction, only in recent years sta- tionrandthepolarangleθ,sothenumberofpointsisenoughto ble numerical schemes for GRMHD in dynamical spacetimes specifythegridineachdirection. 9 N.BucciantiniandL.DelZanna:GRMHDinaxisymmetricdynamicalspacetimes:theX-ECHOcode hot (ρ 10 7, p/ρ 0.2) atmosphere in hydrostatic equilib- − ≃ ≃ rium(αh = const),inpressurebalanceatthesurfaceoftheNS. Contrarytoprevioustreatments,wheretheatmospherewasreset ateverytimesteptokeepitstationary,weleaveitfreetoevolve (collapse or expand) in response to the NS oscillations. Given its low density, the atmosphere has negligible feedback on the star. The simulation is performed using 625 grid points in the radialdomainr = [0,20],correspondingtoastarresolvedover 250 points. The evolution is followed for a time t = 1500 max correspondinginphysicalunitsto 7.5ms. ≃ Fig. 1 shows a comparison between the initial density pro- file and that at t , together with a plot of the central density max ρ asafunctionoftime.RelativevariationsofdensityintheNS c interior are of order of 10 4, with major deviations only at the − contactdiscontinuityoftheNSsurface,duetodiffusionoverthe muchlowerdensityatmosphere.Thistriggersthenaturalvibra- tionmodesoftheNS,thataretheobservedfluctuations,which arethenaturaloutcomeforthisphysicalsystem.Noticethatthe centraldensity,plottedintheinsertofFig.1,showsfluctuations oforderof10 3atmost,butnosignofanyseculartrend.Thisis − duetothelargenumberofpointsoverwhichthestarisresolved. Theslowdampingoftheoscillations,from10 3 toafew10 4, − − isduetothethermaldissipationassociatedtotheuseofanideal gasEoSandto thenumericalviscosityofthe scheme,whereas it is notpresentif a polytropicEoSis used (see nextSect. 5.2 foracomparisonbetweenthetwoEoS).Atthebeginningofthe simulation we observe a relaxation of the central density to a valuewhichis 2 10 4lowerthantheinitialcondition,prob- − ∼ × ablybecauseofdiscretizationerrors.Howevertheaveragevalue seems to remain constant at later times. This shows the ability of the code to maintain a stable equilibrium, even for several ( 10)soundcrossingtimes. Italsoshowsthatthe presenceof ∼ a dynamicalatmosphere causes no problem for the stability of theTOVsolution.Onthecontrary,theatmosphereitselfseems tobe quitestable, asshownbythefactthatitsdensitychanges byafactorsmallerthan1%. In the bottom panel of Fig. 1 we plot a Fourier transform ofthecentraldensityintime.Markersindicatethepositionsof theknowneigenmodes.Thisisatestoftheperformanceofthe Fig.1.Evolutionof astableTOVsolution inspherical symmetryand code in handling a dynamical spacetime, at least in the linear isotropic coordinates. The upper panel shows a comparison between densityintheinitialsolution(solidline)andtheresultatt = 1500 regime,for small perturbations.The values of the eigenmodes, max (diamonds).Forclaritytheresultat isshownevery5points.Thein- thefundamentalinparticular,arequitedifferentintheCowling max sertshowstheresiduals.Thespikeatr 8isduetodiffusiverelaxation approximationwherethemetriciskeptfixedintime(Fontetal. ≈ attheNSboundary. Themiddlepanelshowstherelativevariationsin 2002).Itisinterestingtonotealsothatnoinitialperturbationhas timeofthecentraldensity.ThelowerpanelshowstheFouriertransform beenintroduced,andthattheoscillationsofthestarareonlydue ofthethecentraldensity. Solidlineanddiamonds indicatethepower totherelaxationoftheinitialconditionsandpossibleroundoff- oftheFourierseriesinarbitraryunits.Theverticalmarkersindicatethe errors. Indeed, the large power present in the higher frequency frequencyofknowneigenmodes.Thefrequencyresolutionofourtime modes suggests that the initial excitation is confined to small seriesis 150Hz. ∼ scales, most likely at the surface of the star. The presence of a freelyevolvingatmospheredoesnotseemtoaffectthefrequency ofthemodes,atleastwithintheaccuracyofourtemporalseries. 5.1.StabilityofaTOVstableradialsolution Our first test consists in the evolution of a stable 1D radial 5.2.MigrationofanunstableTOVradialsolution NS configuration. We adopt, as initial condition, a solution of the TOV equations in isotropic coordinates, corresponding to Agenuinelynon-trivialsituationinthefullynon-linearregime, a polytropic gas with K = 100,n = 1, and central density involvinglargevariationsofthemetricandfluidstructure,isthe ρ = 1.280 10 3, a model also known as A0 (AU0) or B0 migrationofanunstable1DTOVsolution.Following(Fontetal. c − × (BU0)intheliterature(Fontetal.2000;Stergioulasetal.2004; 2002; Cordero-Carriońet al. 2009;Bernuzzi & Hilditch 2010) Dimmelmeieretal. 2006).Thiscorrespondstoa NSextending we select a solution of the TOV equations in isotropic coordi- to a radius r = 8.13. It is possible to show that this star lies nates,correspondingtoapolytropicgaswithagainK =100,n= on the stable part of the mass-radius curve, so we expect the 1,andcentraldensityρ = 7.993 10 3.Thiscorrespondstoa c − × codeto beabletomaintainthisconfigurationsfortimeslonger star extendingto a radiusr = 4.26, on the unstable partof the that their typical sound crossing time ( 0.5 ms). Outside the mass-radiuscurve.Theexternalatmosphereissetasinthepre- ∼ star we assume, at the beginningof the run, a low density and vioustest. Due to truncationerrorsand the initialrelaxationof 10