Computational General Relativistic Force-Free Electrodynamics: I. Multi-Coordinate Implementation and Testing
AAstronomy & Astrophysics manuscript no. A_code_paper c (cid:13)
ESO 2020July 16, 2020
Computational General Relativistic Force-Free Electrodynamics: I.Multi-Coordinate Implementation and Testing
J. F. Mahlmann (cid:63) , M. A. Aloy , V. Mewes , , P. Cerdá-Durán Departament d’Astronomia i Astrofísica, Universitat de València, 46100, Burjassot (Valencia), Spain National Center for Computational Sciences, Oak Ridge National Laboratory, P.O. Box 2008, Oak Ridge, TN 37831-6164, USA Physics Division, Oak Ridge National Laboratory, P.O. Box 2008, Oak Ridge, TN 37831-6354, USAReceived
Month Day, Year ; accepted
Month Day, Year
ABSTRACT
General relativistic force-free electrodynamics is one possible plasma-limit employed to analyze energetic outflows in which strongmagnetic fields are dominant over all inertial phenomena. The amazing images of black hole shadows from the galactic center and theM87 galaxy provide a first direct glimpse into the physics of accretion flows in the most extreme environments of the universe. Thee ffi cient extraction of energy in the form of collimated outflows or jets from a rotating BH is directly linked to the topology of thesurrounding magnetic field. We aim at providing a tool to numerically model the dynamics of such fields in magnetospheres aroundcompact objects, such as black holes and neutron stars. By this, we probe their role in the formation of high energy phenomenasuch as magnetar flares and the highly variable teraelectronvolt emission of some active galactic nuclei. In this work, we presentnumerical strategies capable of modeling fully dynamical force-free magnetospheres of compact astrophysical objects. We provideimplementation details and extensive testing of our implementation of general relativistic force-free electrodynamics in Cartesianand spherical coordinates using the infrastructure of the E instein T oolkit . The employed hyperbolic / parabolic cleaning of numericalerrors with full general relativistic compatibility allows for fast advection of numerical errors in dynamical spacetimes. Such fastadvection of divergence errors significantly improves the stability of the general relativistic force-free electrodynamics modeling ofblack hole magnetospheres. Key words.
Magnetic fields - Methods: numerical - Plasmas
1. Introduction
Relativistic, magnetically dominated plasma is a basic elementof the environments around neutron stars (Goldreich & Ju-lian 1969; Michel 1973; Scharlemann & Wagoner 1973; Con-topoulos et al. 1999), especially around magnetars (Lamb 1982;Thompson & Duncan 1995; Beloborodov & Thompson 2007),black holes (BHs; Blandford & Znajek 1977) including the ac-cretion disks surrounding them (MacDonald & Thorne 1982;Thorne et al. 1986; Beskin 1997), and their outflows (Takahashiet al. 1990; Lee et al. 2000; Punsly 2001; Lyutikov 2009). Theinterest for the environments surrounding neutron stars and BHshas sparked recently due to the overwhelming amount of newobservations available, e.g., in supermassive BHs (Event Hori-zon Telescope Collaboration et al. 2019a,b) or magnetars (Tur-olla et al. 2015; Kaspi & Beloborodov 2017). If magnetic fieldsdominate the dynamics, all plasma inertial and thermal contri-butions may be neglected. Thus, the only role of the plasma isto support the electromagnetic fields. Magnetic fields govern allplasma dynamics; the currents are not merely induced by thedrift of the matter distribution but completely determined by theelectromagnetic fields. Under the aforementioned conditions theplasma becomes force-free (e.g. Uchida 1997).The correct interpretation of recent breakthrough observa-tions requires building up a solid theoretical understanding ofthe astrophysical scenarios mentioned above. Due to their com-plexity and system size, they are well-suited for numerical ap-proaches. To model astrophysical plasma numerically under (cid:63) [email protected] force-free conditions, two principal formulations have emergedin recent years (see a detailed review in Paschalidis & Shapiro2013). Komissarov (2004) suggests the time-evolution of thefull set of Maxwell’s equations, where the magnetic inductionand displacement encode the general relativistic spacetime ge-ometry as non-vacuum e ff ects. This formulation has also beenemployed in an implementation relying on spectral methods(P hedra , Parfrey et al. 2015, 2017). Furthermore, Palenzuelaet al. (2010); Carrasco & Reula (2017) carried out simulationsin spherical geometries and higher-order finite di ff erence ap-proximations. McKinney (2006) introduces a formulation that isbased on an adaptation of general relativistic magnetohydrody-namics (GRMHD) to evolve the magnetic field as well as Poynt-ing fluxes in time. As such, it was implemented, for example,in the GRMHD code H arm (Gammie et al. 2003), and in theG i R a FFE code (Etienne et al. 2017) provided in the E instein T oolkit (Lö ffl er et al. 2012). For this project, we have imple-mented the Maxwell’s equations evolution system in general rel-ativity using the infrastructure of the E instein T oolkit . To thisend, we developed a new code for the numerical integration ofthe equations of general relativistic force-free electrodynamics(GRFEE) in dynamically evolving spacetimes. This tool has al-ready been applied to the study of the dynamics in magneto-spheres around compact objects (Mahlmann et al. 2019, 2020b).In this series of papers, we will extensively review implementa-tion details and characterize the numerical properties of our newcode. http: // a r X i v : . [ phy s i c s . c o m p - ph ] J u l & A proofs: manuscript no. A_code_paper
The E instein T oolkit infrastructure, which was originallydesigned for Cartesian coordinates, has recently been adapted tosupport spherical coordinates (Mewes et al. 2018, 2020). For cer-tain applications in the realm of astrophysical compact objects,it is beneficial to exploit coordinates that reflect the approximatesymmetries these systems’ possess. Spherical coordinates pro-vide such a coordinate system that enhances the accuracy of theemployed method. Starting from the so-called reference metricapproach, we write the evolution equations such that they corre-spond to conservation laws in a conformally related metric. Wethen use suitable finite volume discretizations (in Cartesian andspherical coordinates, Cerdá-Durán et al. 2008) for the integra-tion of intercell fluxes in our time-marching scheme. Alternativeapproaches have been employed, e.g., by Montero et al. (2014);Mewes et al. (2020) wherein all information about the underly-ing coordinate system is encoded in geometrical source terms.This work is organized as a series of papers. This manuscript(Paper I) reviews the general theory of GRFFE and our imple-mentation using the infrastructure of the E instein T oolkit . Thesecond paper of this series (Mahlmann et al. 2020a, Paper II) fo-cuses on the characterization of the numerical di ff usivity of ouralgorithm. Sections 2.1 and 2.2 lay out the theory backgroundof GRFFE and introduce an augmented conservative system ofpartial di ff erential equations (PDEs). We discuss the implemen-tation of this system in our scientific code in the E instein T oolkit in Sect. 3. Di ff erent finite volume integrators used for full sup-port of both Cartesian and spherical coordinates are reviewedin Sects. 3.1.1 and 3.1.2. We discuss two key ingredients forthe successful numerical integration of GRFFE, namely the con-servation of force-free constraints and the cleaning of numer-ical errors, in Sects. 3.3 and 3.5. Section 4 assembles a suiteof tests for the numerical calibration and characterization of ourcode. We demonstrate its ability to reproduce the basic dynamicsof force-free configurations (Sect. 4.1). We further demonstratethe code’s potential for the modeling of astrophysical plasma inmagnetar and BH magnetospheres in Sect. 5. Finally, we outlinedistinct features of the presented methods as well as implicationsfor GRFFE schemes in general in Sect. 6.
2. General Relativistic Force-Free Electrodynamics
The following sections as well as the code implementation inthe E instein T oolkit employ units where M (cid:12) = G = c = M (cid:12) ≡ . × − s ≡ .
98 m. This unit system is a variation ofthe so-called system of geometrized units (as introduced in ap-pendix F of Wald 2010), with the additional normalization of themass to 1 M (cid:12) (see also Mahlmann et al. 2019, on unit conversionin the E instein T oolkit ). In the following, Latin indices denotespatial indices, running from 1 to 3; Greek indices denote space-time indices, running from 0 to 3 (0 is the time coordinate). TheEinstein summation convention is used. To numerically solve the field equations of general relativity,a fully covariant formulation is not optimal. Instead, to arriveat a Cauchy initial value problem that can be evolved forwardin time, it is common to introduce a 3 + g µν , is foliated with a set ofnon-intersecting timelike hypersurfaces Σ t , i.e., level surfaces ofthe scalar field t (denoting the time coordinate). We denote the future-pointing, timelike normal on Σ t as n µ . It is defined throughthe constituting relation α n µ ∇ µ t = , (1)where ∇ µ denotes the spacetime covariant derivative built from g µν . The lapse function α indicates the separation in proper timebetween two hypersurfaces. The spatial three-metric γ i j is theprojection of the spacetime metric g µν onto Σ t : γ i j = ( g µ i + n µ n i )( g ν j + n ν n j ) g µν = g i j . (2)Trajectories of constant spatial coordinates across di ff erent hy-persurfaces Σ t define the time vector along them: t µ = α n µ + β µ . (3) β µ is the shift four-vector, which denotes the spacelike displace-ment in the direction perpendicular to n µ , required to reach theoriginal base coordinate in a hypersurface Σ t (cid:48) after leaving Σ t .The shift vector satisfies β µ n µ = t µ = (1 , , , n µ and its (metric) dual n µ (assuming the metric signatureis +
1) can be expressed in terms of lapse and shift as follows: n µ = (cid:32) α , − β i α (cid:33) , n µ = ( − α, , , . (4)The line element of the spacetime may be written in terms of thelapse, shift, and spatial metric in the 3 + ds = − α dt + γ i j (cid:16) dx i + β i dt (cid:17) (cid:16) dx j + β j dt (cid:17) . (5)The spacetime metric g µν is given by: g µν = (cid:32) − α + β i β i β j β i γ i j (cid:33) . (6)In this foliation, the Einstein equations can be cast into a setof evolution and constraint equations (see, e.g. Alcubierre 2008;Baumgarte & Shapiro 2010; Gourgoulhon 2012, for textbook in-troductions). One of the most widely used formulations to nu-merically solve the Einstein equations in 3 + γ i j and the conformalfactor e φ , which are related by¯ γ i j = e − φ γ i j , (7)where the conformal factor e φ can be written as e φ = ( γ/ ¯ γ ) . (8) γ and ¯ γ are the determinants of the spatial and conformallyrelated metric, respectively. The BSSN formalism also intro-duces the conformally related extrinsic curvature and the confor-mal connection functions as evolved variables. Throughout thiswork, we fix ¯ γ to be constant in time (the so-called Lagrangianchoice, Brown 2005): ∂ t ¯ γ = , (9)Thus, the time dependence of the spatial metric determinant isencoded only in the corformal factor, as √ γ = e φ √ ¯ γ . Keeping¯ γ fixed to its initial value simplifies expressions in the GRFFEevolution equations. This choice is particularly useful when inte-grating GRFFE in spherical coordinates, as we elaborate below. Article number, page 2 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing
The evolution of the full set of Maxwell’s equations is one possi-bility to deal with electrodynamics in general relativity (Komis-sarov 2004): ∇ ν F µν = I µ , (10) ∇ ν ∗ F µν = , (11)where F µν is the Maxwell tensor and ∗ F µν is its dual, defined as: ∗ F µν ≡ η µνλζ F λζ , (12)where η µνλζ = − ( − g ) − / [ µνλζ ] , η µνλζ = ( − g ) / [ µνλζ ] . (13)[ µνλζ ] is the completely antisymmetric Levi-Civita symbol with[0123] = + g the determinant of the spacetime metric g µν . I µ is the electric current four-vector associated with the chargedensity ρ = − n µ I µ = α I t , and the current three vector J i = α I i .A covariant definition of the current four-vector is (Komissarov2004) J µ = I [ ν t µ ] n ν , (14)where we use the standard terminology I [ ν t µ ] = ( I ν t µ − I µ t ν ).Note that ρ is the charge density as measured by the normal (Eu-lerian) observer defined by n µ . The current density 3-vector asmeasured by the Eulerian observer is the projection of I µ ontothe hypersurface Σ t : j i = ( g i µ + n i n µ ) I µ = α − ( J i + ρβ i ) . (15)Komissarov (2004) introduces a set of vector fields, which areanalog to the electric and magnetic fields, E and B , and electricdisplacement and magnetic induction, D and B , of the electro-dynamic theory of continuous media (see, e.g., Jackson 1999,§I.4). They have the following spatial components in a 3 + E i = F it , (16) B i = e i jk F jk , (17) D i = e i jk ∗ F jk , (18) H i = ∗ F ti . (19)Following the convention in, e.g., Baumgarte & Shapiro (2003),the four-dimensional volume element induces a volume elementon the hypersurfaces of the foliation: e abc = n µ η µ abc = − αη abc . (20)In the previous expression, e abc (cid:44) e i jk = − αη i jk = [ i jk ] / √ γ .Using the definitions (17) and (18) in the time componentsof Eqs. (10) and (11), one obtains the familiar constraintsdiv D = ρ, (21)div B = . (22) Another evolution scheme, developing energy fluxes rather than elec-tric fields, was employed by e.g., McKinney (2006) or Etienne et al.(2017).
We separately evolve the charge continuity equation, which isobtained from the covariant derivative of Eq. (10), ∇ ν I ν = , (23)to ensure the conservation of the (total) electric charge in thecomputational domain, as well as the compatibility of the chargedistribution obtained numerically, with the currents that theygenerate. If this is not done, the di ff erence | div D − ρ | may growunbounded with time due to the accumulation of numerical er-rors (Munz et al. 1999).Embodied in the definitions (16) to (19) one finds the follow-ing vacuum constitutive relations (Komissarov 2004): E = α D + β × B , (24) H = α B − β × D . (25)We may now write the Maxwell tensor as measured by the Eule-rian observer in terms of the electric, D µ , and magnetic, B µ , fieldfour-vectors (cf. McKinney 2006; Antón et al. 2006): F µν = n µ D ν − D µ n ν − η µνλζ B λ n ζ , (26) ∗ F µν = − n µ B ν + B µ n ν − η µνλζ D λ n ζ , (27)which satisfy D µ = F µν n ν = (cid:16) , D i (cid:17) , (28) B µ = ∗ F νµ n ν = (cid:16) , B i (cid:17) . (29)For later reference, we provide two Lorentz invariants of theFaraday tensor, namely: ∗ F µν F µν = D µ B µ , (30) F µν F µν = B − D ) . (31)In order to build up a stationary magnetic configuration (as, e.g.,in the magnetosphere around a compact object), it is necessary toguarantee that there are either no forces acting on the system or,more generally, that the forces of the system are in equilibrium.Except along current sheets the latter condition implies that theelectric four-current I µ satisfies the force-free condition (Bland-ford & Znajek 1977): F µν I ν = . (32)Eq. (32) is equivalent to a vanishing Lorentz force density f µ onthe charges measured by the Eulerian observer: ∇ ν T νµ = − F µν I ν = − f µ ≡ . (33)Also, Eq. (32) can be seen as a system of linear equations, thenon-trivial solution of which demands that the determinant of F µν vanishes. Since det F µν = ( ∗ F µν F µν ) / = ( D µ B µ ) , theforce-free condition (32) reduces to ∗ F µν F µν = , (34)or, equivalently (see Eq. 30) D µ B µ = . (35)Hence, the component of the electric field parallel to the mag-netic always vanishes. Since det F µν =
0, the rank of F µν (re-garded as a 4 × F µν has non-vanishingcomponents. If a µ is a zero eigenvector of F µν , i.e., F µν a µ = a µ is b µ = ∗ F µν a ν ,and the Faraday tensor can be expressed as F µν = η µνλδ a λ b δ (cf. Article number, page 3 of 19 & A proofs: manuscript no. A_code_paper
Komissarov 2002). Hence, it admits a two-dimensional spaceof eigenvectors associated with the null eigenvalue (cf. Uchida1997). These zero eigenvectors are time-like if the Lorentz in-variant F µν F µν is positive (Uchida 1997). The sign of the invari-ant F µν F µν is not unanimously defined for generic electromag-netic four-vectors B µ and D µ . To chose the sign of the invariant, itis useful to consider the force-free approximation as a low inertialimit of relativistic MHD. This means that a physical force-freeelectromagnetic field should be compatible with the existence ofa velocity field of the plasma. Recalling that the plasma four-velocity u µ is a unit time-like vector ( u µ u µ = − f µ ∝ F µν u ν , a physical force-free electromag-netic field ( f µ =
0) should satisfy F µν u ν = F µν F µν (see Eq. 31) should consistently bepositive, i.e., F µν F µν = B − D ) > . (36)In the introduced language of the full system of Maxwell’s equa-tions in 3 + D · B = , (37) B − D ≥ . (38)Condition (38) implies that the magnetic field is always strongerthan the electric field. Equivalently, one can classify the degen-erate electromagnetic tensor as magnetic, since condition (36)guarantees that there exists a frame in which an observer at restmeasures zero electric field (cf. Uchida 1997). This observer isthe comoving observer with four-velocity u µ in the ideal MHDlimit.The challenge of maintaining the physical constraints ofdiv B = D = ρ in numerical simulations has beenreviewed throughout the literature (e.g., Dedner et al. 2002;Mignone & Tzeferacos 2010), and applied to GRFFE, e.g., byKomissarov (2004) and the relativistic MHD regime, e.g., byPalenzuela et al. (2009); Miranda-Aranguren et al. (2018). Fol-lowing Palenzuela et al. (2009, 2010) as well as Mignone & Tze-feracos (2010), we suggest to modify the system of Maxwell’sequations (Eqs. 10 and 11) in the following way (cf. Alic et al.2012): ∇ ν ( F µν + g µν Φ ) = I µ + t µ κ Φ Φ , (39) ∇ ν ( ∗ F µν + s µν Ψ ) = t µ κ Ψ Ψ . (40)Here, the definition of t µ is given in Eq. (3), and we define s µν ≡ c h γ µν − n µ n ν . c h corresponds to a speed of propagationof the divergence cleaning errors; κ Φ and κ Ψ are adjustable con-stants that control the parabolic damping of the aforementionednumerical errors. The scalar potentials Ψ and Φ are ancillaryvariables employed to control the errors in the elliptic constrainsdiv B = D = ρ , respectively. This is implemented inpractice by augmenting the system of Maxwell’s equations withextra evolution equations for Φ and Ψ (see Sect. 3.5). The aug-mented system of Maxwell Equations (Eqs. 39 and 40), can bewritten as a system of balance laws of the form ∂ t C + ¯ ∇ j F j = S n + S s , (41)where ¯ ∇ is the covariant derivative with respect to the confor-mally related metric, ¯ γ (Eq. 7). C denotes the vector of conservedvariables, F j the flux vectors, S n the geometrical and current-induced source terms, and finally S s are the potentially sti ff source terms (cf. Komissarov 2004, App. C2). Note that each of these quantities consists of elements in a multidimensionalspace. In general, the conserved variables are derived from theso-called primitive variables; primitive variables are usually thequantities measured by the Eulerian observer, namely ρ , B , and D , as well as the numerical cleaning potentials Ψ and Φ . Adapt-ing the notation used by Cerdá-Durán et al. (2008) and Monteroet al. (2014), we specify the components of Eq. (41) in terms ofthe determinant of a reference metric ˆ γ . We define the conservedvariables as C ≡
LQP b i d i = e φ (cid:114) ¯ γ ˆ γ ρ Ψ α Φ α B i + Ψ α β i D i − Φ α β i , (42)with their corresponding fluxes F j ≡ e φ (cid:114) ¯ γ ˆ γ α J j B j − Ψ α β j − (cid:16) D j + Φ α β j (cid:17) e i jk E k + α (cid:16) c h γ i j − n i n j (cid:17) Ψ − (cid:16) e i jk H k + αg i j Φ (cid:17) . (43)For the source terms, the split according to equation (41) yieldsthe source terms S n , and the potentially sti ff source terms S s : S n ≡ e φ (cid:114) ¯ γ ˆ γ α ΨΓ t µν s µν α ΦΓ t µν g µν − ρ − α ΨΓ i µν s µν α ΦΓ i µν g µν − J i , (44) S s ≡ e φ (cid:114) ¯ γ ˆ γ − ακ Ψ Ψ − ακ Φ Φ . (45)In the previous expressions, Γ αβγ are the Christo ff el symbols ofthe Levi-Civita connection associated with the spacetime metic g µν . In both Cartesian and spherical coordinates, we always makethe initial choice ¯ γ ( t = = ˆ γ , so that, due to the Eq. (9), the ratio¯ γ/ ˆ γ = γ i j and iscontinuously enforced in the spacetime evolution by making thereplacement¯ γ i j → (ˆ γ/ ¯ γ ) ¯ γ i j (46)at every sub-step of the time integration. Article number, page 4 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing
In force-free electrodynamics there is no uniquely defined restframe for the fluid motion (e.g., Uchida 1997; McKinney 2006;Paschalidis & Shapiro 2013; Shibata 2015); the electromag-netic current I µ cannot be determined by tracking the veloc-ity of charges throughout the domain. Rather, the enforcementof the force-free conditions (37), and (38) determines a suit-able current. The conservation condition (implicitly embodiedin Maxwell’s equations) L n ( D · B ) = n µ ∇ µ ( D · B ) = , (47)where L n is the Lie derivative with respect to n µ is equivalent to ∂ t ( D · B ) =
0. Together with conditions (37) and (38), it can becombined to obtain an explicit expression for the so-called force-free current I µ ff (cf. McKinney 2006; Komissarov 2011; Parfreyet al. 2017): I µ ff = ρ n µ + ρ B η νµαβ n ν D α B β + B µ B η αβλσ n σ (cid:16) B λ ; β B α − D λ ; β D α (cid:17) . (48)In practice, the combination of the force-free current (48) as asource-term to Eq. (10) - or Eq. (39) if we consider the aug-mented system of equations - with numerically enforcing con-ditions (37) and (38) restricts the evolution to the force-freeregime. The discussion of techniques to ensure a physical (cf.McKinney 2006) evolution of numerical force-free codes is a re-current topic that can be found throughout the literature (e.g.,Lyutikov 2003; Komissarov 2004; Palenzuela et al. 2010; Alicet al. 2012; Paschalidis & Shapiro 2013; Carrasco & Reula 2017;Parfrey et al. 2017; Mahlmann et al. 2019). We review one ofthese techniques in Sect. 3.3.
3. Numerical methodology
Our GRFFE method uses the framework of the E instein T oolkit (Lö ffl er et al. 2012). The E instein T oolkit is an open-sourcesoftware package utilizing the modularity of the C actus code(Goodale et al. 2003), which enables the user to specify so-called thorns in order to set up customized simulations andprovides (basic) adaptive mesh refinement (AMR) via the C ar - pet driver Goodale et al. (2003); Schnetter et al. (2004). Thespacetime evolution is performed using the M ac L achlan thorn(Brown et al. 2009) as an implementation of the BSSN for-malism. Recently, numerical relativity in spherical grids hasbeen successfully enabled on the traditionally Cartesian E instein T oolkit by the new implementation of S pherical NR (Meweset al. 2018, 2020), which is built upon a reference-metric formu-lation of the BSSN equations (Brown 2009; Montero & Cordero-Carrión 2012; Baumgarte et al. 2013). We make use of a varietyof open-source thorns within the E instein T oolkit , such as theapparent horizon finder AHF inder D irect (Thornburg 2004), theextraction of quasilocal quantities Q uasi L ocal M easures (Dreyeret al. 2003), and the e ffi cient S ummation B y P arts thorn (Dieneret al. 2007).In our code, the time update of the system of partial di ff er-ential equations (see Eq. 41) is done applying the method-of-lines (e.g., LeVeque 2007) implemented in the E instein T oolkit thorn M o L. For the numerical test shown in this paper we make http: // https: // bitbucket.org / eschnett / carpet / src / master / http: // / ~eschnett / McLachlan / use of the fourth-order accurate (not strictly TVD) Runge-Kuttamethod implemented in the thorn M o L.To ensure the conservation properties of the algorithm, itis critical to employ refluxing techniques correcting numericalfluxes across di ff erent levels of mesh refinement (e.g., Collinset al. 2010). Specifically, we make use of the thorn R efluxing in combination with a cell-centered refinement structure (cf. Shi-bata 2015). We highlight the fact that employing the refluxingalgorithm makes the numerical code 2 − We solve Eq. (41) by discretizing its integral over the volume V of a cell of our numerical mesh (cf. LeVeque 2007; Mignone2014; Martí & Müller 2003), ∂ t (cid:104)C(cid:105) + V (cid:90) ∂ V d A · F = (cid:104)S n (cid:105) + (cid:104)S s (cid:105) . (49)Here, (cid:104)(cid:105) denotes the volume average of a quantity. The diver-gence term ¯ ∇ j F j appearing in Eq. (41) is integrated by apply-ing Stoke’s theorem and summing up the reconstructed fluxes F through the cell interfaces with respective area elements d A .In practice, we approximate volume averages by cell-centered values for each grid element. We identify each of theseelements by the indices ( i , j , k ) which correspond to the locations x i = x + i ∆ x , y j = y + j ∆ y , and z k = z + k ∆ x . ∆ x , ∆ y , and ∆ z represent the (uniform) numerical grid spacing in each coordi-nate direction. The quantities ( x , y , z ) denote the coordinatesof an arbitrary reference point in 3D. Face-centered quantitiesare indicated by the subscript of a half-step added to the respec-tive index, i.e., the subscript i + / i , j , k ) and ( i + , j , k ). If nosubscript is provided, we refer to the cell-centered values. The system of Eqs. (42) to (45) is specified to its application inCartesian coordinate systems ( x , y, z ) by setting (cid:112) ˆ γ =
1. In thiscase, the cell volume is V = ∆ x × ∆ y × ∆ z , (50)and the area elements are denoted byd A = ( ∆ y × ∆ z , ∆ x × ∆ z , ∆ x × ∆ y ) . (51)Eq. (49) is approximated by evaluating the fluxes F as recon-structed averages at cell interfaces:1 V (cid:90) ∂ V d A · F ≈ F xi + / − F xi − / ∆ x + F y j + / − F y j − / ∆ y + F zk + / − F zk − / ∆ z . (52) Refluxing at mesh refinement interfaces by Erik Schnetter: https: // svn.cct.lsu.edu / repos / numrel / LSUThorns / Refluxing / trunkArticle number, page 5 of 19 & A proofs: manuscript no. A_code_paper
In spherical coordinates ( r , θ, φ ), (cid:112) ˆ γ = r sin θ , and the cell vol-ume is V = − ∆ r × ∆ cos θ × ∆ φ, (53)where ∆ r = r i + / − r i − / and ∆ cos θ = cos θ j + / − cos θ j − / .The numerical stability of the spacetime integral in Eq. (49) crit-ically depends on the balancing of coordinate singularities, e.g.,at the polar axis and the origin of the spherical coordinate sys-tem. We guarantee an exact evaluation of metric contributionsat the location of the cell-interfaces by transforming the recon-structed fluxes F to an orthonormal basis. The area elements inan orthonormal basis are denoted byd ˆ A = (cid:16) r sin θ × ∆ θ × ∆ φ, r sin θ × ∆ r × ∆ φ, r × ∆ r × ∆ θ (cid:17) . (54)Eq. (49) is approximated by evaluating the fluxes F as recon-structed averages in an orthonormal basis at cell interfaces:1 V (cid:90) ∂ V d ˆ A · ˆ F ≈ r i + / ˆ F ri + / − r i − / ˆ F ri − / ∆ r − ∆ r ∆ r sin θ j + / ˆ F θ j + / − sin θ j − / ˆ F θ j − / ∆ cos θ − ∆ r ∆ r ∆ θ ∆ cos θ ˆ F φ k + / − ˆ F φ k − / ∆ φ . (55)In analogy to the above, we use ∆ r = r i + / − r i − / . The recon-structed fluxes F (coordinate basis) are related to their orthonor-mal counterparts ˆ F by the following relations:ˆ F r = F r , ˆ F θ = r × F θ , ˆ F φ = r × sin θ × F φ . (56) We employ an approximate (HLL) Riemann solver (e.g., Rez-zolla & Zanotti 2013) to derive the numerical fluxes at the cellinterfaces: F j = λ + F j ( U − ) − λ − F j ( U + ) + λ + λ − ( U + − U − ) λ + − λ − . (57) U + , and U − correspond to the reconstructed (conserved) vari-ables at the cell interfaces. λ ± are given by the minimal or maxi-mal wave speeds: λ + = max (0 , w ) , λ − = min (0 , w ) . (58)In flat space, the propagation speeds for the conservative schemederived from Eqs. (10) and (11) are λ + = λ − = −
1. Charac-teristic speeds of the force-free electrodynamics equations havebeen obtained, e.g., by Komissarov (2002) and Carrasco & Reula(2016). For the GRFFE system of Eqs. (42) to (45), the charac-teristic speeds w are w = − β i ± α (cid:112) γ ii m = − β i ± c h α (cid:112) γ ii m = w q m = . (59) Here, we do not employ the summation convention; by m wedenote the multiplicity of the respective eigenvalues. The speedsEVI correspond to the coordinate velocity of light as defined byCordero-Carrión et al. (2008). The other two eigenspeeds (EVII)account for the propagation of the divergence cleaning potentialsat speed c h . Finally, EVIII corresponds to the wavespeed inducedby the continuity equation of charge conservation, which is atmost the coordinate velocity of light (EVI). Across the literature (e.g., Komissarov 2004; Alic et al. 2012;Parfrey et al. 2017) we find various modifications in the defini-tion of I µ to drive the numerical solution of the system of PDEs(Eqs. 10 and 11) toward a state which fulfills the magnetic dom-inance condition (38) by introducing a suitable cross-field con-ductivity. This e ff ectively augments condition (47) used to deter-mine expression (48) by a recipe to avoid (numerically) buildingup violations of the perpendicularity condition (37).A straightforward way to guarantee the preservation of the D · B = D ) is projected onto the direction perpendicular tothe magnetic field ( B ) in every point of the numerical mesh (e.g.,Palenzuela et al. 2010): D i → D k (cid:32) δ i k − B k B i B (cid:33) . (60)Alternatively, dissipative currents (induced by so-called driverterms) may ensure the evolution of the electromagnetic fieldstowards physically allowed (force-free) configurations. Usingdriver terms, the numerical evolution does not guarantee thatthe electromagnetic fields are exactly force-free after every time-step. However, force-free constraint violations are significantlyreduced. While Komissarov (2004, 2011), and Alic et al. (2012)introduce a modified Ohm’s law with a suitably chosen cross-field dissipation, Parfrey et al. (2017) modify the force-free cur-rents in Eq. (48) with additional dissipation driving the evolutiontowards a target ( D · B = L n ( B · D ) = κ I (cid:16) α − η J · B − D · B (cid:17) . (61)Here, κ I is the decay rate driving the left-hand side of Eq. (47)toward the target value and η is a dissipation coe ffi cient for theelectric field, which is parallel to the current.As for the preservation of the B − D ≥ D ) is rescaled in every point of the nu-merical mesh to the length of the magnetic field ( B ) in a qualita-tively similar manner as in Eq. (60): D i → D i (cid:32) − Θ ( χ ) + | B || D | Θ ( χ ) (cid:33) , (62)where Θ is the Heaviside function, and χ = D − B . Again,an alternative is employed by Komissarov (2011), and Alic et al.(2012), introducing driver terms for additional dissipative cur-rents, also for the conservation of the B − D ≥ Article number, page 6 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing
Our GRFFE scheme employs, by default, the algebraic cor-rection of electric fields in every (intermediate) step of the timeevolution as given by Eqs. (60) and (62). However, in Mahlmannet al. (2019) we resorted to a suitably chosen resistivity model (inanalogy to Komissarov 2004) replacing the instant algebraic cut-back of the electric displacement field by a gradual relaxation offorce-free violations. For a review on the interpretation of con-straint violations in GRFFE, we refer to Mahlmann et al. (2019).
The last term in Eq. (48) is the component of the current parallelto the magnetic field, with the spatial projection j || = B · ( ∇ × B ) − D · ( ∇ × D ) B B . (63)We have empirically found (see Paper II), that the discretiza-tion of the parallel current is one of the main sources of nu-merical di ff usivity in our code in certain tests. Indeed, the pres-ence of derivatives of conserved quantities in the parallel current(Eq. 48) makes the practical evaluation of this term cumbersomein numerical simulations. This has brought some authors (e.g.,Yu 2011) to ignore it completely (i.e., assuming j || = ), and re-moving the accumulated parallel component of the electric fieldemploying an algebraic procedure akin to that of Eq. (60). How-ever, this specific strategy of Yu (2011) did not yield satisfactoryresults when employed with the numerical framework describedin this paper.The order of the discretization of the derivatives has to becomparable with the order of accuracy of the spatial reconstruc-tion. Otherwise, the global order of accuracy of the scheme de-creases (see Sect. 2.1 of Paper II). In the previous applicationsof our code (Mahlmann et al. 2019, 2020b) and independentlyof the order of the spatial reconstruction, we have employed afourth-order accurate discretization of the partial derivatives inEq. (48). In case of the (exemplary) sweep in the x direction,the finite di ff erence discretization is of the following form (foruniform grids): (cid:34) ∂ A ∂ x (cid:35) thi ≈ A i − − A i − + A i + − A i + × ∆ x , (64)where A denotes a quantity on the numerical mesh (e.g., D or B )and the respective locations are labeled as in Sect. 3.1. In PaperII, we will evaluate the improvements by changing the discretiza-tion of j || according to the sixth-order and eighth-order accurateformulae (cid:34) ∂ A ∂ x (cid:35) thi ≈ − A i − + A i − − A i − + A i + − A i + + A i + × ∆ x , (65) (cid:34) ∂ A ∂ x (cid:35) thi ≈ A i − − A i − + A i − − A i − × ∆ x + A i + − A i + + A i + − A i + × ∆ x . (66) We extend the augmented evolution equations by a hyper-bolic / parabolic divergence error cleaning with the possibility ofhaving a hyperbolic advection speed, c h > ∇ µ yields for the simplified case of stationary spacetimes(cf. Komissarov 2004): − κ Ψ ∇ t Ψ = ∇ µ ∇ ν (cid:16) ∗ F µν + (cid:16) c h γ µν − n µ n ν (cid:17) Ψ (cid:17) = ∇ µ ∇ ν (cid:16) c h γ µν − n µ n ν (cid:17) Ψ= c h ∇ i ∇ i Ψ − ∇ µ ∇ ν n µ g να n α Ψ= c h ∇ i ∇ i Ψ + ∇ t ∇ t Ψ . (67)This compares to telegrapher equations, used for example to de-scribe signal propagation in lossy wires. In this analogy, κ Ψ , and c h are the parameters controlling the damping and advectionof numerical errors (Mignone & Tzeferacos 2010). We stressthe correspondence of c h with a finite propagation speed fordivergence errors (Mignone & Tzeferacos 2010) and their de-cay according to the damping factor κ Ψ . For c h chosen equalto the speed of light, Eq. (11) reduces to the evolution systemgiven in Palenzuela et al. (2009). In order to minimize viola-tions of div B = ≤ c h ≤
2. In practice, a propagationspeed within this interval does not limit the time step strongly,since the numerical evolution of the BSSN equations usuallydemands Courant–Friedrichs–Lewy (CFL) factors significantlysmaller than unity (say f cfl ∼ . − .
3) and, often, choosing c h > c h =
2) to significantlyreduce the numerical noise. We employ the same scheme with c h = Ψ . κ Ψ and κ Φ are damping rates, introducing time scales whichact in addition to the advection time scales of the hyperbolic con-servation laws of the augmented GRFFE system (Eqs. 42 to 45).In order to deal with the potential sti ff ness introduced by theparabolic damping numerically, we employ a time step splittingtechnique (Strang splitting, see, e.g., LeVeque 2007), which hasbeen applied previously to GRFFE by Komissarov (2004). Priorto and after the method-of-lines time integration , we evaluatethe equations P ( t ) = P exp (cid:104) − α κ Φ c h t (cid:105) , (68) Q ( t ) = Q exp (cid:104) − α κ Ψ t (cid:105) , (69)for a time t = ∆ t /
2. We find it beneficial to choose a largevalue for κ Φ , in some cases ∼ ff ectively dissipating chargeconservation errors on very short time scales (and justifying thetime-splitting approach). As for the divergence cleaning, we con-ducted a series of tests, optimizing κ Ψ to yield stable and con-verging evolution for all shown resolutions, ultimately resortingto κ Ψ = . − .
25 (see also Mahlmann et al. 2019).
4. Numerical Tests
We present several tests with results that specifically depend onthe various numerical methods (e.g., reconstruction, cleaning ofnumerical errors) available in our new code. Since the code isgenuinely 3D, in 1D and 2D simulations, the surplus dimensionsare condensed to the extension of one cell by applying appro-priate boundary conditions to them. Section 4.1 reviews the 1D Specifically, we add the exact evaluation of sti ff source terms beforethe scheduling bin M o L_S tep and before M o L_P ost S tep . The latter hasto be restricted to the last intermediate step of the method-of-lines inte-gration. Article number, page 7 of 19 & A proofs: manuscript no. A_code_paper tests of signal propagation and stability in GRFFE following thework by Komissarov (2004) and Yu (2011) closely. In Sect. 4.2we probe the correct representation of force-free plasma wave in-teractions by reproducing key results by Li et al. (2019). All testsin this section are performed in a fixed background Minkowskispacetime. The initial value of the charge density is computed as ρ = div D and the cleaning potentials are set to Ψ = Φ = GRFFE allows two modes of plasma waves (Komissarov 2002;Punsly 2003; Li & Beloborodov 2015; Li et al. 2019): Alfvénwaves which transport energy, charge, and current along mag-netic field lines and fast waves which correspond to the linearlypolarized waves of vacuum electrodynamics (see also Sect. 3.2).The following set of 1D problems is selected to demonstrate(a) the correct propagation of fast waves, (b) the formation ofa current-sheet when magnetic dominance breaks down and (c)the correct modeling of stationary Alfvén waves which do nottransport energy across magnetic field lines (cf. Li et al. 2019).The latter can only di ff use due to a finite numerical resistivity ifthe force-free constraints are not preserved (see Paper II). Thenumerical solution to all these problems critically depends onthe employed reconstruction algorithms. Since our code employsnumerical reconstruction in 1D sweeps across all dimensions, weconsider the following suite of 1D tests a fundamental measurefor the performance of our GRFFE scheme. We verify (in thesense of Roache 1997) the correct implementation of the recon-struction methods evaluating the convergence order from severaldata-sets with increasing resolution. Specifically, we evaluate the(global) di ff erence measure (cf. Antón et al. 2010) (cid:15) ab = N × (cid:88) i (cid:12)(cid:12)(cid:12) u ai − u bi (cid:12)(cid:12)(cid:12) , (70)where u a and u b are the one-dimensional vectors (of N elements)of the considered evolved quantity at di ff erent levels of resolu-tion, a , b ∈ [1 , , ∆ x , ∆ x , ∆ x , where ∆ x / ∆ x = ∆ x / ∆ x . With level1 being the level with the finest resolution, the (empirical) orderof convergence is then defined as (see also Bona et al. 1998): O = log (cid:16) (cid:15) /(cid:15) (cid:17) log ( ∆ x / ∆ x ) . (71)Unless stated otherwise, in the following tests we employ afourth-order accurate discretization of the parallel current j || . Komissarov (2004) examines two variations of a current sheetproblem, one of which has a solution in force-free electrodynam-ics, while the other violates the force-free constraints (Eqs. 37,and 38). The tests for physical current sheets (Fig. 1) and degen-erate current sheets (Fig. 2) are initialized by the following setof data: D = , B = (cid:16) . , B y , . (cid:17) , B y = (cid:40) B x < − B x > . (72)If B <
1, there exists a force-free solution given by two fastwaves traveling at the speed of light (see Fig. 1, also Fig. C2 inKomissarov 2004). For B > B − D to zero in acurrent sheet located at x =
0. At this location, the conservationof the force-free constraints (Sect. 3.3) becomes important forthe field dynamics, i.e., it changes the structure of the propagat-ing waves. The states bounded by the fast waves are terminatedat the current sheet and a standing field reversal remains (seeFig. 2, cf. Fig. C2 in Komissarov 2004). We take advantage ofthis test to compare the performance of two di ff erent reconstruc-tion schemes: The second-order accurate, slope limited TVD re-construction with a monotonized central (MC) limiter (van Leer1977), and the seventh-order accurate monotonicity preserving(MP7, Suresh & Huynh 1997) reconstruction. From the resultsof the presented tests (Figs. 1 and 2), we draw the following con-clusions: – Fast electromagnetic waves propagate correctly at the speedof light. – For a resolution similar to the one employed in Komissarov(2004), where ∆ x = . – For resolutions below the lowest presented resolution (i.e.,for ∆ x > .
05) the wave structure of the presented testquickly smears out. – Conservation of force-free constraints in the degenerate cur-rent sheet test is working well and agrees with similar teststhroughout the literature. – While monotonicity preserving reconstruction is slightlymore oscillatory than, e.g., monotonized central flux limiters,the higher-order schemes allow a steeper resolution of wave-fronts and current sheets. While the order of convergence ofthe (more di ff usive) MC reconstruction approaches the for-mal theoretical order of convergence ( O = – Although some degradation of the order of convergence isexpected in non-smooth regions of the flow (e.g., the discon-tinuities associated with fast or Alfvén waves), the algebraicenforcement of the violated force-free constraints seems tohave a large impact on the computed value of O . Very likely,the latter procedure is the main source of deviation from thetheoretical expectations regarding the order of convergence.Given the previous statements, the developed GRFFE codepasses the 1D (degenerate) current sheet test. Komissarov (2002), Yu (2011) and Paschalidis & Shapiro (2013)suggest the three-wave problem (or a variant thereof, see Fig. 3)as a test for force-free electrodynamics. The initial discontinu-ity at x = ff ectively combines the previously in-troduced test of Sect. 4.1.1 with the standing Alfvén wave testwhich was also employed by Komissarov (2004). The initialelectromagnetic field is given by (Paschalidis & Shapiro 2013): B = (1 . , . , . , D = ( − . , − . , . , if x < B = (1 . , . , . , D = ( − . , . , − . , if x > . (73)We evolve this setup in time and present the results for dif-ferent resolutions in Figs. 3 and 4. Around the standing Alfvénwave at x =
0, high-order reconstructions tend to develop small-scale oscillations, especially visible in the plots of D x , restrictedto the region delimited by the fast waves (at x = ± t = Article number, page 8 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing Δ x = Δ x = Δ x = - - - - - B y ≈ - - - - - - - - D z ≈ Δ x = Δ x = Δ x = - - - - - B y ≈ - - - - - - - - D z ≈ Fig. 1.
Current sheet test (Komissarov 2004; Yu 2011) as described bythe initial data in Eq. (72) on a x ∈ [ − ,
2] grid ( f cfl = .
25) at t = . B = . ff erent resolutions. Two fast waves emerge from theoriginal discontinuity and propagate outwards with the speed of light(analytical position of the waves are indicated by dashed vertical lines).The order of convergence, O is indicated according to Eq. (71). Top:
MP7 reconstruction.
Bottom:
MC reconstruction. Δ x = Δ x = Δ x = - - - - - B y ≈ - - - - - - - - - D z ≈ Δ x = Δ x = Δ x = - - - - - B y ≈ - - - - - - - - - D z ≈ Fig. 2.
Degenerate current sheet test (Komissarov 2004; Yu 2011) as asdescribed by the initial data (Eq. 72) with B = .
0. Two fast wavesemerge from the original discontinuity and propagate outwards with thespeed of light. The cross field conductivity (induced by conserving con-ditions 37, and 38) terminates the fast waves in the breakdown-zone.
Top:
MP7 reconstruction.
Bottom:
MC reconstruction.
Oscillations around this discontinuity can also be observed (forhigher resolutions) in part of the literature (specifically, Fig. 4in Yu 2011). The order of convergence is slightly reduced whencompared to the results shown in the previous section, proba-bly due to the specific challenges of resolving stationary Alfvén Δ x = Δ x = Δ x = - - - B y ≈ - - - - D y ≈ - - - B z ≈ - - - - - - - - - - D x ≈ Fig. 3.
Three-wave problem (Paschalidis & Shapiro 2013) as describedby the setup in Eq. (73). Numerical setup and labels are the same as inFig. 1. The initial discontinuity at x = Δ x = Δ x = Δ x = - - - B y ≈ - - - - D y ≈ - - - B z ≈ - - - - - - - - - - D x ≈ Fig. 4.
Same as Fig. 3 but employing MC reconstruction. waves, not only in GRFFE but also in relativistic MHD (see, e.g.Antón et al. 2010).The empirical order of convergence grows employing MP7in combination with a sixth-order accurate discretization of theparallel current (65). This growth increases from, e.g.
O ≈ . O ≈ .
8, for D x , but is negligible in other variables. In anycase, the overall numerical solution does not significantly changemodifying the order of accuracy of the calculation of j || in the 1Dproblems involving discontinuities.Komissarov (2004) achieves high accuracy maintaining asingle standing Alfvén wave stationary during evolution for reso-lutions comparable to the highest one shown in Figs. 3 and 4. The Article number, page 9 of 19 & A proofs: manuscript no. A_code_paper Δ x = Δ x = Δ x = - - B z ≈ - - B z ≈ Fig. 5.
Stationary Alfvén wave problem (Komissarov 2004), same nu-merical setup as in Fig. 1. The analytic solution (Eq. 74) is indicated bya gray line. numerical techniques of Komissarov (2004) are slightly di ff er-ent from ours, employing, for example, a linear Riemann solverwhich makes use of the full spectral decomposition of the FFEequations. The latter distinguishes all physical, and non-physicalwave speeds and may provide additional accuracy at critical lo-cations (in the context of GRFFE, e.g., at current sheets). Ad-ditionally, Komissarov (2004) employs a di ff erent form of thecurrent in Faraday’s Eq. (10) based on a specific (numerical) re-sistivity model to drive electromagnetic fields towards a force-free state throughout the evolution. Although one could suspectthat this di ff erent treatment of the currents may alter the numer-ical solution significantly in this test (which is dominated by thenumerical di ff usivity of the standing wave), we find that our re-sults are quite similar to the ones of Komissarov (2004) - but seea more detailed analysis in Paper II.Next, we consider the analytical solution of a standingAlfvén wave as initial data in the following: B = (1 , , B z ) , D = ( − B z , , , B z = x ≤ + .
15 [1 + sin [5 π ( x − . < x ≤ . . x > . . (74)We present the results of the Alfvén stationarity test in Fig. 5.With resolutions comparable to the one employed in Komissarov(2004), i.e., ∆ x ≈ . ≈ (cid:46) x (cid:46) .
2. Asmentioned in the previous tests, standing Alfvén waves seem tointroduce severe degradation of the order of convergence in MPmethods (we have also verified these results with MP5). This isvery likely related to the preservation of the D · B = ≤ x ≤ . B z is not uniform.In that region, the cutback of the electric displacement gener-ates numerical errors which accumulate mostly close to its lowerboundary (see the behavior of D x in − . (cid:46) x (cid:46) We perform a test (explored in extensive detail and high-resolution by Li et al. 2019) of the interaction between collid-ing Alfvén modes in suitably chosen 2D and 3D computationalboxes. In this section, we intend to reproduce the most basic re-sults of energy cascades from Alfvén wave interactions to showour GRFFE scheme’s ability to explore such phenomena in fur-ther detail in the future.On the respective numerical meshes, one initializes counter-propagating Gaussian 2D or 3D wave packets traveling along a s = = = U / U ( ) ( ) ( ) ( ) ( ) ( ) / τ [ Number of collisions ] U / U dU / dt ≃ - ⨯ - U / τ / Fig. 6.
Free energy U (normalized to its initial value U ) during the col-lision of Alfvén wave packets on numerical meshes (2D / ff erent line styles). Top:
Evolution of U / U for the set of 2D models. Slopes for the asymptotic linear relation be-tween U / U and ln t are indicated by dashed / dotted lines, comparingto Fig. 7 of Li et al. (2019). Bottom:
Free-energy evolution normalizedto U comparing to Figs. 2 and 5 from Li et al. (2019). The asymptoticslope for 3D models found by Li et al. (2019) is indicated by a graydashed line for reference. uniform guide field B y = B . Periodic boundary conditions fa-cilitate the recurring superpositions and interaction of the wavepackets, eventually triggering an energy cascade of rapid dissi-pation. The 3D Gaussian wave packets are initialized as B = B ˆy + B ∇ × ( φ ˆy ) , (75)where ˆy is the unit vector in the y -direction and the scalar field φ ( r ) = ξ l (cid:88) i = , exp (cid:32) − | r − r i | l (cid:33) . (76)In this section, ξ denotes the perturbation strength, l the width ofthe wave packet, with centers are located at r and r . We followLi et al. (2019) in choosing ξ = . l = . r = (0 . , . , . r = (0 . , . , .
5) for the 3D wave packets. With thissetup, the field perturbation is purely azimuthal with respect tothe y -axis. On a reduced 2D mesh, we initialize Gaussian wavepackets with magnetic fields B = B ˆy + B z ˆz , (77)where ˆz is the unit vector in the z -direction and B z = B ξ (cid:88) i = , exp (cid:32) − | r − r i | l (cid:33) . (78) Article number, page 10 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing t = τ t = τ t = τ t = τ t = τ - - P ( k ⊥ ) k ⊥ - ( ) - - k ⊥ P ( k ⊥ ) k ⊥ - Fig. 7.
Spectrum evolution for the 2D simulation of Alfvén interactionsinitialized according to Eq. (78). We show the spectral energy distribu-tion for wavenumbers k ⊥ (perpendicular to the guide field) in analogyto Fig. 6, Li et al. (2019). Top:
Spectral energy distribution at di ff erenttimes for the resolution 512 . Bottom:
Spectral energy distribution forselected times and wavenumbers (color code as in top panel) and dif-ferent resolutions, indicated by dashed (128 ), dotted (256 ) and solid(512 ) lines. No visible convergence is reached. We employ ξ = . l = . r = (0 . , .
25) and r = (0 . , . D × B / B , that results form an initial electric field D = ± ˆy × B , (79)with opposite signs for each Gaussian wave packet.After the initialization of the electromagnetic fields, thebounding box of length L = τ ( f cfl = . τ = . t /τ the number ofcollisions. For these tests, we employ MP7 reconstruction (witha fourth-order discretization of j (cid:107) ). Following Li et al. (2019),we employ the free energy U as the measure of total electromag-netic energy of the system e tot under removal of the backgroundmagnetic field B : U = e tot − (cid:90) d V B . (80)Fig. 6 shows the free energy for the collision of the wavepackets defined in Eqs. (76) and (78). The wave packets arespherical and - due to their curvature - prone to redistribute t = τ t = τ t = τ t = τ t = τ - - - - - - P ( k ⊥ ) k ⊥ - ( ) - - - - - - k ⊥ P ( k ⊥ ) k ⊥ - Fig. 8.
Spectrum evolution for the 3D simulation of Alfvén interactionsinitialized according to Eq. (75). We show the spectral energy distribu-tion for wavenumbers k ⊥ (perpendicular to the guide field) in analogyto Fig. 2, Li et al. (2019). Top:
Spectral energy distribution at di ff erenttimes for the resolution 384 . Bottom:
Spectral energy distribution forselected times and wavenumbers (color code as in top panel) and dif-ferent resolutions, indicated by dashed (128 ), dotted (256 ) and solid(384 ) lines. For wavenumbers k ⊥ (cid:46)
40, convergence is reached for theshown high resolution cases. energy across wave modes and rapid dissipation in cascade-like processes (e.g., Howes & Nielson 2013; Nielson et al.2013). Such processes are likely to be found along curved guide-fields, for example in magnetar magnetospheres. Collisions ex-cite waves of higher frequency than initially setup and eventuallytrigger the rapid decay of the wave free-energy. In order to makea more quantitative comparison of the results, we also computethe spectral distribution of free energy according to componentsof the propagation wave vector which are parallel ( k || ) and per-pendicular ( k ⊥ ) to the guide field, following the same prescrip-tion as in Li et al. (2019) (see Figs. 8 and 7). We compare our 2Dand 3D results with the reference work of Li et al. (2019) below. Our GRFFE code is able to reproduce the dissipation patternsof free electromagnetic energy presented in Fig. 5 of Li et al.(2019) for the 2D setup of Eq. (78). We note that our tests corre-spond to the lowest three mesh resolutions employed by Li et al.(2019). In the top panel of Fig. 6, we display the evolution ofthe inverted free energy ( U / U ) along with the slopes for their Article number, page 11 of 19 & A proofs: manuscript no. A_code_paper decay in 2D. Contrasting the findings by Li et al. (2019), the de-cay of U initially proceeds at the same rate ( s ≈ .
35) for all ofthe analyzed 2D models, independent of the chosen resolution.Only at later times, the slopes deviate and (roughly) approach thenumerical values given in Fig. 7 of Li et al. (2019). The redistri-bution of spectral energy happens at all times from the smaller tothe larger values of k ⊥ (Fig. 7), and an approximate k − ⊥ spectraldependence is observed for intermediate values 10 (cid:46) k ⊥ (cid:46) k ⊥ , k ⊥ , max , is resolution dependent (thefiner the resolution, the larger k ⊥ , max ). Hence, there is no evi-dence of spectral convergence in 2D (in agreement with Li et al.2019, see Fig. 7 lower panel). For the respective resolution, thedecay of U begins later than in Li et al. (2019), suggesting thatthe numerical di ff usivity in our method (combining MP7 recon-struction, a fourth-order accurate Runge-Kutta time-integratorand no additional driving terms in j || ) is smaller (compared toa fifth-order spatial WENO reconstruction, a third order accu-rate Runge-Kutta time-integrator, and an extra dissipation termin j || ). As a result, our models with a grid of 512 zones displayan evolution of U trending roughly in between the curves cor-responding to ∼ and 2048 in Li et al. (2019). A morethorough characterization of the (numerical) dissipation of ouralgorithm is considered in Paper II. For tests in 3D, we are limited to the two lower resolutions of thecorresponding model in Li et al. (2019) to stay within the com-putational costs that are reasonable for a test setup. In spite of thereduced resolution we employ compared to the literature, we finda remarkable agreement with the reference results. For instance,we observe a faster onset of the energy cascade than in 2D, andthe same asymptotic value of the free energy ( U / U ≈ . t onset after which dissi-pation commences in 3D models. After this time there is a rela-tively sharp drop of the free energy, which tends to level o ff forsu ffi ciently long times. However, the feature that unanimouslysets t onset (according to the definition of Li et al. 2019) is found inthe evolution of the spectral energy distribution. Before t onset , thecolliding Gaussian packets shu ffl e energy towards smaller scales(larger values of the wave number k ) until a maximum k = k max is reached. However, after t = t onset there exists a redistributionof energy from the smaller to the larger scales, which manifestitself as an increasing spectral power at small values of k . For Liet al. (2019), t onset ≈ τ . The conclusions of Li et al. (2019) arebased upon 3D models with finer resolution than the ones em-ployed in this section (e.g., models with 512 and 768 zones).With a smaller resolution, we also find a time after which there isa steep drop o ff the free energy, which begins at t (cid:38) τ (Fig. 8top panel). Nevertheless, we do not clearly see an increase in thespectral power at small values of k ⊥ , but we observe a decreaseof power for k ⊥ (cid:38)
30 for a time 30 τ (cid:46) t onset (cid:46) τ . At about thistime, the fast drop o ff in U takes place. Since we have evolvedour models longer in time, we note that at t = τ there is adecrease of spectral energy at intermediate values 10 (cid:46) k ⊥ (cid:46) ff erence that we observe some (small) excess of powerin the range 10 (cid:46) k || (cid:46)
50. As in 2D, most of the (small) quantita-
Stability ( ) Stability ( ) Relaxation ( ) Relaxation ( ) M a g n e t o s p h e r i c e n e r g y e / e d Cartesian ( ) Spherical ( ) - - - ( ms ) | e / e d | Cartesian ( ) ( ms ) Spherical ( ) Fig. 9.
Stability and relaxation test of magnetar magnetospheres en-dowed with an analytic dipole field structure for di ff erent resolutions.Upper panels show the total magnetospheric energy normalized to theenergy of the corresponding dipole. The lower panels show the timederivative of the energy normalized to the dipole. 16 and 32 pointsper stellar radius in a 3D Cartesian C arpet grid ( left ) correspond tothe setup in Mahlmann et al. (2019). Axially symmetric simulationsin spherical coordinates ( right ) are set up with the indicated resolutionat the stellar surface. Black (solid and dashed) lines in the upper-rightpanel correspond to a simulation on a smaller domain, extending up to r = . M (cid:12) − ∆ r , in which the pulses of the initial relaxation havesu ffi cient time to leave the domain. tive di ff erences observed in the comparison with Li et al. (2019)results can be attributed to the di ff erent order of accuracy of ourcodes and the di ff erent terms entering in the current parallel tothe magnetic field.
5. Astrophysically Motivated Tests
The magnetospheres of magnetars are a well-suited laboratoryfor numerical methods dealing with force-free plasma, and wehave explored their dynamics in Mahlmann et al. (2019). Priorto those numerical simulations of a potentially very dynamicscenario, we performed numerical tests to assess the ability ofour GRFFE code to maintain the structural stability of a mag-netosphere around a spherical, non-rotating neutron star with adipolar magnetic field. We use a spherical mask to cut out theneutron star interior in order to avoid dealing with the equation
Article number, page 12 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing of state of nuclear matter, the di ff erent phases of matter whichmay occur inside of the neutron star, and the solid structure ofthe stellar crust. This is achieved by setting an internal bound-ary in a 3D Cartesian grid (i.e., stair-stepping along the sphericalboundary mask) inside of which the evolution is frozen (see be-low). We note that 3D Cartesian coordinates are neither adaptedto the spherical shape of the neutron star nor the axial symme-try of the magnetospheric dipole. Therefore, we expect signifi-cantly improved results when employing GRFFE in a sphericalcoordinate system. We compare the results of simulations in 3DCartesian coordinates to axisymmetric (2D) tests in spherical co-ordinates, with the magnetic axis aligned to the symmetry axisof the initial data.It is straightforward to specify the employed initial data inspherical coordinates ( r , θ, φ ) and, subsequently, map it to thecomputational grid. The analytically derived equilibrium dipolarmagnetic field in (the coordinate basis of) spherical coordinatesreads: B = (cid:32) θ r , sin θ r , (cid:33) , D = (0 , , . (81)The Cartesian simulations are conducted in a 3D box with di-mensions [4741 . M (cid:12) × . M (cid:12) × . M (cid:12) ] with a gridspacing of ∆ x ,y, z = . M (cid:12) on the coarsest grid level. Forthe chosen magnetar model of radius R ∗ = . M (cid:12) this cor-responds to a [512 R ∗ × R ∗ × R ∗ ] box with a grid spac-ing of ∆ x ,y, z = R ∗ . For the low and high-resolution tests weemploy seven and eight additional levels of mesh refinement,each increasing the resolution by a factor of two and encompass-ing the central object, respectively. This means that the finestresolution of our models (close to the magnetar surface) are ∆ min x ,y, z = . × R ∗ = . M (cid:12) and ∆ min x ,y, z = . × R ∗ = . M (cid:12) for the low and high-resolution models. This corre-sponds to 16 and 32 points per R ∗ , respectively. The sphericalsimulations are conducted in axial symmetry enclosing the vol-ume [9 . M (cid:12) , . M (cid:12) − ∆ r ] × [0 , π ]. In order to issue aresolution which is comparable to the Cartesian setup, we em-ploy ∆ r ∈ [ R ∗ / , R ∗ /
32] and ∆ θ ∈ [ π/ , π/ t = . M (cid:12) (cid:39) .
84 ms. We provideextensive details on the internal boundary conditions (frozenelectromagnetic fields but balanced radial current) in Mahlmannet al. (2019). For this section we employ f CFL = .
2, the MP7reconstruction and a fourth-order accurate discretization of j || .The stability test initializes the dipole structure throughoutthe entire computational domain and tracks the stability duringa dynamical evolution. The relaxation test is even more chal-lenging than the stability test since it requires the time evolu-tion toward the physical topology set by the boundary condi-tions. Precisely, in a relaxation test we fix the dipolar structureinside of the star, but fill the magnetosphere with a purely radialfield at the start of the simulation. Once initialized, the energyof the dipole magnetosphere (cf. Mahlmann et al. 2019) is wellconserved (stability) or else gradually approaches the dipole en-ergy (relaxation) once all initially introduced perturbations leavethe domain. Figure 9 shows stability and relaxation tests of thedipole magnetosphere for di ff erent resolutions in both, Cartesian(3D) and spherical (2D, axial symmetry) coordinates.The initial spike of the relaxation model can be attributed to asurge of electromagnetic energy during a rapid rearrangement inthe early phase. The excited energy pulses propagate as plasmawaves through the magnetosphere. A part of these pulses is con-fined to closed field lines in the vicinity of the central object. The rest of this energy propagates outwards through the domain. Asthe dissipation of electromagnetic energy in collisions of force-free waves strongly depends on the employed resolution (seeSect. 4.2, and Paper II) and grid geometry, the confined energypulses remain within the domain longer for higher resolution andspherical coordinates. This is visible best in the relaxation testpresented in Fig. 9. For di ff erent resolutions, the asymptotic en-ergy di ff ers by < ff ects, associated to the total simulated time, canbe partly addressed by considering a computational domain witha reduced outer radial boundary. In the top right panel of Fig. 9,we show the time evolution in a reduced computational domainwith black lines. The abrupt drop-o ff the magnetospheric energyis due to the desertion of the initial perturbation through the outerradial boundary. Remarkably, the energy level to which thesemodels evolve is the same as the corresponding stability testswith the corresponding numerical resolution. In this section, we explore the full 3D capabilities of our newlydeveloped code in spherical coordinates ( r , θ, φ ) by consideringthe stability test from the previous section, but tilting the mag-netic dipole axis by an angle α with respect to the spherical po-lar axis along θ =
0. For this, we carry out the transformation B t = R − α B ( ˜r ), where ˜r = R α r , B corresponds to the initial datagiven in Eq. (81), and R − α = cos α sin θ − sin α cos θ sin φχ sin α cos φ − sin α csc θ cos φχ cos α − sin α cot θ sin φ ,χ = (cid:113) sin θ cos φ + (sin α cos θ − cos α sin θ sin φ ) . (82)We choose α = ◦ for simulations that are exclu-sively conducted in a spherical domain with dimensions[9 . M (cid:12) , . M (cid:12) − ∆ r ] × [0 , π ] × [0 , π ]. In order to use aresolution which is comparable to Sect. 5.1.1, we employ ∆ r ∈ [ R ∗ / , R ∗ / , R ∗ /
32] and ∆ φ = ∆ θ ∈ [ π/ , π/ , π/ t = M (cid:12) (cid:39) .
48 ms with f CFL = .
25, using MP7 spatial reconstruction, and the defaultfourth-order discretization of j || .Besides the measurement of the total energy in the magne-tosphere, in order to quantify the deviation of the numerical so-lution B from the analytical (initial) configuration B , we definethe l error norm of the magnetic field as (cid:107) (cid:15) (cid:107) = N (cid:115)(cid:88) i , j , k (cid:104) B ( r i , θ j , φ k ) − B ( r i , θ j , φ k ) (cid:105) , (83)where indices i , j , and k extend over all the computational cells.Fig. 10 shows the evolution in time of the error norm in Eq. (83)normalized to the magnetic field strength at the pole of the neu-tron star, B p . Increasing the mesh resolution by a factor of 2 (i.e.from ∆ r = R ∗ /
16 to R ∗ /
32) reduces the error by roughly thesame factor. At the same time, the total magnetospheric energyslightly increases throughout the simulation time ( (cid:46)
Article number, page 13 of 19 & A proofs: manuscript no. A_code_paper × - × - × - × - ϵ / B p ( ms ) T o t a l e n e r g y e / e d Fig. 10.
Stability test of magnetar magnetospheres endowed with a tilted (30 degrees) analytic dipole field structure for di ff erent resolutions. Left:
Evolution of the l -norm of the error with respect to the analytical (initial) configuration normalized to the magnetic field strength B p atthe magnetar pole ( top ), and of the total magnetospheric energy content ( bottom ). Right:
3D impression of the final state ( t ≈ . ∆ r = R ∗ / global structure of the 3D tilted dipole is conserved throughoutthe simulation (Fig. 10), with only slight kinks arising aroundthe polar axis of the spherical coordinates.Solving hyperbolic PDEs in spherical coordinates su ff erfrom very small timesteps due to the fact that cell volumes arenot constant in space and get smaller as the polar axis or theorigin of the coordinate system are approached. The timestepis proportional to r × sin( ∆ θ/ × ∆ φ and hence becomes pro-hibitively small for full 3D simulations with high angular reso-lutions, this is the reason for the shorter simulation time com-pared to the axially symmetric models in Sect. 5.1.1). In order tomitigate this limitation for simulations in spherical coordinates,additional grid coarsening or filtering approaches will be nec-essary when considering computationally feasible long-term 3Devolution (cf. Obergaulinger & Aloy 2017; Mewes et al. 2020;Zlochower et. al. 2020; Obergaulinger & Aloy 2020). The im-pact of the tilt across the axis on the magnetospheric energy (asobserved in Fig. 10) may be further diminished by such tech-niques.In this test we employ the full 3D capacities of our GRFFEmethod in spherical coordinates. The tilted magnetar magneto-spheres maintain their topology stable for ∼
32 light-crossingtimes of the central object. Extrapolating the deviation of themagnetospheric energy to longer times is uncertain due to thenon-monotonic evolution. However, we foresee that stationarytilted magnetospheres may be maintained approximately stablesu ffi ciently for more than a few hundred light-crossing times ofthe central object. These longer periods of evolution may su ffi ceto address numerically dynamical phenomenae in the magneto-sphere (e.g. Carrasco et al. 2019; Mahlmann et al. 2019). Be-sides, our results show that increasing the resolution decreasesthe global error, hence, if needed, finer grids may be used to ad-dress longer evolutionary times. Blandford & Znajek (1977) presented analytic equilibrium solu-tions of BH magnetospheres by applying perturbation techniquesto the Grad-Shafranov equation (GSE) which match the Znajekcondition (Znajek 1977) at the BH horizon and the flat spacesolution of Michel (1973) at infinity. One of these results is amonopole-like magnetic field, which is often adapted to the so-called split monopole by mirroring the field quantities across theequatorial plane. The latter is a necessary step to avoid diver-gences of the magnetic field. In this section, however, we followKomissarov (2004) in considering the monopole field structureto avoid the challenge of resolving a current sheet at the equator.The monopole electromagnetic fields for slowly spinning BHs( a ∗ (cid:28)
1) as derived in Blandford & Znajek (1977) can be writtenin the spatial components of vectors in Boyer-Lindquist coordi-nates ( r , θ, φ ) as follows: B = B (cid:32) − sin θ √ γ , , − a ∗ sin θ αg φφ (cid:33) , D = B (cid:32) , − Ω F + β φ αg θθ sin θ, (cid:33) . (84)Here, Ω F is the field line angular velocity as defined for axi-ally symmetric equilibrium solutions, and we employ B = M (cid:12) × M (cid:12) × M (cid:12) ] with a grid spacing of ∆ x ,y, z = M (cid:12) on the coarsest grid level. We employ eight ad-ditional levels of mesh refinement, each increasing the resolu-tion by a factor of two and encompassing the central object, re-spectively. This means that the finest resolution of our Carte-sian models is ∆ min x ,y, z = . M (cid:12) . The spherical simulationsare conducted on a 2D slab (axial symmetry) with extensions Article number, page 14 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing ( M ⊙ ) M a g n e t o s p h e r i c e n e r g y e / e o Spherical ( ) Cartesian ( )
40 60 80 100 120Time ( M ⊙ ) ( r g ) Ω F / Ω B H Cartesian ( ) ( r g ) Spherical ( ) Fig. 11.
Time evolution of the Schwarzschild monopole ( Ω F = Ω BH /
2) of a slowly spinning Kerr BH ( a ∗ = . M =
1) in Cartesian (3D C arpet gridwith nine refinement levels, with the highest resolution of 0 . M (cid:12) completely enclosing the central object) and spherical (2D, axsymmetric)coordinates. The spacetime metric is dynamically evolving. Left:
Evolution of the total magnetospheric (electromagnetic) energy normalized tothe initial value, e . Right:
Field line angular velocity along the equatorial plane. The final value is shown in a strong color. Intermediate statesthroughout the simulation are depicted by lighter colored lines (the strength of the color increasing with simulation time).
Fig. 12.
Simulations of Wald magnetospheres for a ∗ = . M = t = M (cid:12) . Left:
3D Cartesian C arpet grid (vacuum spacetime modeledby M ac L achlan ). Right:
2D (axially symmetric) spherical grid (vacuumspacetime modeled by S pherical
BSSN). The poloidal field is indicatedby streamlines, the toroidal field by red and blue colors (color scalecoincides for all panels) indicating whether the toroidal field leaves orenters into the displayed plane, respectively. The BH ergosphere is de-noted by a solid white line. [0 , M (cid:12) ] × [0 , π ]. In order to use a resolution which is com-parable to the Cartesian setup, we employ ∆ r = .
032 and ∆ θ = π/
64. The setup is evolved for a period of t = M (cid:12) . Fig. 13.
Simulations of a Wald magnetosphere for a ∗ = . t = M (cid:12) . Left:
Same format and color scale as inFig. 12.
Right:
3D magnetic field line impression of the Wald magne-tosphere. Darker colors indicate a stronger twist of the magnetic field.
For this section we employ f CFL = .
25, the MP7 reconstructionand a fourth-order accurate discretization of j || .Fig. 11 summarizes the time evolution of the monopole fieldfor a dynamically evolving spacetime metric for both Cartesian(3D) and spherical (2D, axisymmetric) meshes. As the magne-tospheres considered in this section (and the following one) are Article number, page 15 of 19 & A proofs: manuscript no. A_code_paper a * = a * = a * = Ω F / Ω B H r = - - θ T o r o i d a l m a g n e t i c f i e l d B ϕ r = r H Fig. 14.
One-dimensional values of the field line angular velocity ( top )and the toroidal magnetic field ( bottom ) for the Wald test (CartesianC arpet grid) at t = M (cid:12) and using di ff erent BH dimensionless ra-pidities (see legends). The interpolation radius (for the extraction in aCartesian grid) is indicated in the respective panel, corresponding to theergosphere radius at the equator ( top ) or the BH horizon radius ( bot-tom ). idealized cases, namely a magnetic monopole and with unboundenergy, we do not couple the field energy to the source terms ofthe BSSN equations. During a transient phase in which the met-ric terms relax to the chosen mesh and gauge, the electromag-netic fields can di ff er significantly from their initial state. Thistest demonstrates that, while the spacetime evolves, the electro-magnetic fields relax towards the equilibrium given by Eq. (84)concurrently. Though the energy evolution shown in Fig. 11 ap-proaches the energy of the initial model rather well, the resultingequilibrium has to be taken with care. The evolution of dynam-ical spacetimes and corresponding GRFFE fields can be subjectto the influence of small changes of the BH mass and spin (dueto finite numerical resolution), as well as an involved array ofgeometric source terms (see Sect. 2.2).The geometric (i.e., spacetime) quantities determined by theinitial data for spinning BHs presented in Liu et al. (2009) re-lax to their equilibrium state depending on the chosen numericalresolution of the mesh and specification of gauge quantities (i.e.,the lapse and shift) during an initialization phase. The choice ofthese quantities is preferably done in a way which causes theleast possible noise across all metric quantities during their evo-lution. As an example, instead of providing the spacetime datawith the analytic lapse function defined in Boyer-Lindquist co-ordinates (Liu et al. 2019), we specify the lapse initially as:˜ α (0) = × (cid:34) + (cid:18) + M r (cid:19) (cid:35) − . (85)With this initialization, the spacetime relaxes swiftly to its equi-librium state during the first ∆ t init ≈ M (cid:12) . The tests presentedin this section give some important hints on the strategies chosento set up BH magnetospheres for our future research. The goal ofthis test was to show that the magnetospheric data is conserved throughout the (dynamic) relaxation of the spacetime induced,e.g., by the BSSN algorithms of the E instein T oolkit . As boththe magnetospheric energy as well as the field line angular ve-locity at the equator, are recovered after ∆ t init , our GRFFE codepasses this test of spacetime-field coupling. The immersion of a BH into a magnetic field which is uniformat infinity was originally suggested by Wald (1974) and then ex-plored throughout the literature, both as a test and as a laboratoryfor force-free plasma (Komissarov 2004; Komissarov & McKin-ney 2007; Carrasco & Reula 2017; Parfrey et al. 2019). In thissection, we reproduce the initial data of the Wald magnetosphereof a Schwarzschild BH in Boyer-Lindquist coordinates (rescaledaccording to the prescription of Liu et al. 2009) and evolve it fordi ff erent BH spins. We therefore, extend the testing of the GRcapacities of our code to rapidly spinning BH (up to a ∗ = . r , θ, φ ) can be ini-tialized as follows: B = B (cid:32) − (cid:114) r + r cos θ, θ √ r + (2 + r ) , (cid:33) , D = (0 , , . (86)We employ B = M (cid:12) × M (cid:12) × M (cid:12) ] with a gridspacing of ∆ x ,y, z = M (cid:12) on the coarsest grid level. We employnine additional levels of mesh refinement, such that the finestresolution of our Cartesian models is ∆ min x ,y, z = . M (cid:12) . Thespherical simulations are conducted in axial symmetry with ex-tensions [0 , M (cid:12) ] × [0 , π ]. In order to use a resolution whichis comparable to the Cartesian setup, we employ ∆ r = .
032 and ∆ θ = π/
64. The setup is evolved for a period of t = M (cid:12) inCartesian coordinates and t = M (cid:12) in spherical coordinates.For this section we employ f CFL = .
25, the MP7 reconstructionand a fourth-order accurate discretization of j || .Figure 12 shows the results from time evolution of thesefields in spacetimes of rotating BHs for a selected case ( a ∗ = . a ∗ = .
9, Fig. 13), preventingthe development of static magnetospheric conditions (cf. Komis-sarov 2004). The overall topology of the magnetic field through-out the BH ergosphere broadly coincides with respective equi-librium solutions of Kerr magnetospheres (as derived, e.g., inNathanail & Contopoulos 2014; Mahlmann et al. 2018).The simulations in spherical coordinates are significantlymore expensive than in Cartesian coordinates, due to the severerestrictions on the timestep imposed by the converging sphericalmesh close to the central singularity. Future code developmentswill include mesh-coarsening strategies close to r = a ∗ (cid:38) . θ , which makes thisnumerical experiment (currently) too expensive as a validationtest of our code. Even in axisymmetric simulations, the timesteprestriction from the θ coordinate, which imposes a timestep pro-portional to r × ∆ θ , is too restrictive when the coordinate ori- Article number, page 16 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing gin is included in the computational domain, which is necessaryfor a spacetime evolution without excision as the one employedhere. In order to alleviate this shortcoming of doing high resolu-tion simulation in spherical coordinates that include the origin,algorithms to circumvent the timestep restrictions imposed byboth the θ and φ coordinates are currently being developed (Zlo-chower et. al. 2020).In Fig. 14 we extract the field line angular velocity andtoroidal magnetic field at di ff erent locations for Cartesian coor-dinates for comparison with Fig. 5 in Komissarov (2004) (com-puted for a BH with a ∗ = . ff erent from the literature in order to represent thecomplete range of BH spins. We find that our GRFFE code re-produces qualitatively the results in the literature, though somedi ff erences remain to be mentioned. Komissarov (2004) usesspherical coordinates and axial symmetry, as opposed to our 3Dsimulations with mesh refinement. More even, the angular reso-lution of 800 cells in the θ -direction is almost ten times the res-olution which we have used on our finest refinement level (theresolution limit is simply imposed by the aim of running numer-ical tests that do not consume disproportionate computationalresources). Without evolving the space time (as in Komissarov2004) the numerical grid may extend outwards in the radial di-rection from the event horizon as a boundary, in practice, excis-ing the central singularity and allowing for significantly largertimesteps. The quantitative di ff erence in the shape of the angularvelocity distribution (V-shape in Fig. 5 of Komissarov 2004, vs.U-shape in Fig. 14) may, hence, be significantly improved withincreasing resolution or by resorting to a GRFFE code in spher-ical coordinates (as we plan to do once the timestep restrictionsare overcome). Also, we point out that we show the toroidal com-ponent of the magnetic field B rather than H . The overall formof the toroidal field for the rapidly rotating case ( a ∗ = .
9) cor-responds well (up to a di ff erence in sign) with Fig. 5 of Komis-sarov (2004). The BH magnetosphere simulations in this sectionhave been repeated on a fixed Kerr-Schild background metric,e ff ectively confirming the results presented in this section. Fora direct comparison of numerical data, comparable resolutionsand exact convergence (longer simulations) are required; this isbeyond the scope of the test presented here.In conclusion, especially field lines threading the ergosphereare gradually twisted by the rotating BH. Due to the broad coin-cidence with Komissarov (2004), and the reproduction of mag-netospheres which resemble respective equilibrium solutions ofthe (Nathanail & Contopoulos 2014; Mahlmann et al. 2018), theWald magnetosphere test is passed.
6. Conclusions
We have developed a new GRFFE code that models magneticallydominated plasma in dynamical spacetimes with support forboth Cartesian and spherical coodinates provided by the C arpet grid of the E instein T oolkit . Our simulation tool combines tech-niques from an array of literature (especially Komissarov 2002;McKinney & Gammie 2004; Palenzuela et al. 2009; Paschalidis& Shapiro 2013; Parfrey et al. 2017) and improves further onnumerical strategies as well as the understanding of their limits: – We explicitly couple the continuity equation of charge to ourconservative scheme (Sect. 2) and, thus, ensure a consistentmodeling of the force-free current (Eq. 48). – We employ a hyperbolic / parabolic cleaning of errors (ex-tending, e.g., the techniques in Palenzuela et al. 2009;Mignone & Tzeferacos 2010, to general relativity). Allowing for arbitrary advection speeds for the cleaning of divergenceerrors significantly improves the conservation of total chargein spacetimes of spinning BHs (see Appendix A). – The current parallel to the magnetic field j (cid:107) is the dominantdriver of resistivity in GRFFE schemes. In case of the force-free current (Eq. 63), high-order discretization allows us tomodel (smooth) force-free plasma waves with nearly theo-retical order (Sect. 2.1, Paper II), and di ff using only due tonumerical resistivity. – Current sheets are genuine regions of significant physical re-sistivity. Conventional GRFFE methods (i.e., schemes thatdo not include phenomenological models of artificial phys-ical resistivity) are unable to properly resolve such resistivelayers, especially for highly accurate reconstruction methods(Paper II). Thus, at discontinuities, the order of convergenceof GRFFE is significantly reduced; a true limit of applicabil-ity of GRFFE is reached.Writing our E instein T oolkit thorn from scratch enabled usto implement suitable finite volume integrators for both Carte-sian and spherical coordinates. Spherical coordinate systemsprove exceptionally valuable for the highly accurate modelingof magnetar magnetospheres (Sect. 5.1). For the simulation ofBH magnetospheres on dynamically evolving spacetimes, ourGRFFE method will benefit from future updates to the supportof spherical coordinates in the E instein T oolkit (Mewes et al.2018, 2020). With the presented numerical code we broadly ex-ploit the modular nature of the E instein T oolkit and implementa cutting-edge tool for the simulation of GRFFE on dynamicalspacetimes in Cartesian and spherical coordinates.
7. Acknowledgments
JM acknowledges a Ph.D. grant of the
Studienstiftungdes Deutschen Volkes . We acknowledge the support fromthe grants AYA2015-66899-C2-1-P, PGC2018-095984-B-I00,PROMETEO-II-2014-069, and PROMETEU / / ff ort of the U.S. De-partment of Energy (DOE) O ffi ce of Science and the NationalNuclear Security Administration. Work at Oak Ridge NationalLaboratory is supported under contract DE-AC05-00OR22725with the U.S. Department of Energy. VM also acknowledges par-tial support from the National Science Foundation (NSF) fromGrant Nos. OAC-1550436, AST-1516150, PHY-1607520, PHY-1305730, PHY-1707946, and PHY-1726215 to Rochester In-stitute of Technology (RIT). The shown numerical simulationshave been conducted on infrastructure of the Red Española deSupercomputación (AECT-2020-1-0014) as well as of the
Uni-versity of Valencia
Tirant and LluisVives supercomputers.
References
Alcubierre, M. 2008, Introduction to 3 + Article number, page 17 of 19 & A proofs: manuscript no. A_code_paper
Baumgarte, T. W. & Shapiro, S. L. 2010, Numerical Relativity: Solving Ein-stein’s Equations on the Computer (Cambridge University Press)Beloborodov, A. M. & Thompson, C. 2007, ApJ, 657, 967Beskin, V. S. 1997, Sov. Phys.-Usp., 40, 659Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433Bona, C., Masso, J., Seidel, E., & Walker, P. 1998, arXiv e-prints, grBrown, D., Diener, P., Sarbach, O., Schnetter, E., & Tiglio, M. 2009, PRD, 79,044023Brown, J. D. 2005, PRD, 71, 104011Brown, J. D. 2009, PRD, 79, 104029Carrasco, F., Viganò, D., Palenzuela, C., & Pons, J. A. 2019, Monthly Notices ofthe Royal Astronomical Society: Letters, 484, L124Carrasco, F. L. & Reula, O. A. 2016, PRD, 93, 085013Carrasco, F. L. & Reula, O. A. 2017, PRD, 96, 063006Cerdá-Durán, P., Font, J. A., Antón, L., & Müller, E. 2008, A&A, 492, 937Collins, D. C., Xu, H., Norman, M. L., Li, H., & Li, S. 2010, ApJSS, 186, 308Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351Cordero-Carrión, I., Ibáñez, J. M., Gourgoulhon, E., Jaramillo, J. L., & Novak,J. 2008, PRD, 77, 084007Darmois, G. 1927, Mémorial des sciences mathématiques, 25, 1Dedner, A., Kemm, F., Kröner, D., et al. 2002, JCP, 175, 645Diener, P., Dorband, E. N., Schnetter, E., & Tiglio, M. 2007, Journal of ScientificComputing, 32, 109Dreyer, O., Krishnan, B., Shoemaker, D., & Schnetter, E. 2003, PRD, 67, 024018Etienne, Z. B., Wan, M.-B., Babiuc, M. C., McWilliams, S. T., & Choudhary, A.2017, CQG, 34, 215001Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a,ApJ, 875, L1Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019b,ApJ, 875, L5Fourès-Bruhat, Y. 1952, Acta Mathematica, 88, 141, english translation at http: // hdl.handle.net / / + ff erence methods for ordinary and partial di ff erentialequations: Steady-state and time-dependent problems, Other Titles in AppliedMathematics (Society for Industrial and Applied Mathematics)Li, X. & Beloborodov, A. M. 2015, ApJ, 815, 25Li, X., Zrake, J., & Beloborodov, A. M. 2019, ApJ, 881, 13Lichnerowicz, A. 1944, J. Math. Pures et Appl., 23, 37Liu, H., Zong, Q. G., Zhang, H., et al. 2019, Nature Communications, 10, 1040Liu, Y. T., Etienne, Z. B., & Shapiro, S. L. 2009, PRD, 80Lö ffl er, F., Faber, J., Bentivegna, E., et al. 2012, CQG, 29, 115001Lyutikov, M. 2003, MNRAS, 346, 540Lyutikov, M. 2009, Monthly Notices of the Royal Astronomical Society, 396,1545MacDonald, D. & Thorne, K. S. 1982, MNRAS, 198, 345Mahlmann, J. F., Akgün, T., Pons, J. A., Aloy, M. A., & Cerdá-Durán, P. 2019,MNRAS, 490, 4858Mahlmann, J. F., Aloy, M. A., Mewes, V., & Cerdá-Durán, P. 2020a, arXiv e-prints, arXiv:2007.06599Mahlmann, J. F., Cerdá-Durán, P., & Aloy, M. A. 2018, MNRAS, 477, 3927Mahlmann, J. F., Levinson, A., & Aloy, M. A. 2020b, MNRASMartí, J. M. & Müller, E. 2003, LRRel, 6, 7McKinney, J. C. 2006, MNRAS, 367, 1797McKinney, J. C. & Gammie, C. F. 2004, ApJ, 611, 977Mewes, V., Zlochower, Y., Campanelli, M., et al. 2020, PRD, 101, 104007Mewes, V., Zlochower, Y., Campanelli, M., et al. 2018, PRD, 97, 084059Michel, F. C. 1973, ApJ, 180, 207Michel, F. C. 1973, ApJ, 180, L133Mignone, A. 2014, JCP, 270, 784Mignone, A. & Tzeferacos, P. 2010, JCP, 229, 2117Miranda-Aranguren, S., Aloy, M. A., & Rembiasz, T. 2018, MNRAS, 476, 3837Montero, P. J., Baumgarte, T. W., & Müller, E. 2014, PRD, 89, 084043Montero, P. J. & Cordero-Carrión, I. 2012, PRD, 85, 124037 Munz, C., Schneider, R., Sonnendrucker, E., & Voss, U. 1999, Academie desSciences Paris Comptes Rendus Serie Sciences Mathematiques, 328, 431Nathanail, A. & Contopoulos, I. 2014, ApJ, 788, 186Nielson, K. D., Howes, G. G., & Dorland, W. 2013, PPl, 20, 072303Obergaulinger, M. & Aloy, M. Á. 2017, MNRAS, 469, L43Obergaulinger, M. & Aloy, M. Á. 2020, MNRAS, 492, 4613Palenzuela, C., Garrett, T., Lehner, L., & Liebling, S. L. 2010, PRD, 82, 044045Palenzuela, C., Lehner, L., Reula, O., & Rezzolla, L. 2009, MNRAS, 394, 1727Parfrey, K., Giannios, D., & Beloborodov, A. M. 2015, MNRAS, 446, L61Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101Parfrey, K., Spitkovsky, A., & Beloborodov, A. M. 2017, MNRAS, 469, 3656Paschalidis, V. & Shapiro, S. L. 2013, PRD, 88, 104031Punsly, B. 2001, Black hole gravitohydromagnetics (Springer)Punsly, B. 2003, ApJ, 583, 842Rezzolla, L. & Zanotti, O. 2013, Relativistic hydrodynamics (OUP Oxford)Roache, P. J. 1997, ARFM, 29, 123Scharlemann, E. T. & Wagoner, R. V. 1973, ApJ, 182, 951Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Classical and Quantum Gravity,21, 1465Shibata, M. 2015, Numerical relativity, 100 Years of General Relativity (WorldScientific Publishing Company)Shibata, M. & Nakamura, T. 1995, PRD, 52, 5428Suresh, A. & Huynh, H. T. 1997, JCP, 136, 83Takahashi, M., Nitta, S., Tatematsu, Y., & Tomimatsu, A. 1990, ApJ, 363, 206Thompson, C. & Duncan, R. C. 1995, MNRAS, 275, 255Thornburg, J. 2004, CQG, 21, 743Thorne, K. S., Price, R. H., & MacDonald, D. A., eds. 1986, Black holes: Themembrane paradigm (Yale University Press)Tondeur, P. 2012, Geometry of foliations, Vol. 90 (Birkhäuser)Turolla, R., Zane, S., & Watts, A. L. 2015, Reports on Progress in Physics, 78,116901Uchida, T. 1997, PRE, 56, 2181van Leer, B. 1977, Journal of Computational Physics, 23, 276Wald, R. 2010, General Relativity (University of Chicago Press)Wald, R. M. 1974, PRD, 10, 1680York, J. W., J. 1979, in Sources of Gravitational Radiation, ed. L. L. Smarr, 83–126Yu, C. 2011, MNRAS, 411, 2461Zlochower et. al., Y. 2020, in preparationZnajek, R. L. 1977, MNRAS, 179, 457 Article number, page 18 of 19. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Multi-Coordinate Implementation and Testing κ Ψ = c h = κ Ψ = c h = κ Ψ = c h = κ Ψ = c h = - T o t a l d i v e r g e n c ee rr o r ϵ ∇ · B ( t ) × - × - × - × - Time ( M ⊙ ) Ψ m a x Fig. A.1.
Time evolution of numerical errors (div B = top panel) andthe corresponding maximum cleaning potential Ψ ( bottom panel). Wepresent combinations of di ff erent κ Ψ and c h for a fixed κ Φ = . Appendix A: Code Performance (3D Cartesian):Cleaning of Errors
We describe the implementation of the generalized Lagrangemultiplier method employed to preserve the electromagnetic dif-ferential constraints div B = D = ρ in Sect. 3.5. Thissection explores the code performance for the cleaning of diver-gence errors in the 3D Cartesian version of the BH monopoletest (with the setup from Sect. 5.2) for di ff erent choices of theparameters governing the numerical cleaning of errors. We mea-sure the numerical errors to the aforementioned conditions byconsidering the global measures: ε ∇· B ( t ) = (cid:90) [ ∇ · B ( t )] d V − (cid:90) [ ∇ · B ( t = V . (A.1)Here, we employ the 3D region outside of the BH horizon as anintegration region, and subtract the initially present errors due tothe discretization of the magnetic field. Figure A.1 shows theevolution of numerical errors and the corresponding cleaningpotentials for di ff erent combinations of the parameters κ Ψ , and c h . The optimization of these parameters may di ff er for di ff erentapplications and can be critical in highly dynamical processeswhere strong numerical violations of the divergence constraintsoccur (e.g., by strong violations of the force-free conditions, seealso the discussion in Mahlmann et al. 2019). For the tests athand, the exact calibration of the parameters of the cleaningmethod may have very small e ff ects (the total magnetosphericenergy presented in Fig. 11 is not notably changed by any of thedi ff erent combinations shown in Fig. A.1). However, their anal-ysis provides crucial information about the code’s performance,and in other applications the proper calibration of the cleaningroutines has a significant impact (Mahlmann et al. 2019).As is lucidly shown in Fig. A.1, the introduction of the su-perluminal advection velocity c h into the augmented system ofequations (Eq. 40) for divergence cleaning reduces the error ε ∇· B (especially in the early and late phase of the evolution) signif-icantly. Furthermore, the maximum magnitude of the cleaningpotential Ψ decreases by two orders of magnitude. Small varia-tions in ε ∇· B are also observed for stronger damping of errors bygreater values for κ Ψ . Though the presented tests for flat back-ground geometries employ c h =
1, we conclude from the re-sults in Fig. A.1 that c h = / parabolic cleaning of numerical errors should and willbe readjusted in the light of future applications of our GRFFEmethod.parabolic cleaning of numerical errors should and willbe readjusted in the light of future applications of our GRFFEmethod.