Computational General Relativistic Force-Free Electrodynamics: II. Characterization of Numerical Diffusivity
AAstronomy & Astrophysics manuscript no. B_code_paper c (cid:13)
ESO 2020July 16, 2020
Computational General Relativistic Force-Free Electrodynamics: II.Characterization of Numerical Diffusivity
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
Scientific codes are an indispensable link between theory and experiment; in (astro-)plasma physics, such numerical tools are onewindow into the universe’s most extreme flows of energy. The discretization of Maxwell’s equations - needed to make highly mag-netized (astro)physical plasma amenable to its numerical modeling - introduces numerical di ff usion. It acts as a source of dissipationindependent of the system’s physical constituents. Understanding the numerical di ff usion of scientific codes is the key to classify theirreliability. It gives specific limits in which the results of numerical experiments are physical. We aim at quantifying and characterizingthe numerical di ff usion properties of our recently developed numerical tool for the simulation of general relativistic force-free elec-trodynamics, by calibrating and comparing it with other strategies found in the literature. Our code correctly models smooth wavesof highly magnetized plasma. We evaluate the limits of general relativistic force-free electrodynamics in the context of current sheetsand tearing mode instabilities. We identify that the current parallel to the magnetic field ( j (cid:107) ), in combination with the break-downof general relativistic force-free electrodynamics across current sheets, impairs the physical modeling of resistive instabilities. Wefind that at least eight numerical cells per characteristic size of interest (e.g. the wavelength in plasma waves or the transverse widthof a current sheet) are needed to find consistency between resistivity of numerical and of physical origins. High-order discretizationof the force-free current allows us to provide almost ideal orders of convergence for (smooth) plasma wave dynamics. The physicalmodeling of resistive layers requires suitable current prescriptions or a sub-grid modeling for the evolution of j (cid:107) . Key words.
Magnetic fields - Methods: numerical - Plasmas
1. Introduction
The numerical modeling of the dissipation, transport and emis-sion of high-energy particles by strongly magnetized plasma isa necessary ingredient to the theoretical interpretation of highlyenergetic astronomical phenomena associated with compact ob-jects such as neutron stars (NS) or black holes (BH). Currentsheets seem to be one important location where such processestake place. In these, reconnection of the magnetic field results inthe acceleration of relativistic particles (see, e.g., Ball et al. 2019;Petropoulou et al. 2019; Kilian et al. 2020, for some recent ex-amples) as well as locations of Fermi-type processes (Guo et al.2019).When the magnetic di ff usivity (or resistivity, η ) is su ffi cientlylarge, plasma dynamics change due to non-ideal processes andthe plasma can no longer be modeled as ideal, i.e., a perfectlyconducting plasma. Astrophysical plasma is typically an envi-ronment of extremely low magnetic di ff usivity. Thus, resortingto ideal (relativistic) magnetohydrodynamics (MHD) or force-free electrodynamics (FFE) is a reliable assumption when mag-netic fields dominate all plasma dynamics. However, a numericaltreatment of the challenges at hand requires introducing a con-trolled (often implicit) amount of numerical di ff usivity. Such dif-fusivity of numerical origin stems from two main sources: First,from the discretization of a set of physical (balance) laws, oftenpartial di ff erential equations. Secondly, from the need of stabi- (cid:63) [email protected] lizing numerical solutions across the various types of disconti-nuities existing in ideal MHD or FFE.The design of numerical codes commonly focuses on min-imizing the amount of numerical di ff usivity across discontinu-ities. Howewer, much less weight is placed on minimizing or(at least) characterizing the numerical di ff usivity resulting fromthe discretization of the governing equations (with considerableexceptions, e.g., Rembiasz et al. 2017, in Eulerian MHD, orObergaulinger & Aloy 2020, in Newtonian hydrodynamics). Wefind that a thorough account of the numerical di ff usivity of al-gorithms designed to solve the equations of FFE or general rel-ativistic FFE (GRFFE) is of utmost importance. We thereforeassess whether numerical di ff usivity behaves as a physical di ff u-sivity and if it introduces pathological biases in the modeling of(astrophysical) plasma. This is a necessary step to grade the qual-ity of numerical results when the FFE approximation reaches itslimit. Such limits are expected specifically at the location of cur-rent sheets, (smooth) Alfvén waves, or Alfvén discontinuities.The linear phase of the dynamics in pinched force-free cur-rent sheets has been modeled in FFE (Komissarov et al. 2007)and compared to particle-in-cell simulations (Lyutikov et al.2018). Still, conventional FFE codes are insu ffi cient to accountfor the non-ideal electric fields that drive complex dynamics. Incontrast to the limits of non-vanishing particle inertia (MHD),energy is dissipated by violations of the algebraic constraints un-der which the FFE approximation is valid. In this second part ofour two-series publication describing our new GRFFE code, weaim to quantify the impact of such non-ideal dissipation in both Article number, page 1 of 14 a r X i v : . [ phy s i c s . c o m p - ph ] J u l & A proofs: manuscript no. B_code_paper smooth plasma waves and current sheets as thin current-carryinglayers across which the magnetic field changes either directionor magnitude (Harra & Mason 2004). With a view to a com-prehensive numerical modeling of (astrophysical) plasma acrossmagnetization regimes we ask the questions: Why do FFE meth-ods fail to resolve the dynamics in current-dominated domains?What are the physical consequences of these failures?This second manuscript in our series of publications is dedi-cated to answer these questions, to evaluate previous results ob-tained by our code, and to provide leverage points at which ki-netic plasma modeling may bridge the limits of FFE. In the con-text of instabilities in magnetar magnetospheres Parfrey et al.(2013); Carrasco et al. (2019); Mahlmann et al. (2019) have en-countered current sheets at the stellar surface, and sensitivityto a (conservative) transport of charges throughout the domain.The longest variability timescales in the TeV emission observed,e.g., in M87 have recently been linked to recurring periods ofe ffi cient Blandford / Znajek (Blandford & Znajek 1977) type out-flows induced by the accretion of magnetic flux tubes (Parfreyet al. 2015; Mahlmann et al. 2020b). Reconnection and plasmoidformation in BH accretion processes are likely to act on muchshorter timescales (studied in the ideal limit by, e.g., Nathanailet al. 2020) and involve relevant physical non-ideal electric fields(analyzed in the resistive limit by, e.g., Ripperda et al. 2020). Alarge array of work makes use of numerical laboratories set up in(GR)FFE in order to simulate the most extreme environments ofthe universe while constantly breaking their own limits. Under-standing the pathology of such breaches, and provide valid testcases for their evaluation has motivated this additional technicalexploration of our recent code development e ff ort.This work is organized as a series of papers. This manuscript(Paper II) characterizes the numerical resistivity of our GRFFEcode (as introduced in Mahlmann et al. 2020a, Paper I) in depthby studying the 1D di ff usion of plasma waves (Sect. 2.1) and thegrowth of 2D tearing modes (Sect. 2.2) under force-free condi-tions. In Sect. 3 we explore a phenomenological current, e ff ec-tively allowing for a modeling beyond ideal FFE, i.e., driving theevolution with a physical resistivity η . We discuss the implica-tions of our results on GRFFE methods in Sect. 4 and presentour conclusions in Sect. 5. We use units with G = c = M (cid:12) =
2. Assessment of Numerical Resistivity
In this section, we quantify the numerical resistivity of theGRFFE code presented in Paper I. This analysis is analogousto the one in Rembiasz et al. (2017) for Eulerian MHD codes,but applied to the case of FFE. For our analysis, we employ twokey techniques: i) Measuring damping rates of plasma waveswhich are captured in a one-dimensional periodic domain fora large number of dynamical timescales (Sect. 2.1). ii) Measur-ing growth rates of tearing mode (TM) perturbations in a two-dimensional force-free current sheet (Sect. 2.2). A similar strat-egy has been followed in Miranda-Aranguren et al. (2018) inresistive relativistic MHD. Since the development of TMs re-quires, at least, a longitudinal current sheet in 2D (though, obvi-ously, they can also develop in 3D), the study of the growth rateof resistive TMs allows us to quantify the numerical resistivityof our code in more than one dimension.
One of the simplest one-dimensional assessments of numericalresistivity is to check the capacity of the code of maintaining standing or propagating waves over significant numbers of light-crossing times (Rembiasz et al. 2017). Hereafter, quantities with-out any subscript refer to the total electromagnetic field, whichwe assume to be composed of a background field (annotated withsubscript 0) and a perturbation (annotated with subscript 1), i.e., B = B + B , D = D + D , ρ = ρ + ρ , (1)which correspond to the magnetic field, the electric field and thecharge density, respectively. These are the main evolved vari-ables in our code, apart from the two potentials that we use topreserve the value of the divergence of the fields (which we al-ways initialize to zero, see Paper I). We perform three series of1D simulations for di ff erent kinds of waves with the goal of mea-suring the damping rate induced by the existence of a finite nu-merical di ff usivity.Series A corresponds to a standing wave in FFE, which isexcited by initialising the variables to B = D = (0 , , , ρ = , B = ( B sin k y, , B cos k y ) , D = (0 , , , (2)where k is the wavenumber. ρ is initialized numerically to en-sure ∇ · D = ρ (here and in all the tests hereafter). B isthe scale-free amplitude of the standing wave, which we set to B = B along the y -direction, whichare excited by the following initial perturbation (Punsly 2003): B = (0 , B , , D = (0 , , , ρ = , B = ( B (cid:15) cos k y, , , D = (0 , , B (cid:15) cos k y ) . (3)Here, the amplitude of the perturbations is controlled by the di-mensionless parameter (cid:15) . We choose the perturbation amplitudesu ffi ciently small, (cid:15) = − , to guarantee that the perturbationcan be regarded as linear with respect the FFE equations. Thiswill allow us to measure the numerical resistivity from the damp-ing rate of the waves (see below).Series C corresponds to Alfvén waves traveling along thex-direction at an angle θ = π/ B = ( B / √ , , B / √ , D = (0 , , , ρ = , B = (0 , B (cid:15) cos kx , , D = B (cid:15) √ kx , , − cos kx ) . (4)where (cid:15) is of the same value as in series B. In all series, the 1Ddomain is discretized in the y -direction with an extent L =
16 inthe domain [ − L / , L / ff erent cases isevolved for 125 light crossing times ( τ = L ) of the domain. Weimpose periodic boundary conditions, so we chose a wavenum-ber k = π m / L , with m ∈ N . Simulations of the cases m = , ff erentnumerical resolutions to study the convergence of the method.We use the monotonicity preserving (MP, Suresh & Huynh 1997)reconstruction with three di ff erent orders: MP5, MP7 and MP9.In all cases we use a CFL factor of f cfl = . D , of Alfvén and fast magnetosonicwaves is proportional to the square of the wavenumber and to thedi ff usivity, ξ , i.e. D = k ξ. (5) Article number, page 2 of 14. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Characterization of Numerical Di ff usivity ■ ■ ■ ■■ ■ ■ ■■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ MP5 ( r = ) MP7 ( r = ) MP9 ( r = ) ■ m = ■ m = ■ m = - - N u m e r i c a l r e s i s t i v i t y ( m ⨯ η * ) Static diffusion ■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■
MP5 ( r = ) MP7 ( r = ) MP9 ( r = ) ■ m = ■ m = ■ m = Fast waves ■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■■ ■ ■ ■
MP5 ( r = ) MP7 ( r = ) MP9 ( r = ) ■ m = ■ m = ■ m = Alfven waves
Fig. 1.
Assessment of numerical resistivity (normalized, see Eq. 14) as a function of grid points per wavelength for a 1D static di ff usion test(series A) as well as propagating 1D plasma waves (series B / C, fast / Alfvén). Numerical fit parameters (Table 1) are obtained for the data pointsof the m = j (cid:107) hassignificant impact on the order of convergence. Thick lines in the background employ standard fourth-order finite di ff erences in the calculation ofreconstruction of j (cid:107) , while solid thin lines display the improved high-order finite di ff erences (Sect. 3.4, Paper I). - - - B r e f z / B ϵ - - B r e f z / B ϵ - - / τ B r e f z / B ϵ ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ m = = = / τ N o r m a l i z e d f r eee n e r g y U / ( B ϵ ) ² N y = ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ m = = = / τ N y = Fig. 2.
Exemplary amplitude evolution of selected modes of the fast wave model (showing m = m =
2, as well as the superposition m = , τ ). Left:
Field amplitude at the point y = m = , m = m = ff erent resolutions are visualized by solid ( N y =
12) and dotted ( N y =
24) lines.
Right:
Evolution of the energy contained in the perturbations,given by Eq. (8), scaled to (cid:15) B (squares) and theoretical model (solid lines) employing the fit parameters from Table 1. Mixed modes evolveaccording to a superposition of two exponentially damped plasma waves with the appropriately scaled damping rates inferred from Fig. 1. Here, the di ff usivity is a linear combination of the shear and bulkviscosity as well as the resistivity, η . Expression (5) strictly holdsin the so-called weak damping approximation, i.e., k ξ / (4 v ) (cid:28) v A is the Alfvén velocity).Komissarov et al. (2007) has shown that the equations ofFFE in the so-called incompressible limit can be cast into a formnearly identical to that of incompressible MHD. In the incom-pressible limit, the drift speed V = D × B / B fulfils ∇· V =
0. Formodels of the series A, V = ∇ · V = V = O ( (cid:15) ) ≈ ∇ · V = O ( (cid:15) ) ≈ (cid:15) . Hence, for our choice of (cid:15) , we can safely make the assumption that the damping of Alfvén and fast modes in in-compressible FFE proceeds analogously to their incompressibleMHD counterparts.Since FFE neglects the thermal and inertial contributions ofthe plasma, its equations are not a ff ected either by bulk or byshear viscosity (neither of physical nor of numerical origin).Thus, the di ff usivity may only come from resistive e ff ects. Pro-vided that in FFE there is no physical resistivity, the di ff usivityin our numerical results comes exclusively from the numericalresistivity, η ∗ . Following Campos (1999) with the notation ofRembiasz et al. (2017), we find ξ = η for all three series of Article number, page 3 of 14 & A proofs: manuscript no. B_code_paper initial data, D = k η ≡ k η ∗ ⇐⇒ η ∗ = k D . (6)Hence, the di ff usion of one-dimensional plasma waves proceedsaccording to B ( x , y, z , t ) = B ( x , y, z , t ) e − D t , (7)where B denotes the magnetic field of the ideal solution ofwave dynamics in the absence of any numerical dissipation. Thedi ff usion of the electric field D can be modeled analogously.In order to evaluate the di ff usive properties of our code, wedefine the electromagnetic energy contained in the wave compo-nents ( B and D ) U = (cid:90) d V (cid:104) B + D (cid:105) ≡ U e − D t . (8)Since the energy U is associated with the perturbations of thebackground electromagnetic field, we may refer to it also asfree energy. In our analysis, we evaluate the right-hand side ofEq. (8), for which we obtain the linear relationln U = − D t + ln U , (9)where U is the initial energy in the respective waves. Eq. (9)allows us to obtain D (and hence, η ∗ ) from the slope of linear fitsof the form ln U vs. t . The values of η ∗ computed from the fitsdepend on the grid spacing with which the model is run or, equiv-alently, on the number of points per wavelength, p = N / m ( N isthe number of points in the domain and m is the number of fullmodes), with which a target mode is resolved in a given setup.Fixing the remaining parts of the numerical method, the value of η ∗ also depends on the cell interface reconstruction employed.Following Rembiasz et al. (2017, Sect. 4.3.3), we assume thatthe numerical resistivity can be written in the form η ∗ = N × V × L × (cid:32) ∆ x L (cid:33) r , (10)where N is a (resolution-independent) numerical coe ffi cient, r the measured order of convergence of the employed scheme, ∆ x = L / N the grid spacing, and L and V are the characteris-tic length and speed of the problem. In the case of FFE, di ff er-ently from incompressible MHD, the characteristic velocity isthe speed of light, thus V = L = λ = π/ k = L / m = p ∆ x . Thus, we canrewrite expression (10) as η ∗ = N × λ × p − r . (11)In view of the previous relation, we define subsets of tests wherethe cell interface reconstruction is fixed and di ff erent number ofpoints per wavelength are considered. For every subset of tests(i.e. reconstruction method), we fit the functionln η ∗ = − r ln p + d , (12)where the fit parameter d corresponds to d = ln [ N × λ ] = ln [ N × L ] − ln m . (13) Table 1.
Assessment of numerical resistivity (numerical fit parametersof the data visualized in Fig. 1, see Eq. 14) for a 1D standing sinewave (series A) as well as propagating 1D plasma waves (series B / C,fast / Alfvén). In parentheses we show the results for the cases that dif-fer significantly when using the standard fourth-order reconstruction forthe parallel current.
Series Reconstruction N rA MP5 34 .
96 4 . .
13 6 . .
03 8 .
79B MP5 34 .
96 4 . .
34 6 . .
37 8 .
38C MP5 36 .
47 4 . .
91 6 . .
52 8 . (cid:2) m × η ∗ (cid:3) = − r ln p + ln [ N V L ] . (14)We therefore combine measurements of di ff erent wavemodes m in the same numerical domain as a measure of nu-merical resistivity. For the entire set of considered cases (seriesA / B / C) we find that at least eight grid points per wave mode arerequired in order to resolve the respective plasma wave with anorder of convergence which approaches the theoretical order ofthe employed reconstruction scheme. Table 1 shows the coe ffi -cients N and r obtained by a linear fit to Eq. (14). For MP5 theirvalues are similar (as expected from Eq. 11) independent of theseries. The numerical di ff usivity in series B is somewhat largercompared to series A and C (the latter with the eighth-order ac-curate discretization of j || ) for the set of models using the MP9reconstruction. We have not been able to identify a single, dom-inant source for the slight reduction of the order (and slightlyincreased di ff usivity) of MP9 when applied to fast waves. How-ever, the (cid:46)
5% reduction of r does not seem to be a great draw-back for applying MP9 in combination with either 4th or eighth-order accurate discretization of j || to numerical models involvingfast waves.In the right panel of Fig. 1 we also show the e ff ects of theorder of discretization of the derivatives appearing in the parallelcurrent (Sect. 3.4, Paper I). The thick (opaque) lines show theresults when a fourth-order accurate discretization is employedin combination with MP7 (blue thick line) or MP9 (black thickline) reconstruction. For comparison, the blue (black) thin linesshow the fits to the resolution dependence of the numerical re-sistivity using a sixth(eighth)-order accurate parallel current dis-cretization in combination with MP7 (MP9). In order to obtainan empirical order of convergence close to the theoreticaly ex-pected values ( r = j || is crucial. We notethat for resolutions of less than 8 points per wavelength, the val- Article number, page 4 of 14. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Characterization of Numerical Di ff usivity ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■■ MC ■ MP5 ■ MP7 ■ MP910 20 50 1005. × - ( p ) per current sheet width ( a = ) G r o w t h r a t e γ T M d = / m = - = / m = - = / m = - = - / m = - - - - - - - N u m e r i c a l r e s i s t i v i t y η * Fig. 3.
Growth rate of the (2D) tearing mode instability (Sect. 2.2) fordi ff erent resolutions and reconstruction schemes. We depict the mea-surements of γ tm from numerical experiments (squares) and the corre-sponding best fit parameters (solid lines, Eq. 24). The scale on the rightshows the numerical resistivity computed from the measurements of thegrowth rate employing model A (Eq. 19). ues of the numerical resistivity are not too di ff erent. This jus-tifies our choice of a fourth-order accurate discretization of j || (Mahlmann et al. 2019, 2020b) in combination with MP7. Theresolution employed in these global 3D models is smaller than p ∼ p = p =
20) grid zones in order to reducethe numerical resistivity below ∼ − ( ∼ − ). Extrapolat-ing these requirements to 3D numerical experiments reveals thehuge computational e ff ort needed to resolve the dynamics of theinteraction of FFE modes with su ffi ciently small numerical re-sistivity.Figure 2 gives an illustrative example of mode damping fordi ff erent resolutions per wavelength, including the superpositionof two selected modes of fast waves with the same initial am-plitude (adding the initial data in Eq. 3 for m = m = m = m = TM instabilities are resistive instabilities, which can develop incurrent sheets, and dissipate magnetic energy into kinetic energyif a plasma fluid is considered. Note, however, that in a pureforce-free approximation, the only existing energy is that of theelectromagnetic field and, thus, the dissipation of magnetic en-ergy may simply result in a sink of the latter. In GRFFE, besides t t - - - × - × - v z ( x = . ) MP9 - - - × - × - v z ( x = . ) MP7 - - - × - × - v z ( x = . ) MP5 - - - - - z v z ( x = . ) MC Fig. 4.
Evolution of the drift velocity V z along the mid point of thecurrent sheet ( x = .
5) in the direction transverse to the current sheet( z -direction) during the linear phase (time-interval [ t , t ] for which thegrowth rate is derived in Fig. 3). The di ff erent panels correspond to dif-ferent reconstruction schemes (see legends) and a resolution of p = . assumed width of the resistive layer L R employed for the analysis in Sect. 2.2 isindicated by the vertical grey lines in the background; each of these linesdenotes the limits of the computational zones in the z − direction . Thegrowth of the velocity component V z is strongly suppressed for higher-order (MP) reconstruction methods and is significantly larger when us-ing the MC reconstruction (note the di ff erent scale for the lower panel). the standard sources of numerical di ff usivity, there is yet anotherone, namely the numerical resistivity induced by the algorithmsused to control the violations of the force-free conditions. Aswe shall see, both the standard sources of numerical resistivityas well as the resistivity produced by the violation of the force-free constraints are especially sensitive to the resolution of thenumerical mesh.In this test, we adapt the relativistic resistive MHD currentsheet setup of Del Zanna et al. (2016) to FFE, employing a two-dimensional domain of ( x , z ) ∈ [ − L / , L / × [ − a , a ]. Here, a = . L = π/ k is its Article number, page 5 of 14 & A proofs: manuscript no. B_code_paper
Table 2.
Estimates of the dimensionless coe ffi cient N and the empir-ical order of convergence r for the two models defined by Eqs. (19)and (21). The parameters are obtained by the linear fit given by Eq. (24).In the last column, we display the relative deviation of the fit parameter m (see Fig. 3) from its analytical value m a (for the latter employing aresolution-dependent length scale L , Eq. 33, and Rembiasz et al. 2017)in the corresponding limit and assuming ideal order of accuracy of thescheme (see Sect. 4). Limit Reconstruction N r ∆ m / m Model A MC 4 . × − .
95 0 . . × − .
31 0 . . × − .
57 0 . . × − .
84 0 . . × − .
34 0 . . × − .
78 0 . . × − .
08 0 . . × − .
41 0 . k = π . Thisvalue allows us to fulfill approximately L (cid:29) a and ka (cid:28)
1, re-quired for the development of TMs. We use periodic boundaryconditions in the x direction (i.e., along the current sheet) andcopy boundary conditions in the z direction (i.e., in the directionperpendicular to the current sheet). The numerical grid is uni-form and consists of N x × N z zones, where N x =
32 in all cases.In order to trigger TM instabilities, the initial magnetic field B x = B tanh ( z / a ) , B y = B sech ( z / a ) , (15)is perturbed by B x = (cid:15) ( ak ) − B sin ( kx ) tanh ( z / a ) sech ( z / a ) , B z = (cid:15) B cos ( kx ) sech ( z / a ) . (16)We set B = (cid:15) = − . We assume that the initial electric field and the charge den-sity are zero. The growth rate of the TMs, γ tm , may be traced,e.g., by examining the growth of the magnetic field compo-nent B z (cf. Rembiasz et al. 2017), which grows exponentially, B z = B z ( t = e γ tm t . To obtain a globally and positively definedquantity for the growth rate, we study the integral of B z over asuitably chosen patch of the computational domain (covering theentire length and width of the current sheet):ln (cid:90) B z d S = γ tm t + ln (cid:90) B z d S . (17)The slope 2 γ tm of the linear relation in Eq. (17) may be obtainedby fitting ln (cid:82) B z d S vs. t . In the fits we disregard the initial (nu-merically dominated) adjustment phase. We measure the growthrate for series of simulations with di ff erent numerical resolu-tions and numerical schemes. Apart from the MP reconstructionschemes used in the previous sections we also test slope lim-ited TVD reconstruction with a monotonized central (MC) lim-iter (van Leer 1977).We aim to provide an estimate of the numerical resistivityas a function of the numerical resolution based on two di ff erent approximations. The first one (model A) assumes that the plasmais inviscid. In this case, the growth rate of the TM mode (for agiven mode k ) depends on the physical resistivity as (Rembiaszet al. 2017, Eq. 147) γ tm = . − / × η / v / A a − / ( ak ) / (cid:32) ak − ak (cid:33) / . (18)In our system of units, the Alfvén speed is v A =
1. For thespecific choice of k and a employed in our setup, we can writeEq. (18) as γ tm ( k = π, a = . ≈ . × η / . (19)Alternatively, in the so-called long-wavelength approximation(model B, characterized by ka (cid:28) γ tm , max ≈ . × a − / v / A η / , (20)or, equivalently, γ tm , max ( a = . ≈ . × η / . (21)The force-free models we run do not include any physical resis-tivity. Hence, the growth of TM modes is induced by the actionof the resolution-dependent numerical resistivity η ∗ . Thus, wemay replace η in Eqs. (19) and (21) by η ∗ .Once the TM growth rate is obtained for each resolution andreconstruction method, we can compute the corresponding nu-merical resistivity (as we have done in Eq. 10, cf. Rembiaszet al. 2017, Sect. 4.3.3). Like in the previous section, the lightspeed ( V =
1) is the only possible choice for the characteristicvelocity of the problem in the force-free regime. The selectionof L is much more involved (see Rembiasz et al. 2017, for adetailed discussion) and prone to accuracy restrictions since typ-ically L (cid:28) a , which can result in a prohibitive numerical reso-lution needed to reliably measure the growth rate. To obtain anorder of magnitude estimate, we may assume that L ∼ . a , (22)but cautiously note that assuming that L is constant, i.e.resolution-independent, is not fully accurate (Rembiasz et al.2017). A partial justification for this choice (but see discussionin Sect. 4) is the shape of the drift velocity, V = D × B / B , dis-played in Fig. 4, which shows that 0 . a (region with grey lines)covers a significant part of the width of the current sheet (esti-mated as the distance between the two extrema).We express the numerical resistivity in terms of the numberof grid zones within the transition layer, i.e., p = a / ∆ z ( ∆ z be-ing the grid spacing in the z direction). Note that p here has adi ff erent meaning as in Sect. 2.1, and the numerical resistivity(Eq. 10) in these tests can be written as (compare with Eq. 11) η ∗ = N × . a × (0 . p ) − r , (23)where N is, again, a (resolution-independent) numerical coef-ficient. Plugging Eq. (23) into either Eq. (19) or Eq. (21) oneobtains a linear relation between ln γ tm and ln p that allows tocompute N and r from the coe ffi cients of the linear fitln γ tm = m ln p + d , (24)where the specific definition of m and d depends on the employedgrowth model. Article number, page 6 of 14. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Characterization of Numerical Di ff usivity Figure 3 visualizes the growth rates γ tm for di ff erent resolu-tions and reconstruction methods. We summarize the numericalfit parameters, i.e. the coe ffi cient N and the estimated order ofthe scheme r in Table 2. The order of convergence estimatedfrom the fits is roughly correct for models using the MC recon-struction, but is significantly smaller than the formal order ofconvergence estimated for MP methods, where values of r = r .The reason for the discrepancy is related to the nature of theTM itself, which involves dissipation in a current sheet. Dissi-pative e ff ects do not exist in FFE (except numerically) and thepresence of these current sheets breaks the force-free conditionitself (see the discussion in Sect. 4). In this regime, it is not sur-prising that our numerical methods do not behave as expectedand, hence, the theoretical order of convergence is not recov-ered. These circumstances are not present in the tests in Sect. 2.1where no current sheets are on hand and no reduction of the ex-pected order of convergence is observed. A possible solution tothis problem is to include a physically consistent model for thedissipation in FFE, such that TMs can develop in a physicallyrealistic way forming resistive layers resolvable by numericalsimulations. Its application to di ff erent astrophysical scenarios(where the resistivity is likely too small to be handled) wouldbe done by studying the limit of decreasing resistivity. The nextsection explores this possibility.The influence of the exact boundary conditions on the growthrate in the direction transverse to the current layer is probablysignificant and deserves further study, something that is beyondthe scope of this publication. We have considered di ff erent val-ues of a , finding that a = . z -direction lessthan ± a induce significant modifications in the growth rate ofTMs. Our choice of a box size extending up to z = ± a has beenvalidated by cross-checking against even larger computationalboxes and finding no appreciable changes in the growth rate. Wehave also considered higher-order discretizations of the deriva-tives used to compute j || (Sect. 3.4, Paper I), finding no signifi-cant influence in this test, likely because of the non-smoothnessof the parallel current at z =
0. The discontinuity of the fieldderivatives makes a high-order discretization of j || less e ff ective(i.e., more noisy). This result is consistent with our finding thathigher-order discretization of j || does not significantly change theresults of the 1D Riemann problems presented in Sect. 4.1 of Pa-per I.
3. Beyond Ideal Force-Free Electrodynamics
FFE is inherently formulated in the limit of infinite conductivity,e ff ectively resulting in the magnetic dominance and ideal fieldscondition (see Paper I and also Lyutikov 2003) B − D ≥ , (25) D · B = . (26) Under these conditions the 4-vector electric current density takesthe form of the force-free current (Eq. 47 in Paper I): I µ = I µ ff = ρ n µ + ρ B η νµαβ n ν D α B β + B µ B η σαβλ n σ (cid:16) B λ ; β B α − D λ ; β D α (cid:17) , (27)where ρ = − n µ I µ = α I t is the charge density, n µ = α − (1 , β i ) isthe vector normal to the hypersurfaces in the 3 + α is the lapse function and β i is the shift vector,respectively. Equation (27) is a direct consequence of L n ( B · D ) = n λ ∇ λ ( B · D ) = , (28)where L n is the Lie derivative with respect to n µ . The cur-rent density 3-vector appearing in the 3 + J i = α I i . Addition-ally one can define the current density 3-vector for the normalobserver ( n µ ) as the projection of I µ onto the hypersurface, i.e. j i = I i − β i I t , which is related to the former by J i = α j i − β i ρ .Any occurring resistivity measured in such ideal schemesmust, hence, be of numerical origin; we quantify this numer-ical resistivity for our method in Sect. 2. The absence of anyphysical resistivity makes resolving genuinely resistive layers atrue limit of FFE. Thus, it is desirable to extend the theory to in-clude a small phenomenological resistivity η and reduce to FFEin the limit η →
0. Several attempts to undertake this task havebeen presented in the literature (e.g., Lyutikov 2003; Gruzinov2007; Li et al. 2012; Parfrey et al. 2017). In spite of the oxy-moron, it is common to refer to these prescriptions as resistive
FFE. All of them have in common that they integrate the fullset of Maxwell’s equations together with a suitable choice ofOhm’s law. In practice, this replaces the force-free electric cur-rent density I µ ff by a more general form explicitly including theresistivity. For this section, we have tested the generalization ofa current density recently introduced by Parfrey et al. (2017). Itscovariant form (cf. Eq. 27) is I µ = ρ n µ + ρ B η νµαβ n ν D α B β + + κ I η B µ B (cid:104) η σαβλ n σ (cid:16) B λ ; β B α − D λ ; β D α (cid:17) + κ I B σ D σ (cid:105) , (29)where κ I is a constant and η is the resistivity. With this new 4-current the evolution of B · D fulfills (cf. Parfrey et al. 2017) L n ( B · D ) = κ I (cid:16) α − η J − D (cid:17) · B = κ I (cid:0) η j || − D || (cid:1) · B − κ I ηα − ρ βββ · B , (30)where || indicates the component of a vector parallel to B . There-fore, 1 /κ I is the timescale over which the parallel electric field D || is driven towards α − η J || (or in for the case of flat spacetime,towards η j || ). Therefore, the form of this current e ff ectively isenforcing a form of Ohm’s law for the parallel current. In thelimit η → B · D → I µ → I µ ff .The current density observed by the normal observer nat-urally splits into components parallel and perpendicular to themagnetic field 3-vector ( j = j ⊥ + j || , cf. Eq. 70 in Komissarov2011): j ⊥ = ρ D × BB , (31) j || = + κ I η B · ( ∇ × B ) − D · ( ∇ × D ) + κ I B · DB B . (32) Article number, page 7 of 14 & A proofs: manuscript no. B_code_paper
In the following subsections, we will investigate the abilityof the current given by Eq. (29) to model resistivity in FFE. Wespecifically address its ability to a) resolve current carrying dis-continuities in 1D plasma wave tests (extending Sect. 4.1 of Pa-per I) as well as b) its capacity to resolve the resistive layer in2D tearing modes (cf. previous Sect. 2.2). Along the way, wecross-validate our procedure to evaluate the numerical resistivityof our algorithm by comparing this prescription for resistive FFEwith its ideal limit.From the numerical point of view, the introduction of re-sistive e ff ects needs to be done with extra care. The choice ofthe driving rate κ I is guided by numerical convenience. Its valueshall be large enough to drive D || towards α − η J || as quickly aspossible., however, it shall not be so large as making the driv-ing time scale become much smaller than the time step of thenumerical model (1 /κ I (cid:28) ∆ t ). In this case, Maxwell’s equa-tions become numerically sti ff and suitable time evolution algo-rithms such as implicit Runge-Kutta methods (Palenzuela et al.2009; Miranda-Aranguren et al. 2014, 2018) would need tobe employed. Numerical experimentation leads us to resort to κ I ∆ t ∼ −
10. Note that if the modified current is employed, theperpendicularity condition D · B = B · D towards zeroin FFE simulations (with η = The three-wave test in Sect. 4.2.1 of Paper I is a well-suitedtestbed for the analysis of resistivity introduced by Ohm’s law.This test consists of an initial discontinuity that splits into twofast discontinuities and one stationary central standing Alfvénwave. The Alfvén wave is a charge carrying discontinuity, whichis e ff ectively stabilized by strong currents. In this section, wecompare the results of the three-wave test using three di ff erentOhm’s laws. In all tests we use the highest resolution employedin Paper I, namely, ∆ x = . η = j (cid:107) = κ I = – The force-free current I µ ff preserves the central charge car-rying discontinuity well. Setting j (cid:107) = x = – The resistive Ohm’s law implementing I µ approaches theforce-free solution for decreasing values of η . For η = − ,the largest considered value, the charge-carrying discontinu-ity shifts to the right of the central point. Decreasing theresistivity to smaller values η = − almost reproducesthe ideal FFE results: Both charge and stabilizing currentsaround the central point gradually increase.The fact that decreasing values of the phenomenological resis-tivity η approach the results from ideal FFE is showing (reas-suringly) that the numerical resistivity in our ideal formulationof FFE is below 10 − even in charge-carrying discontinuities. The results presented in Fig. 5 also include somewhat counter-intuitive observations: Setting j (cid:107) = η →
0, Yu 2011) showsa similar behavior as the case of η = − (hence, η (cid:29)
0) inthe resistive modeling. We attribute this contradictory behaviorto the feature of enhanced charge conservation, which we ensurein our code (see Paper I). Abandoning part of the current ( j (cid:107) ) andonly enforcing D · B = ρ and the electric field D . As we evolve thecontinuity equation of the charge density (and do not reconstruct ρ = div D in every timestep, as most other FFE codes do) thealgebraic reset of the electric fields to enforce D · B = In this section, we repeat the TM tests from Sect. 2.2 for the re-sistive current I µ from Eq. (29). In all cases we use κ I = η , set in the expressionfor the current, and the numerical resistivity η ∗ , depending onthe numerical scheme and resolution (as seen in Sect. 2.2). Forsu ffi ciently high resolution, one should obtain convergent resultsdepending only on η . If the resolution is decreased, the point atwhich the growth rate γ tm starts di ff ering from the converged re-sult should correspond to η ∗ ∼ η , and marks the limit at whiche ff ects of numerical resistivity become dominant. This allowsfor an independent estimation of the numerical resistivity of thecode, to be compared with those in Sect. 2.2.Figure 6 assembles growth rates for di ff erent values of theresistivity parameter η . If the imposed physical resistivity iswell above the numerical resistivity quantified in Sect. 2.2, i.e., η ∗ < η , the measured growth rates γ tm converge to a single valuealmost independently of the employed resolution and reconstruc-tion method. As η is reduced approaching the limit set by thenumerical resistivity, the growth rates become more scatteredaround the corresponding converged value. In cases where nu-merical and imposed resistivity are comparable, η ∗ ≈ η , (e.g.for MC reconstruction with the lowest η ) the growth rate can beexplained by the combined e ff ect of numerical and physical re-sistivity. Finally, we expect that in case of η ∗ > η , the evolutionis fully dominated by numerical resistivity.The solution in the regime with finite resistivity di ff ers inqualitative ways with respect to the ideal FFE case. Accordingto (Low 1973) the velocity component V z is driven to a finitevalue η/ a asymptotically by the employed isotropic resistivitymodel. We observe this feature in Fig. 7, where the drift veloc-ity approaches asymptotically the dashed gray line representing η/ a . This behavior di ff ers from the ideal case studied in Sect. 2.2,where V = ffi ciently high numerical res-olution, as seen in Fig. 4). The key di ff erence between both casesis that the numerical resistivity measured in Sect. 2.2 acts mainlyin the region of the current sheet, while in the resistive simula-tions in this section η is the same over the whole domain. There-fore, the (physical) boundary conditions of both sets of simula- Article number, page 8 of 14. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Characterization of Numerical Di ff usivity j ∥ : Force - Free j ∥ = B y MP7 - - - - - - ρ j ∥ j ⊥ η = - η = - η = - η = - η = - - - B y MP7 - - - - - - - -
20 x ρ - - j ∥ - - j ⊥ Fig. 5.
Three-wave test ( ∆ x = . ff erent Ohm’s laws. We display the magnetic field component B y , the electric charge ρ , aswell as the magnitudes of the parallel ( j (cid:107) ) and perpendicular ( j ⊥ ) projections of the current. Top panels:
Conventional force-free current I FF (blackcolor, Eq. 27) and j (cid:107) = Bottom panels:
FFE current (Eq. 29) for di ff erentvalues of η . tions di ff er slightly (although in the limit η → ff er from the FFEresults. In Sect. 2.2 (see Fig. 4) we found that the size of theresistive layer was e ff ectively of the size of a few numerical cells.However, in the resistive simulations (see Fig. 7) the widths ofthe resistive layers (estimated as the extent between the extremaof the drift velocity profile) extends well beyond a few grid zoneso ff the central current sheet, especially for η > η ∗ . In the lattercase, the layer is well resolved. The di ff erence between the twoclasses of simulations is that in resistive FFE the width of theresistive layer is set by the value of η while in the ideal FFEcase its width is set by the numerical resistivity and shrinks withincreasing resolution. Therefore, the width of the shear layer inFFE is a moving target as we increase resolution and, dependingon the numerical scheme used, it may not be possible to resolveit with increasing resolution. For cases with η ≈ η ∗ , additional(local) maxima emerge corresponding to a thinner resistive layerassociated to numerical resistivity (see bottom panels of Fig. 7).The fact that the region of convergent results of η > η ∗ roughly coincides with the limits established in our previoustests further supports our findings from Sect. 2.2: For intermedi-ate resolutions of p ≈
10, our GRFFE tool has a numerical resis-tivity below η ∗ ≈ − . Furthermore, phenomenological Ohm’slaws as the one introduced in Eq. (29) are apt to model resistivee ff ects in FFE. Although we cannot expect, for the time being,to be able to work with the values of resistivity expected in mostastrophysical scenarios, at least we can aim at performing sim-ulations with controlled values of the resistivity and study itsbehavior in the limit of vanishing η . Notwithstanding these en-couraging results, using resistive FFE brings its own problems,such as the appearance of a non zero asymptotic drift velocitymoving away from the resistive layer in TMs. This indicatesthe need for further fine-tuning when employing resistive Ohm’slaws in FFE. Non-isotropic resistivities (as employed, e.g., in Komissarov 2004) or current-sheet-capturing models (as used inParfrey et al. 2017) are candidates for driving the asymptotic dy-namics of current sheet instabilities further towards the physicalsolution. Such fine tuning is very likely to depend on the em-ployed method and the problem at hand, i.e., the correct deter-mination of the threshold η ≈ η ∗ . A full implementation of dif-ferent resistive FFE models with non-constant or non-isotropicresistivities will be further explored in the future. Correctly mod-eling phenomenological Ohm’s laws could be a valuable assetwhen bridging codes of di ff erent plasma regimes as they allowfor physical transition layers (and consistent signal propagation).
4. Discussion
The discretization of any system of partial di ff erential equations(PDEs) on a finite mesh of cells to make them amenable for theirnumerical integration unavoidably introduces numerical di ff u-sion and dispersion. The equations of GRFFE are no exception tothis rule. Furthermore, the necessity of enforcing the constraintequations specific to GRFFE ( D · B = B − D >
0) isan additional source of numerical di ff usion. The use of mono-tonic and consistent numerical methods guarantees the conver-gence to the physical non-dissipative solution of a problem inthe limit of vanishing cell size. In this limit (which is unreachablein practice), the e ff ects of numerical di ff usivity are unimportant.However, for any finite cell size of practical use (especially in3D models), numerical di ff usion does not necessarily mimic thee ff ects of physical dissipation, particularly when the numericaldiscretization of the PDEs is performed employing high-ordermethods. It is therefore sound to assess the e ff ects of numeri-cal di ff usion as a source of dissipation specifically in (GR)FFE,which should, in theory, not incorporate any physical dissipation.We have performed a number of tests to quantify the amountof numerical dissipation as a function of resolution in our nu-merical algorithm. For that purpose, we have employed tests in Article number, page 9 of 14 & A proofs: manuscript no. B_code_paper ■ ■ ■ ■ ■■ ■ ■ ■ ■■ ■ ■ ■ ■■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆◆ ◆ ◆ ◆ ◆◆ ◆ ◆ ◆ ◆◆ ◆ ◆ ◆ ◆▲ ▲ ▲ ▲ ▲▲ ▲ ▲ ▲ ▲▲ ▲ ▲ ▲ ▲▲ ▲ ▲ ▲ ▲
10 15 20 250.10.20.5 Points ( p ) per current sheet width ( a = ) G r o w t h r a t e γ T M η = - η = - η = - η = - MP5MP7 MC
Fig. 6.
Measured growth rates for simulations using di ff erent resistivityparameters η in the resistive current density, Eq. (29), and for di ff er-ent reconstruction schemes (MC: magenta, MP5: red, MP7: blue, MP9:black) as a function of the numerical resolution. Thick lines in the back-ground show the fits obtained in Sect. 2.2 for numerical simulationswithout any added resistivity (see Fig. 3). As η is lowered, the e ff ectof numerical resistivity emerges (especially in the bottom left corner ofthe panel). Gray dashed lines mark the values of the growth rates cor-responding to the values of η used in the di ff erent series of simulations,computed using Eq. (19).
1D and 2D, but not 3D tests. This is, firstly, due to the lack ofstraightforward genuinely 3D analytic solutions including resis-tivity in the literature. Secondly, 3D tests are computationallymuch more expensive than 1D or 2D numerical experiments.Nevertheless, we assume that the multidimensional nature of 2Dtests presented should also shed light on the numerical dissipa-tion in 3D. Numerical experience tells us that the dissipation in-troduced by algorithms based upon a directional (split) integra-tion of the equations is only (very) weakly dependent on dimen-sionality. Backing up these assumption, we have shown here inthis work that the level of numerical resistivity depends on thenumber of zones per characteristic size to be resolved in everydimension (compare the values of η ∗ in Sect. 2.1 for 1D, Sect. 2.2for 2D), and not on the dimensionality of the problem.In the following sections, we address the similarities and dif-ferences between numerical and physical resistivity in FFE. Theaim of this discussion is to provide support for the interpretationof the numerical dissipation found in astrophysical applicationsof our code (Mahlmann et al. 2019, 2020b) as a consistent modelof physical dissipation. Along the way, we may also contribute toasses whether suitably extended FFE models (including physicaldissipation, albeit not at the astrophysically expected levels) maybe used to study a number of dissipative processes in relativisticastrophysics, e.g. resistive solutions for pulsar magnetospheres(Li et al. 2012), the dissipation induced by the kink instabilityin relativistic jets (Bromberg et al. 2019), or the conversion be-tween fast and Alfvén modes (Li et al. 2019). t t - - v z ( x = . ) η = - MP7 - - v z ( x = . ) η = - - - v z ( x = . ) η = - - - - - v z ( x = . ) η = - Fig. 7.
As Fig. 4 but showing one specific reconstruction method (MP7)for di ff erent resistivity parameters η . The asymptotic velocity is drivento a magnitude η/ a (gray dashed lines), di ff erent from zero, by theisotropic phenomenological resistivity employed in I µ . In physically driven TMs, the resistive e ff ects are restricted to theso-called resistive layer, with width L R (cid:28) a for inviscid, incom-pressible MHD plasma (e.g. Furth et al. 1963). Resolving theresistive layer numerically is challenging, due to the hierarchyof scales L R (cid:28) a (cid:28) L that must be spanned in these experi-ments (Rembiasz et al. 2017). L R depends on resistivity and TMgrowth rate (Rembiasz et al. 2017): L R ≈ . γ tm η ∗ a k v / = . a ( γ tm η ∗ ) / ( ka ) / , (33)where we have replaced the physical resistivity by the numericalone in account of the fact that the TMs in our models employ-ing the force-free current I µ ff are induced by numerical resistivity.Also, we have used that v A = L R on the product ( γ tm η ∗ ) / . Article number, page 10 of 14. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Characterization of Numerical Di ff usivity MCMP5MP7MP90 20 40 60 80 - - - - -
10 Time t l n ∫ B z ⅆ S ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ●
10 20 50 1000.20.512 Points ( p ) per sheet width L R / Δ z ■ Model A ● Model B
Fig. 8.
Left:
Growth of the (2D) tearing mode instability for di ff erenttreatments of the parallel current and a fixed number of points per cur-rent sheet width p =
6. We display our standard method for treating J || (i.e., with a fourth-order accurate discretization of the derivatives, solidlines), extreme cases in which J || = (dotted lines, following Yu 2011),and the case in which the algebraic constraint D · B = ff erent colors correspond to distinct reconstructions(see legends). In the case of no enforcement of the D · B =
0, the resultsare indistinguishable for the three MP reconstructions.
Right:
Size ofthe resistive layer L R normalized to the grid spacing ∆ z in the directionperpendicular to the current sheet for all the reconstruction methods (seelegends) and models (A, circles and B, squares). γ tm and the physical resistivity η depend on each other(Eqs. 19 and 21). However, the relation between γ tm and the nu-merical resistivity η ∗ is, a priori, unknown. This is the reason totest more than a single possible limit in Sect. 2.2 (i.e. to considermodels A and B). In both, γ tm → η ∗ → γ tm η ∗ in Eq. (33) tends towards a finite (non-zero) valueas the grid spacing is decreased (independently of the order ofspatial convergence of the method).The width of the resistive layer in our numerical experimentscan be approximately obtained by the half-distance between thetwo extrema of the component of the drift velocity associatedwith the growth of the TMs ( V z in our setup, see Fig. 4; cf. Rem-biasz et al. 2017, Fig. 10). We locate the position of these veloc-ity extrema (measured at x = .
5) within | z | (cid:46) ∆ z , and typically,at z = ± ∆ z . In practice, these results imply that for every resolu-tion we use, the numerical value L R (cid:46) ∆ z and, hence, there is nopossibility of finding a convergent behavior. The finer the gridspacing, the smaller L R . A finer grid spacing yields smaller η ∗ ,which translates into a shorter timescale of energy-momentumloss (see below). This prevents the resistive layer to fully reachthe conditions under which the growth of TMs would developtheoretically.In order to cross-validate these findings, we assume thatthe numerical TMs develop in regimes where the approxima-tions that hold for model A or model B are correct. PluggingEq. (19) or Eq. (21) into Eq. (33), one obtains two expres-sions for L R ( γ tm ), which implicitly depend on the grid spac-ing. For the two models A and B, the right panel of Fig. 8 shows the ratio L R ( γ tm ) / ∆ z for di ff erent resolutions, and di ff er-ent methods of cell interface reconstruction. The obtained values L R ( γ tm ) / ∆ z < direct measurement of L R (seeabove). Interestingly, the estimated value of L R ( γ tm ) / ∆ z dependson the spatial order of accuracy of the method. The highest-order reconstruction (MP9) yields the smallest γ tm η ∗ , while thesecond-order accurate MC method allows to resolve the resistivelayer width very marginally, as L R ( γ tm ) / ∆ z (cid:38) . p =
20 grid zones span the current width a .In the light of the previous considerations, we can now dis-cuss our assumption of a fixed characteristic length, indepen-dent of resolution (22). For that, we have also analyzed the re-sults following more closely the procedure of Rembiasz et al.(2017): We use L = L R , substitute expression (33) into Eq. (10),and isolate η ∗ as a function of γ tm , ∆ z , and the remaining pa-rameters of the problem. This expression η ∗ ( γ tm , ∆ z ) is thenprovided to Eq. (18), which yields a power-law dependence γ tm ∝ ∆ z r / (3 + r ) for model A. Analogously, for model B one ob-tains γ tm ∝ ∆ z r / (5 + r ) . Thus, the slope of the relation in Eq. (24)is related with the order of convergence r through m = − r / (3 + r ) ⇔ r = − m / (3 + m ) (model A) , (34) m = − r / (5 + r ) ⇔ r = − m / (4 − m ) (model B) . (35)Evaluating r using Eqs. (34) or (35) yields inaccurate results,since the assumption that L R is given by Eq. (33) is only approx-imated (as shown above). The degree of inaccuracy is roughlyquantified in the last column of Table 2, where we show the de-viation of the fit parameter m with respect to its analytical value, m a , assuming that r is equal to the formal order of convergenceof the respective method (i.e., we calculate ∆ m / m : = − m a / m ).The empirical order of convergence is very sensitive to the exactvalue of m . This sensitivity is rooted on the dependence of L R on ∆ z . We evaluate this dependence in our results in the rightpanel of Fig. 8. Repeating the exercise of the previous para-graph but isolating L R , one finds a dependence with resolution L R / ∆ z ∝ ∆ z m L , where m L = − / (3 + r ) (model A) , (36) m L = − / (5 + r ) (model B) . (37)Therefore, m L (cid:46) r . A direct comparisonwith the data in the right panel of Fig. 8 suggests a small positivevalue of the exponent m L , i.e., m L (cid:38)
0. We also highlight the factthat a value m L ≈ r from the slope of the linear fit in Eq. (24). Therefore, assumingthat the physical relation of Eq. (33) holds is not appropriate inour models.Admittedly, the assumption made in Eq. (22) is simplistic,but it serves for a double purpose. First, it states that L (cid:28) a .Second, it provides an order of magnitude estimate of the nu-merical resistivity. A direct inspection of Fig. 3 shows that η ∗ isbelow 10 − (2 × − ) even using p = ff usion of 1Dplasma waves (Sect. 2.1) is not straightforward due to the di ff er-ent meaning of p . For 1D plasma waves p expresses the numberof numerical zones per characteristic length of the problem (i.e.,the wavelength of the waves). Here it is the number of zones percurrent sheet width, a scale that is (significantly) larger than thecharacteristic scale of the problem at hand (i.e., L R ). Thus, theestimate of η ∗ via TMs is less accurate than it is in the 1D case.It is a lower bound to the true numerical resistivity of our methodas verified with the results of Sect. 3.2. Article number, page 11 of 14 & A proofs: manuscript no. B_code_paper
MP5MP7MP95 10 2010 - - - N o r m a l i z e d n u m e r i c a l r e s i s t i v i t y η * / ℒ Plasma waves
Rembiasz et al. ( ) ℒ / Δ z Tearing modes
Fig. 9.
Comparison of normalized numerical resistivity η ∗ / ( LV ) for se-lected cases between the presented FFE method (solid lines) and the 3DEulerian MHD code AENUS (Obergaulinger 2008) as characterized byRembiasz et al. (2017, dot-dashed lines). Left:
Comparison of the (nor-malized) numerical resistivity as derived for (smooth) plasma waves inSect. 2.1.
Right:
Contrasting approximately the development of TMs asanalyzed in Sect. 2.2.
The described situation is specific to the numerical modeling ofFFE: The dissipation of energy in current sheets proceeds nearlyinstantaneously when the force-free constraints are algebraicallyenforced (Sect. 3.3 of Paper I). Energy-momentum may leak outof the system. Contrasting this, the electromagnetic energy inresistive MHD may be stored in other forms of energy (kineticor thermal) as resistive processes develop within the system. In-deed, even when employing an ideal MHD numerical modeling,the numerical resistivity may find channels to convert magneticenergy-momentum into other dynamical components of the sys-tem (e.g., Rembiasz et al. 2017). In Fig. 9 we directly comparethe estimates of numerical resistivity determined in Sects. 2.1and 2.2, for our FFE code, with the characterization of the 3DEulerian MHD code AENUS (Obergaulinger 2008; Rembiaszet al. 2017). The contrasting results for plasma waves and tear-ing modes lucidly illustrate a key feature of FFE methods: Whilesmooth plasma waves are resolved with very competitive accu-racy, comparable to the reference MHD code, the limits of idealFFE are reached in resistive layers and discontinuities (e.g., cur-rent sheets). The presented comparison of tearing modes (rightpanel of Fig. 9) is merely an approximation due to the physicaldi ff erences between the FFE and MHD regime. The fact that ourFFE method shows less numerical di ff usion than the MHD ref-erence at small resolutions must be taken with care. Rather thanan improvement of the method, this signature highlights the factthat, in the presence of current sheets, numerical resistivity inFFE does not necessarily behave as the physical one.The physical conditions under which FFE is valid maymimic, to some extend, plasma regimes were fast (i.e., e ff ec-tive) cooling acts (cf. Komissarov 2006): The enforcement ofthe constraints of FFE is an intrinsically non-conservative pro- cess, acting almost instantaneously (see Sect. 3.3 of Paper I). Asa direct consequence, energy-momentum may leak out of the nu-merical domain over timescales as small as the simulation timestep. Di ff erently from (relativistic) MHD, the numerical resistiv-ity of an FFE code does not only depend on the discretization ofthe partial di ff erential equations, but also on the degree by whichthe FFE constraints are violated during the time evolution. Theamount by which FFE constraints are not preserved does notnecessarily decrease with increasing resolution. Thus, the nu-merical resistivity e ff ectively acts as a timescale for the numeri-cal dissipation of the electromagnetic energy-momentum, wheresmaller values of η ∗ yield shorter cooling timescales in the sensethat ideal (perfectly conducting) FFE most e ffi ciently dissipatesnon-ideal electric fields. However, the e ff ects of cooling are notspecifically accounted for in the derivation of the dispersion re-lation of TMs (in our case Eqs. 18 and 19). We may interpret thatthe standard TM dispersion relation (which strictly holds in adi-abatic , incompressible MHD) is adequate for conditions of ex-tremely slow cooling rather than in the limit of e ff ective coolingthat the FFE constraints impose. As η ∗ is predominantly dictatedby the order of the method, only FFE methods of relatively low-order may yield numerical dissipation timescales long enoughthat conditions of (su ffi ciently) slow cooling hold. However, thisfeature significantly di ff erentiates FFE from (relativistic) MHD,where η ∗ cannot be interpreted as a cooling timescale at all. j || in FFE codes The part of the force-free current which is along the magneticfield ( j || ) is closely related to the force free constraints: The cur-rent I µ ff conserves the perpendicularity of the electromagneticfields. This can be done either by combining I µ ff , i.e., the condi-tion L n ( B · D ) = ∂ t ( D · B ) =
0, withan algebraic enforcement of perpendicularity. Alternatively, thealgebraic enforcement can be replaced by suitable Ohm’s laws inthe limit of low resistivity, as the one employed by Parfrey et al.(2017), which we probed in Sect. 3.As commented above, not enforcing an instantaneous ful-filment of the FFE constraints is equivalent to introducing a fi-nite timescale for (resistive) dissipation of energy-momentum inthe system of equations. If this finite timescale is long enough,the e ff ective cooling may be relatively small during a sizablefraction of the linear phase of TM development. Hence, condi-tions of relatively slow (ine ffi cient) cooling can be phenomeno-logically mimicked in this way. However, beyond the timescalementioned above, energy-momentum yet continues to leak outof the system. This energy is not converted to other dynamicalforms, again, di ff ering from the physical mechanisms of energy-momentum transformation operating in resistive MHD.The left panel of Fig. 8 shows the e ff ect of di ff erent meth-ods for the discretization of j || on the measured TM growth rates.The solid lines correspond to our default fourth-order accuratediscretization of j || in combination with di ff erent methods of cellinterface reconstruction and a fixed number of points per cur-rent sheet width p =
6. The dashed lines correspond to tak-ing j || = (following Yu 2011). Clearly, employing j || = yields a significantly smaller growth rate and, hence a signif-icantly smaller numerical resistivity. However, employing thisprocedure, the three-wave Riemann problem renders nonphysi-cal results, as the reduced numerical di ff usivity of the methodallows for the breakup of the fast waves into a pair of discon-tinuities (one of which should not exist; see top left panel ofFig. 5). We have further probed this specific pathology in thecontext of non-ideal FFE in Sect. 3.1. We attribute the inaccurate Article number, page 12 of 14. F. Mahlmann, M. A. Aloy, V. Mewes, P. Cerdá-Durán: Computational GRFFE: Characterization of Numerical Di ff usivity modeling of charge carrying current-layers in this particular case( j || = , but enforcing the force-free constraints algebraically) tothe continuity equation of charge which we evolve in our method(Mahlmann et al. 2020a). Not supplying a current which is con-sistent with the plasma dynamics, i.e., altering the electric fieldalgebraically without feedback to the currents and charges theysupport, introduces errors into the numerical solution that willeventually lead to nonphysical behavior. For reference, we alsodisplay our default fourth-order discretization of j || (dash-dottedlines). However, in this case, we do not enforce the algebraicfulfillment of the D · B = ff erent cell interface recon-structions are indistinguishable and all the lines are located ontop of each other. From these results, we conclude that withoutresorting to phenomenological approaches to introduce resistivedissipation in the current parallel to the magnetic field, the limitsof FFE as a zeroth-order approximation of relativistic magneto-hydrodynamics are reached (and breached) in current sheets.At this point, it is important to highlight that alternative ex-pressions of the force-free current I µ ff (Eq. 27) may help to in-crease the time scale over which electromagnetic dissipation inFFE is generated. These alternative expressions of the spatialcurrent encode di ff erent forms of Ohm’s law in order to derivethe functional form of j || (see, e.g., Sect. 2.3 of Alic et al. 2012,for a comprehensive overview and also Lyutikov 2003; Li et al.2012; Parfrey et al. 2017). Such alternative Ohm’s laws relax thestrict fulfillment of the FFE constraints temporarily. Instead, theFFE constraints are asymptotically enforced with suitable driv-ing terms (see Sect. 3.3 of Paper I) or by adopting some phe-nomenological form of anisotropic electric conductivity, similarto the approach we have introduced in Sect. 3. There, we em-ploy and test the isotropic resistivity model suggested in Parfreyet al. (2017), which induces a simple Ohm’s law of the form D || → α − η J || . Several other choices of suitable Ohm’s laws areimaginable and found at least to some extend throughout the lit-erature, e.g., with di ff erent values along and across the magneticfield lines (Komissarov 2004). By employing such techniques,Komissarov et al. (2007) reproduce the theoretical dispersionrelation for TMs roughly, employing a second-order accuratemethod in combination with a (phenomenological) anisotropicconductivity. As such, it corresponds to a plasma with seeminglylarge resistivity ( η = − − − ), and even qualitatively repro-duces the linear phase of collapse of stressed current sheets ascompared with PIC simulations (Lyutikov et al. 2017). Futureupgrades of our methodology will explore the impact of di ff er-ent forms of the Ohm’s law in the numerical di ff usivity of thealgorithm more quantitatively; in Sect. 3 we have characterizedone specific prescription of resistive FFE in this context.
5. Conclusions
We have very carefully assessed the numerical resistivity of ournew GRFFE code. This topic is of the utmost importance to in-terpret the dissipation of magnetic flux and energy in magnet-ically dominated astrophysical scenarios. The intrinsic numeri-cal resistivity of our code depends on the number of numericalcells per characteristic size of the structure we aim to resolve. Wehave considered two types of tests. 1) The damping of 1D plasmawaves. Such ideal (smooth) waves allow for a direct comparisonof the damping rates with existing analytic results in the litera-ture. 2) The growth of TMs, induced by the action of the numer-ical resistivity, in a simple 2D setup. The growth rate of the TMs depends (non-linearly) on the numerical resistivity. We have per-formed an extensive suite of tests with di ff erent resolutions, spa-tial reconstruction methods and di ff erent forms of Ohm’s law(i.e., di ff erent closure relations for the current as a source of theelectromagnetic field and flux of charge conservation). With thisstrategy, we have assessed the signatures of numerical resistiv-ity in our specific method, and FFE codes in general. We haverelated these findings to the employed numerical techniques aswell as to the grid spacing employed for the discretization of theplasma continuum.Resolving a characteristic length L with at least ∼ η ∗ ∼ − ( η ∗ ∼ − ). In practical applications ofour code to, e.g., global 3D models of magnetar magnetospheres(Mahlmann et al. 2019) we employ a typical resolution of 32zones per radius of the magnetar, R ∗ . Hence, at scales of 1 R ∗ our magnetosphere models incorporate a numerical resistivity of ∼ − − − (MP5 vs. MP9 reconstruction). In the same mod-els, we observe that there exists dissipation of magnetic energyon structures with sizes as small as ∼ . R ∗ (Mahlmann et al.2019). At these scales, the numerical resistivity can be as largeas η ∗ ∼ − . Similar comments are in place for global 3D mod-els of BH magnetospheres, where we have also used 32 cells pergravitational radius in the finest grid surrounding the BH (but16 zones in the extended region surrounding the BH Mahlmannet al. 2020b). Finding dissipation on scales of the order of onegravitational radius in the vicinity of the BH (induced by the gen-eration of local current sheets) means that the numerical resistiv-ity of our code was (cid:46) − in these (small) regions of the wholecomputational domain. We highlight the non-uniform characterof numerical resistivity in practical applications. For su ffi cientlysmall scales, a GRFFE code may become extremely resistive.This justifies the implementation of very high-order spatial re-construction methods in GRFFE, even at the cost of some par-allel performance loss due to the increased number of (internal)boundary zones that need to be synchronised. As we have shown,for resolutions at reach in global 3D numerical simulations (sayof 5 −
10 cells per length scale to be resolved) using a seventh-order accurate method as MP7 reduces the numerical resistivityby more than two orders of magnitude compared to, e.g. employ-ing a second-order accurate method such as MC reconstruction.The tested resistive FFE using a simple Ohm’s law (Eq. 29,and Parfrey et al. 2017) is able to induce the expected phe-nomenological behavior in a current-carrying discontinuity(Sect. 3.1): With increasing (physical) resistivity η , the current(and associated charge density) sustaining a standing Alfvénwave becomes weaker and the layer di ff uses (Fig. 6). Eventually,the solution becomes non-physical and sharp transition layersemerge. On the other hand, the force-free solution is recoveredfor low values of η . In the case of TMs, the considered Ohm’slaw converges towards growth rates comparable to the ones ex-pected by the dynamics driven by a physical resistivity whenendowed with η above the threshold of numerical resistivity.We are, in conclusion, able to give three answers to theinitially posed questions: a) With an accurate modeling of thecurrent j (cid:107) (including a su ffi ciently high-order discretization ofit), our finite volume ideal FFE code resolves smooth current-carrying plasma waves with close-to-exact accuracy. b) In force-free current sheets a true limit of ideal FFE is reached; resistivelayers cannot be resolved on distances of the order of the grid-scale. c) When reducing the numerical resistivity by increasingresolution or the numerical accuracy (increasing the order of thespatial reconstruction), growth rates of the considered dynami-cal instability (tearing modes) in FFE are not fully comparable Article number, page 13 of 14 & A proofs: manuscript no. B_code_paper with their counterpart in numerical experiments, e.g., in resis-tive MHD. The estimated order of accuracy of the FFE methodis reduced in such regions. The action of numerical resistivity inthe limit of validity of FFE does not exactly match the action ofphysical resistivity. In other words, numerical resistivity does notnecessarily mimic physical resistivity in tearing-unstable currentsheets. Alternative Ohm’s laws that incorporate a physical resis-tivity are promising candidates for a correct physical modelingof genuinely resistive e ff ects in FFE. Thus, they may be used asan algorithmic strategy to bridge between plasma regimes in hy-brid simulation tools combining, e.g., FFE with particle in cellsimulations (Parfrey et al. 2019). The numerical techniques pre-sented in this series of publications are able to resolve the globaldynamics of force-free plasma with a great consistency of itsphysical constituents and achieves competitive accuracy.
6. Acknowledgments
We thank B. Ripperda and A. Philippov for their perspectiveson plasma (codes). We acknowledge the fruitful discussionswith T. Rembiasz, as well as his critical revision of the re-sults related to the characterization of the numerical di ff usiv-ity of our code. JM acknowledges a Ph.D. grant of the Studi-enstiftung des Deutschen Volkes . We acknowledge the supportfrom the 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