A Simflowny-based high-performance 3D code for the generalized induction equation
Daniele Viganò, David Martínez-Gómez, José A. Pons, Carlos Palenzuela, Federico Carrasco, Borja Miñano, Antoni Arbona, Carles Bona, Joan Massó
AA Simflowny-based high-performance 3D code for the generalized induction equation
Daniele Vigan`o a,b , David Mart´ınez-G´omez a,b , Jos´e A. Pons c , Carlos Palenzuela a,b , Federico Carrasco a,b ,Borja Mi˜nano a,b , Antoni Arbona a,b , Carles Bona a,b , Joan Mass´o a,b a Departament de F´ısica, Universitat de les Illes Balears and Institut d’Estudis Espacials de Catalunya, Palma de Mallorca,Baleares E-07122, Spain b Institut d’Aplicacions Computacionals de Codi Comunitari (IAC3), Universitat de les Illes Balears, Palma de Mallorca,Baleares E-07122, Spain c Departament de F´ısica Aplicada, Universitat d’Alacant, Ap. Correus 99, 03080 Alacant, Spain
Abstract
In the interior of neutron stars, the induction equation regulates the long-term evolution of the magneticfields by means of resistivity, Hall dynamics and ambipolar diffusion. Despite the apparent simplicity andcompactness of the equation, the dynamics it describes is not trivial and its understanding relies on accuratenumerical simulations. While a few works in 2D have reached a mature stage and a consensus on the generaldynamics at least for some simple initial data, only few attempts have been performed in 3D, due to thecomputational costs and the need for a proper numerical treatment of the intrinsic non-linearity of the equation.Here, we carefully analyze the general induction equation, studying its characteristic structure, and we presenta new Cartesian 3D code, generated by the user-friendly, publicly available
Simflowny platform. The codeuses high-order numerical schemes for the time and spatial discretization, and relies on the highly-scalable
SAMRAI architecture for the adaptive mesh refinement. We present the application of the code to severalbenchmark tests, showing the high order of convergence and accuracy achieved and the capabilities in termsof magnetic shock resolution and three-dimensionality. This paper paves the way for the applications to arealistic, 3D long-term evolution of neutron stars interior and, possibly, of other astrophysical sources.
Keywords:
MHD; neutron stars; Simflowny; magnetic fields; high resolution shock capturing; AMR
1. Introduction
The induction equation plays a key role in manyphysical problems involving magnetized plasma flows.In the most common and simple model for thestudy of plasmas, namely ideal magnetohydrodynam-ics (MHD), this equation arises from the combina-tion of Faraday’s law and the ideal Ohm’s law for aperfect conductor, which simply relates the electricfield to the cross product of the fluid velocity and themagnetic field. The two basic assumptions in idealMHD are that all the components of the plasma arestrongly coupled (a single fluid), and that the resis-tive timescales are much longer than the dynamicaltimescales of interest. If some of the previous as-sumptions are relaxed, more complex theoretical de-scriptions are required. For example, relaxing theassumption of negligible resistivity leads to resistiveMHD , in which a term proportional to the electriccurrent, is introduced in Ohm’s law. In addition, each constituent of the plasma mayhave a different hydrodynamical velocity. Thus,multi-fluid models with a separate set of equationsfor each component would be required to properlydescribe the dynamics of the plasma [1]. Althoughthey are more accurate and have a larger range ofapplicability, their numerical cost is significant. Nev-ertheless, the additional complexity can often be syn-thetized in simpler descriptions that combine theequations of all, or several of, the constituents of theplasma and rely on the use of a generalized Ohm’slaw.An important MHD extension of practical inter-est is that of
Hall MHD , a limit in which onecharged component (usually electrons) is free tomove, while the other charged particles (typicallyions) have a restricted mobility, or in other words,their gyro-frequency is much longer than other rel-evant timescales. Under this condition, the electro-
Preprint submitted to Computer Physics Communications December 12, 2018 a r X i v : . [ a s t r o - ph . I M ] D ec agnetic field evolves according to a non-linear in-duction equation with an electric field corrected byan extra term proportional to the Lorentz force. HallMHD finds direct applications in plasma physics andastrophysics [2]. For instance, it is used to describeopening switches in plasma laboratories, magnetic re-connections in the magneto-pause and magneto-tail(see, e.g. [3, 4, 5]), interstellar magnetic field dynam-ics, plasma expansion models, small-scale dynamos,and proto-planetary disks [6, 7, 8]. The limit in whichions are considered completely fixed, and only the dy-namics of one charged fluid (electrons) is relevant isalso known as electron MHD (eMHD). A comprehen-sive review of Hall MHD, including its characteristicmodes and the most suitable numerical methods todeal with it is described in [9, 10].A different approximation holds if we assumethat both positive and negatively charged fluids arestrongly coupled to each other, but they are weaklycoupled to the neutral components, so that the differ-ence of their velocities is smaller than their average(but not zero). In this case, there is a drift of the cou-pled charged components with respect to the neutralfluid. This is usually called ambipolar diffusion, andhas been studied in a variety of astrophysical scenar-ios, like solar physics [11, 12], star formation [13],accretion disks and molecular clouds [14, 15].Focusing on neutron stars, all the above describedeffects have been found to be important at differ-ent stages of their long-term evolution (see eg. [16]).These compact stellar objects are thought to be con-stituted by a highly (or very likely super-) conduc-tive core with a radius of about 10 km, where mat-ter is extremely packed under unique and only par-tially understood conditions of high pressure anddensity for which the nuclear forces dominate; a ∼ km-thick crust, resembling an ultra-high-density(10 − g cm − ) metal; a thin liquid and/orgaseous envelope (sometimes called ocean), and usu-ally a ∼ cm-thick atmosphere. In the crust, fullypressure-ionized electrons are free to move in a latticeof heavy ions, therefore being well represented by theeMHD limit. The evolution of magnetic fields in thecrust of neutron stars has been thoroughly studied[17, 18, 19, 20, 21, 22, 23] with more or less realis-tic 2D simulations over relevant timescales (typically,Myr), including the resistive (or Ohmic) and Hallterms and sometimes coupling temperature and mag-netic field evolution. Independently, a set of papers[24, 25, 26] have confirmed the general picture of the Hall dynamics found in the works mentioned above.Those works used a code that includes only the mag-netic field evolution and pointed out new interestingeffects. They find a Hall-attractor solution, whichcorresponds to a quasi-stationary state that slowlyevolves due to the simultaneous Hall and resistiveterms. Recently, the first 3D attempts to the wholeproblem have been presented in [27] and [28], in whichthe authors use a mixed spectral-finite difference codeto simulate a neutron star crust coupled with externalvacuum. Their results show new dynamics and thecreation and persistence of km-size magnetic struc-tures over long timescales. In the core, however,the situation is less clear and there is increasing evi-dence that ambipolar diffusion and/or fluid motionsare key to describe the magnetic field evolution. Sev-eral recent papers have started to address this prob-lem with numerical simulations or novel theoreticalideas [29, 30, 31, 32, 33, 34, 35, 36].The long-term evolution of magnetic fields, givenby the apparently simple induction equation, is anexcellent benchmark to test several theoretical issues.These are related to the interdependence with theevolving temperature, the composition and the stateof the core (which affects the relative importance ofone or another term and the dynamical timescales),the possible presence of nuclear pasta phase in thecrust/core transition layer, the chemical composi-tion of the atmosphere, the initial data (the mag-netic topology at birth is basically unknown), andthe treatment of the boundary condition at the star’ssurface. Furthermore, numerical issues, regardingthe non-linearity of the Hall and ambipolar terms,have represented a major technical obstacle to a sys-tematic study of the various scenarios. The three-dimensionality is thought to be a crucial factor tostudy realistic cases (in analogy to solar magnetism)and, ultimately, to compare the theoretical predic-tions on the observables with the observations, espe-cially for magnetars, which are young and stronglymagnetized neutron stars. Thus, the need for High-Performance-Computing (HPC) is clear. For this rea-son, in this paper we present a numerical code anda series of testbed simulations of a generalized formof the induction equation. This code was built byusing Simflowny [37, 38], a versatile platform ableto automatically generate parallelized codes for par-tial differential equations. It employs the AdaptiveMesh Refinement (AMR) libraries of
SAMRAI [39],and a graphical user interface which easily allows to2mplement equations and to choose among differenttime and space discretization schemes. The choice of
SAMRAI is mainly due to its proven high scalabil-ity [40], which is an important requirement to dealwith 3D simulations of non-linear dynamics. Thepaper is structured as follows. In § § Simflowny . In §
2. Generalized induction equation
Faraday’s law, which describes the temporal evolu-tion of the magnetic field in terms of the spatial vari-ations of the electric field, is given by (in SI units): ∂ (cid:126)B∂t = − (cid:126) ∇ × (cid:126)E. (1)In a fluid description of plasmas, Faraday’s lawis accompanied by a set of hydrodynamic equations,Amp`ere’s law, and Ohm’s law. The latter allows toclose the system of equations by relating the electricfield to other variables such as velocities and the mag-netic field. Solving this full system of equations canbe complex and, in most situations, one must assumesome simplifications. In the simplest case, when mat-ter is described as a single fluid, the electric field inthe reference frame comoving with matter is simplyrelated to the electrical conductivity, σ , and the elec-trical current density, (cid:126)j , by (cid:126)j = σ (cid:126)E, (2)which, together with Ampere’s law (neglecting thedisplacement currents) (cid:126)j = (cid:126) ∇ × (cid:126)B, (3)and Faraday’s law, results in an induction equationthat adopts the form of a vectorial diffusion equation: ∂ (cid:126)B∂t + (cid:126) ∇ × (cid:16) η (cid:126) ∇ × (cid:126)B (cid:17) = 0 , (4)where we have defined the magnetic diffusivity η = σ . In a more general case, the situation is more com-plex and different terms arise, as described in detailin [16] for neutron stars. However, we can build ageneral induction equation to study its mathematicalproperties and possible strategies for numerical so-lutions. The electric field, as any other vector, canbe expressed in any convenient basis. Let us write ageneral form as (cid:126)E = f d (cid:126)j + f h ( (cid:126)j × (cid:126)B ) − f a ( (cid:126)j × (cid:126)B ) × (cid:126)B (5)Therefore, the general vectorial equation that weaim to solve is ∂ (cid:126)B∂t + (cid:126) ∇ × (cid:16) f d (cid:126)j + f h ( (cid:126)j × (cid:126)B ) − f a ( (cid:126)j × (cid:126)B ) × (cid:126)B (cid:17) = 0(6)with the current given by eq. (3).Physically, the coefficients f d = 4 πη/c and f a aredetermined from collision terms in the Euler equa-tions and are positive-defined, while f h depends onthe charge density of the charges carrying the current:it is positive if negative charges are faster than thepositive charges (as in the eMHD limit), and negativeotherwise. In general, all these coefficients depend onthe physical conditions (density, temperature, chemi-cal composition), and one or another term may dom-inate in different regimes. Other possible terms con-tributing to the electric field are due to the thermo-electric effect [41] and the chemical potential imbal-ance [16], wherever the gradient of temperature andof chemical potential, respectively, induce a move-ment of the charges. Here we do not consider suchterms, since they require the coupled evolution withthe temperature and the chemical composition.Let us analyze each term appearing in the electricfield equation. The first term on the r.h.s. in eq. (5)is the purely resistive (Ohmic) term. The second andthird terms (the Hall term and the ambipolar diffu-sion term) can be combined in a vectorial advectionterm ( (cid:126)v × (cid:126)B ), in which (cid:126)v = − f h (cid:126)j + f a ( (cid:126)j × (cid:126)B ) is inter-preted as a weighted average velocity of the chargedfluid that implicitly depends on B .The Hall term is proportional to the current, itis quadratic in (cid:126)B , and hides a Burgers-like behav-ior, leading to the creation of localized current sheets(see e.g. [42] and test § (cid:126)B ), except for surface terms (Poynt-ing flux), or noting that the dissipated energy is (cid:126)E · (cid:126)j ∝ (cid:126)j · ( (cid:126)j × (cid:126)B ) = 0.The ambipolar diffusion term, which is cubic in (cid:126)B , is related to a velocity with the same directionof the total magnetic (Lorentz) force acting on thecharged particles, given by ( (cid:126)j × (cid:126)B ). Its effect isto dissipate the currents perpendicular to (cid:126)B , havingno effect in the current flowing along magnetic fieldlines. In other words, this term acts to align mag-netic and current fields, thus bringing the system intoa force-free configuration, characterized by definitionby (cid:126)j × (cid:126)B = 0. Note that a formally similar term hasbeen used to find configurations of twisted force-freemagnetosphere, under the name of magneto-frictionalterm [43, 44]. In this paper, we will present a simpleapplication in § The characteristic analysis captures the behavior ofperturbations around a background field at a partic-ular point in space (and time), generally related withthe study of local stability for the given evolution sys-tem and the concept of well-posedness. Moreover, itmight also provide useful information for numericalapplications, in particular how to impose boundaryconditions or to treat interfaces. With this purpose,let us consider perturbations of equation (6) over asmooth background solution (cid:126)B o . We shall look thenfor plane-wave solutions of the form, (cid:126)B = (cid:126)B o + (cid:126)B e i ( (cid:126)k · (cid:126)x − ωt ) (7)and take the high-frequency limit, in which nei-ther the amplitude (cid:126)B nor the coefficients ( f d , f h , f a )would vary in space and it has been assumed that | (cid:126)B | (cid:28) | (cid:126)B o | . Hence, ∂∂t → − iω ; (cid:126) ∇ → i(cid:126)k (8)Before continuing, let us introduce first some conve-nient notational abbreviations: A k := ˆ k · (cid:126)A ; (cid:126)A p := (cid:126)A − A k ˆ k ; (cid:126)A q = ˆ k × (cid:126)A for any vector (cid:126)A . Notice (cid:126)A p and (cid:126)A q are both or-thogonal to the propagation direction ˆ k := (cid:126)k/k (with k ≡ | (cid:126)k | ) and also perpendicular among them. Withthis notation: (cid:126) ∇ × (cid:126)B → i k (cid:126)B q ( (cid:126) ∇ × (cid:126)B ) × (cid:126)B → i k ( (cid:126)B q × (cid:126)B o ) (cid:16) ( (cid:126) ∇ × (cid:126)B ) × (cid:126)B (cid:17) × (cid:126)B → ik (cid:104) ( (cid:126)B o · (cid:126)B q ) (cid:126)B o − B o (cid:126)B q (cid:105) . Therefore, the terms in eq. (6) expand as follows: (cid:126) ∇ × (cid:126)j → k (cid:126)B p (cid:126) ∇ × (cid:16) (cid:126)j × (cid:126)B (cid:17) → − k B ok (cid:126)B q (cid:126) ∇ × (cid:104)(cid:16) (cid:126)j × (cid:126)B (cid:17) × (cid:126)B (cid:105) → − k (cid:104) ( (cid:126)B o · (cid:126)B q ) (cid:126)B oq + B o (cid:126)B p (cid:105) and the linearized system finally reads: i ωk (cid:126)B = (cid:0) f d + f a B o (cid:1) (cid:126)B p + − f h B ok (cid:126)B q + f a ( (cid:126)B o · (cid:126)B q ) (cid:126)B oq (9)Contracting the equation above with (cid:126)k reveals thatthere is a zero-mode ( ω = 0) associated with thesolenoidal constraint (cid:126) ∇ · (cid:126)B = 0, with its eigenvec-tor along the propagation direction, (cid:126)B ∝ (cid:126)k . Thismight be interpreted by saying that any deviationfrom the physical constraint (cid:126) ∇ · (cid:126)B = 0 will not prop-agate. Thus, hereafter we focus only on the physicalperturbations, which are transversal to (cid:126)k . The eigen-values are given by: i ω ± k = f d + 12 f a ( B o + B ok ) ± (cid:113) f a B op − f h B ok (10)with their respective eigenvectors , (cid:126)B ± = 2 B ok (cid:126)B op − (cid:104) f r B op ∓ (cid:113) f r B op − B ok (cid:105) (cid:126)B oq (11)where we have defined the ratio, f r := f a /f h .In the particular situation where (cid:126)B o (cid:107) (cid:126)k (i.e (cid:126)B op =0), the solution can be re-written as: (cid:126)B ± = ˆ p ± i ˆ q (12) iω ± = k (cid:8) f d + f a B ok ± if h | B ok | (cid:9) (13)being ˆ p and ˆ q any two transversal unit vectors satis-fying ˆ k × ˆ p = sign ( B ok ) ˆ q . These solutions represent Note that the particular case ( f h = 0) cannot be directlyrecovered from this expression. In this case, the eigenvectorsmust be obtained directly from eq. (9). f d = f a = 0), the system be-comes hyperbolic and the characteristic modes are nolonger damped. Namely, (cid:126)B ± = ˆ p ± i ˆ q (14) iω ± = ± ik f h | B ok | (15)Another solution for the ideal case can be obtainedby relaxing the assumption of high-frequency limitand considering a non-zero gradient of f h . In thisscenario, the so-called Hall drift wave is recovered(see e.g. [9]). This is a single mode that propagatesalong the direction ∇ f h × (cid:126)B o , with a frequency givenby ω = k |∇ f h × (cid:126)B o | . See Appendix B for a generalderivation.Finally, in the case with f a = f h = 0 and constant f d (cid:54) = 0, the characteristic modes are purely dissipa-tive and are given by any two transversal (linearly in-dependent) vectors, with a frequency iω = k f d . In-terestingly, the induction equation admits here globalsolutions such that (cid:126) ∇ × (cid:126)B = µ (cid:126)B (with µ constant),which decays in a self-similar way, i.e. (cid:126)B = (cid:126)B in e − t/τ d (16)with τ d = 4 π/f d c µ . When f d or µ varies withthe position (which is the realistic case for a neu-tron star), the dissipation is not homogeneous, thefield configuration is not self-similar, and the decaywill be enhanced in the regions with highest values of f d and (cid:126) ∇ × (cid:126)B . Previous works used numerical algorithms that pre-serve the (cid:126) ∇ · (cid:126)B = 0 constraint by construction, likethe use of staggered grids for finite difference schemes[47]. Since we employ a more general and flexible codewith a standard grid (all components of the fields aredefined and evolved at every grid node), this con-straint is not guaranteed. Therefore, we adopt a di-vergence cleaning method [48], extensively used in MHD, consisting of the addition of a scalar field asfollows: ∂ (cid:126)B∂t + (cid:126) ∇ × (cid:126)E + (cid:126) ∇ Φ = 0 ∂ Φ ∂t + c h (cid:126) ∇ · (cid:126)B = − κ Φ (17)where Φ is a scalar field that allows the propa-gation and damping of divergence errors. Noticethat c h is the propagation speed of the constraint-violating modes Φ (cid:54) = 0, which decay exponentiallyon a timescale 1 /κ . We implement this equation inour code to ensure that the constraint is always pre-served, except in the case where the treatment is notneeded because the constraint (cid:126) ∇ · (cid:126)B = 0 is numer-ically guaranteed by construction (tests where onlyone component of (cid:126)B is non-zero in § § § c h = κ . However,if κ is too large, the system of equations becomesstiff and it is difficult to evolve with explicit numer-ical schemes. We set c h = κ = 4, unless otherwisespecified (see also § C , with the magnetic energy (set-ting the magnetic permeability µ = 1 for simplicity): C = (cid:90) V || (cid:126) ∇ · (cid:126)B || (∆ x ) d V (18) E m = (cid:90) V B V (19)where, in the case of applying AMR, ∆ x is the gridsize corresponding to the finest resolution.
3. Numerical methods
The code presented here has been generated by us-ing
Simflowny [37, 38] together with the infrastruc-ture
SAMRAI [39, 40].
Simflowny is an open-sourceand user-friendly platform developed by the IAC3group since 2008 to facilitate the use of HPC infras-tructures to non-specialist scientists. It allows to eas-ily implement scientific dynamical models, by meansof a Domain Specific Language, and a web-based in-tegrated development environment, which automati-cally generates efficient parallel code for simulationframeworks.
Simflowny splits the physical models5nd problems from the numerical techniques. Theautomatic generation of the simulating code allowsto properly include the parallelization features, whichin this case rely on the
SAMRAI infrastructure [39] . SAMRAI is a patch-based structured AMR developedover more than 15 years by the Center for AppliedScientific Computing at the Lawrence Livermore Na-tional Laboratory. The latest upgrades on the AMRalgorithms allow to improve the performance andreach a good scaling on up to 1.5M cores and 2M MPItasks [40], at least for some specific problems. Thecombination of these two platforms provides a finalcode with a good balance of speed, accuracy, scala-bility, ability to switch physical models (flexibility),and the capacity to run in different infrastructures(portability).
Our system of equations can be written formally inconservation law form, namely ∂ t U + ∂ k F k ( U ) = 0 (20)where U = { (cid:126)B, Φ } is the list of evolved fields and F k ( U ) their corresponding fluxes which might benon-linear but depend only on the fields and not ontheir derivatives. Notice then that this is not com-pletely true in the eMHD system due to the presencein the fluxes of the current (cid:126)j , defined by Eq.(3).The discretization of the continuum equations isperformed by using the Method of Lines, which al-lows to address separately the time and the spacediscretization. Space finite difference schemes basedon Taylor expansions, which are extremely suitablefor smooth solutions, might be a good option to cal-culate the current, which might be computed in twodifferent ways: either by a direct finite difference dis-cretization, or invoking the Stokes theorem to inte-grate the circulation of the magnetic field on a closedpath. Interestingly, at second order accuracy bothapproaches lead to the same discretization scheme,which is our preferred choice here (see Appendix Afor a discussion about other choices), and are simplycomposed by terms like D x B yi = B yi +1 − B yi − x See also the website https://computation.llnl.gov/project/SAMRAI/ where i identifies a specific discretized cell positionon the grid along a direction perpendicular.However, these simple finite difference operatorsare not the optimal choice for the spatial discretiza-tion of fluxes in intrinsically non-linear systems likethe eMHD. In that case, it is advisable to use High-Resolution-Shock-Capturing (HRSC) methods [49],to deal with the possible appearance of shocks andto take advantage of the existence of weak solutionsin the equations. We will therefore use a conservativescheme to discretize the fluxes, which in one dimen-sion reads: ∂ t U i = − x (cid:16) F xi +1 / − F xi − / (cid:17) where F xi ± / are the set of fluxes along the x -direction evaluated at the interfaces between twoneighboring cells, located at x i ± / . The crucial issuein HRSC methods is how to approximately solve theRiemann problem, by reconstructing the fluxes at theinterfaces such that no spurious oscillations appear inthe solutions. This calculation consists in two steps: • We consider the following combination of thefluxes and the fields, at each node i : F ± i = 12 ( F i ± λU i ) (21)where λ is the maximum propagation speed ofthe system in the neighboring points. Then, wereconstruct the flux at the left (right) of eachinterface, e.g. F Li +1 / ( F Ri +1 / ), using the val-ues { F + } ( { F − } ) from the neighboring nodes { x i − n , .., x i +1+ n } . The number 2( n + 1) of suchneighbors used in the reconstruction proceduredepends on the order of the method. Simflowny already incorporates some commonly used re-constructions, like PPM [50] and MP5 [51], aswell as other implementations like the FDOCfamilies [52] which are almost as fast as cen-tered finite difference schemes at the cost of somebounded spikes near the shock regions. Herewe mostly focused on the Weighted-Essentially-Non-Oscillatory (WENO) reconstructions [53,54], which is our preferred choice for their flex-ibility (i.e., they can achieve any order of accu-racy) and robustness. The detailed implementa-tion of the WENO flavors used here can be foundin [55].6
We use a flux formula to compute the final fluxat each interface, e.g.: F i +1 / = F Li +1 / + F Ri +1 / (22)Note that our reconstruction method does notrequire the characteristic decomposition of thesystem of equations (i.e., the full spectrum ofcharacteristic velocities). At the lowest order re-construction F Li +1 / = F + i and F Ri +1 / = F − i +1 , sothat the flux formula (22) reduces to the popularand robust Local-Lax-Friedrichs flux [49].The time integration of the resulting semi-discreteequations is performed by using a 4 th -order Runge-Kutta scheme, which ensures the stability and conver-gence of the solution for a small enough time step ∆ t .We defer the reader to [55] for further details on thenumerical schemes and an extensive analysis of theperformance with different discretization schemes fordifferent problems, including MHD.We note that, in previous works, numerical codesdesigned to handle Hall dynamics show numerical sta-bility issues caused by the the quadratic dispersionrelation of the whistler waves present on the system.Recent work by [56], for instance, include stabilizingtechniques introduced in [57] for the time advance ofthe non-linear terms. These techniques, namely theSuper Time-Stepping and the Hall Diffusion Schemes,allow to maintain the stability and efficiently speedup the time evolution when the ambipolar or theHall term dominate and the equations become es-sentially a parabolic system. Another usual, butless elegant technique in eMHD to control the Hall-generated instabilities is the use of high-order dissi-pation (also called hyper-resistivity in this context)[10, 21], and/or special recipes for the time advance(e.g., a semi-implicit time advance in the toroidal fieldcomponent, which acts as an effective high-order dis-sipative term, [21]). A similar effect (i.e., a filter ofthe high-frequency modes which cannot be accuratelyresolved by our numerical grids) can be achieved byapplying artificial Kreiss-Oliger (KO) dissipation toour evolved fields { (cid:126)B, Φ } along each coordinate di-rection [58]. A sixth-order derivative dissipation op-erator, which does not spoil the convergence order ofour numerical scheme, can be written in 1D as Q xd U i = η ko x (cid:0) U i − − U i − + 15 U i − − U i + 15 U i +1 − U i +2 + U i +3 (cid:1) . where η ko is a positive, adjustable parameter control-ling the amount of dissipation added. The instabilities are especially significant when us-ing 5 th -order methods (i.e., FDOC5, WENO5 andMP5), which need a high dissipation factor to be sta-bilized, with a potential loss of accuracy. For this rea-son, in this paper we restrict ourselves only 3 rd orderschemes, that do not require any additional artificialKO dissipation. However, we do not discard usinghigher-order schemes in the future with an adequateKO dissipation. Unless otherwise indicated, hereafterwe will employ WENO3YC [59] without KO dissipa-tion.In order to test the convergence order for a givensolution u ( x, t ) (where u is a scalar or vector compo-nent), we consider the numerical discretized solutions U δ obtained with three different resolutions ∆ x δ , with δ = 0 , ,
2, such that r = ∆ x / ∆ x = ∆ x / ∆ x > r = 2). The convergence order n c of thenumerical solution at a given time can be defined as n c ( t ) = log r (cid:18) || U − U |||| U − U || (cid:19) (23)where || U m − U n || = Σ (cid:126)i | U (cid:126)im − U (cid:126)in | is the L1-norm of thedifference of the two discretized fields, and the sumis performed over all indices i k , k=1...N, identifyingthe N-dimensional grid with the lowest resolution.
4. Benchmark Tests
We now discuss the results for a battery of simula-tions performed to test individual terms of the induc-tion equation, before addressing the general problem.All simulations are run in 2D or 3D Cartesian grid, inrectangular/square or cubic dominion, respectively.The list of tests include those performed in [21] fora finite-difference 2D axisymmetric code, plus a newone designed to test the behavior of the ambipolarterm.
We first consider a smooth solution for the induc-tion equation with only the Hall term, setting f h = 1and f d = f a = 0. We consider a whistler wave inthe linear regime, discussed in § Note also that FDOC methods can be actually interpretedas a variant of KO with a local (instead of global) dissipationfactor. igure 1: Whistler wave solution for the case with 256 × t = { , , } (in units of τ ). WENO3YC, which is our standard choice, gives almostindistinguishable results. The component B z is shown in blueand red, while black arrows represent the ( j x , j y ) vector field. perturbation of a homogeneous background magneticfield (cid:126)B = B ˆ x . As shown in Appendix B.1, wecan find particular circularly polarized solutions withan amplitude | (cid:126)B | (cid:28) B , and a generic propagationvector (cid:126)k , as long as it is not aligned with (cid:126)B . Weconsider the case (cid:126)k = k (ˆ x + ˆ y ) / √
2, corresponding tothe following initial values of the magnetic field: B x = B + B cos( ky ) cos( kx ) B y = B sin( ky ) sin( kx ) (24) B z = √ B sin( ky ) cos( kx ) Figure 2: Whistler wave solution: convergence order for thetwo methods of 3 rd -order. As seen above, this solution will propagate along the x direction with velocity v w = √ f h kB (25)We set up a periodic 2D box x ∈ [0 , L ], y ∈ [ − L, L ],with L = 1, B = 1, B = 10 − , k = nπ/L , similarlyto [21]. We run a simulation following several crossingtimes. In Fig. 1 we show for the resolution (256 × t = 6 . × − . the solutions at three different times, givenin units of the Hall timescale, which is defined as τ = f h L /B . We have checked that the propagationspeed is correct, and the numerical solutions convergeto the analytical one for increasing resolutions. Fig. 2shows that both WENO3 and FDOC3 methods have3 rd -order convergence, calculated comparing the B y components in simulations with (64 × × × (cid:126) ∇· (cid:126)B = 0 constraint is maintainednumerically at the round-off error without any needfor divergence cleaning. Next, we consider a non-zero vertical gradient ofthe Hall coefficient, f h = f h (1 + βy ), with again f d = f a = 0 in the same 2D rectangular, periodicdomain, x ∈ [0 , L ], y ∈ [ − L, L ]. We set the initialvalues of the magnetic field to: B z = B + B cos( kx ) B x = B y = 0 , (26)8ith B = 1, B = 10 − , f h = 1, β = 0 . k = π/ (2 L ). This solution, as shown in § (cid:126) ∇ f h × (cid:126)B (i.e., positive x-direction in our case) withvelocity v d = βL τ = βf h B = 0 . . (27)In Fig. 3 we show results using WENO3, with L = 1.The top panel compares four different resolutions at t = 100 τ , after the wave has crossed the domain5 times (the domain crossing time is 4 L/v d = 20 τ ).The bottom panel shows the evolution at three differ-ent times with a resolution given by 256 ×
128 points,with a time-step ∆ t = 6 . × − . We checked thatthe 3 rd -order convergence is maintained. Note thatwe have kept B small to be able to follow severalperiods in the nearly linear regime, before the non-linear evolution leads to the steepening of the sinu-soidal wave and the formation of a discontinuity, in-trinsic to the Burgers-like mathematical form of thisterm. Both in this and in the next test, the (cid:126) ∇ · (cid:126)B = 0constraint is zero by construction. For the third test, also with only the Hall termactivated, we adopt the same domain and the samestratification of the Hall coefficient that we have usedin the previous test, f h ( y ) = f h (1 + βy ), and thefollowing initial conditions: B z = B cos( kx ) B x = B y = 0 (28)with k = π/L . In this case, we can write theinduction equation as a Burgers equation for B z : ∂ t B z = βB z ∂ x B z (see the detailed analysis in [21]).Any smooth solution will form discontinuities. Forour initial conditions, the expected breaking time (atwhich the solution reaches its maximum steepness be-fore starting dissipating) is τ b = λ/ (4 β ) = 2 .
5, where λ = 2 π/k = 2 is the wavelength of the solution.This test is more challenging, because of the fullnon-linearity and the formation of a current sheet af-ter the initial sinusoidal profile steepens. In Fig. 4 wecompare the different resolutions (different N shown)and we note that FDOC3 is able to follow the test,but it does not accurately reproduce the shock profile,due to the presence of a single-point spike. WENO3,instead, follows it, even though some numerical dis-sipation eventually damps the shock. Of course, the Figure 3: Hall drift wave solution. Top panel: Comparison ofdifferent resolutions with WENO3 at t = 100. Bottom panel:snapshots at three different times for the resolution 256x128.Lines in the bottom panel represent the analytical solutionswhile symbols represent the numerical solutions. All times arein units of τ . presence of a shock implies a downgrading of the or-der of convergence. The breaking time in the sim-ulations converges toward the theoretical one for in-creasing resolutions. For t = 2 (shown in the bottompanel of Fig. 4), the low- N cases have already formedand started to dissipate the shock due to the lack ofspatial resolution. We turn now to the study of purely resistive modes,by setting f d = 1, f a = f h = 0. Simulations are runin a cubic 3D Cartesian domain [ − , . We con-sider the particular case of the self-similar solution,eq. (16). In axial symmetry, the resistive eigenfunc-tions are provided by the spherical Bessel functions of9 igure 4: Hall current sheet solution. Comparison of the be-havior for FDOC3 and WENO3 at different times (top), andsolution at t = 2 for different resolutions with WENO3 (bot-tom). Times are given in units of τ . Note that the spikes inthe FDOC3 case correspond to a single-point deviation fromthe solution. index l (indicating the multipole). The initial valuesfor the dipolar l = 1 case in spherical coordinates are: B r = B ξ (cid:18) sin ξξ − cos ξ (cid:19) cos θ ,B θ = B ξ (cid:18) sin ξξ − cos ξ − ξ sin ξ (cid:19) sin θ , (29) B ϕ = B ξ µ (cid:18) sin ξξ − cos ξ (cid:19) sin θ , where ξ = µr , B is a normalization factor, andthe free parameter µ controls the amount of twistof the solution. We map this initial model to ourCartesian grid by the standard transformation, r = (cid:112) x + y + z , θ = cos − ( z/r ), ϕ = tan − ( y/x ),with the azimuthal direction contained in the ( x, y ) plane: B x = B r sin θ cos ϕ + B θ cos θ cos ϕ − B ϕ sin ϕB y = B r sin θ sin ϕ + B θ cos θ sin ϕ + B ϕ cos ϕB z = B r cos θ − B θ sin θ (30)In the limit ξ = µr →
0, we recover the solutioncorresponding to a homogeneous vertical field (cid:126)B → ( B / z .At the boundaries, we impose the exp( − t/τ d ) ana-lytical solution of eq. 16. We follow the evolution ofthe modes on several Ohmic timescales, until the fieldhas decreased to a few orders of magnitude below theinitial values. In Fig. 5, we show the evolution of the B = 1, µ = 1 case for N = 64 at three differenttimes (top), and the solution at t = τ d for differentresolutions (bottom). Despite the relatively low res-olution, the errors inherent to mapping solutions inspherical coordinates to Cartesian grids, we find anexcellent agreement with the analytical solution. Wealso found that the convergence order is about twofor 3 rd -order and 5 th -order schemes here tested. Thedifference with the nominal convergence reached inthe smooth cases above, according to our interpreta-tion, lies on the numerical errors introduced in thediscretized version of current. While for a smooth,small perturbation (whistler, Hall drift waves), theleading error is the spatial discretization, for a generalcase, where there is no background magnetic field, thesource of error is likely to come from the calculationof the current. In this respect, see Appendix A for adiscussion of different discretized expressions of cur-rent. In this case, the divergence of magnetic fieldtends to slowly grow, so that the divergence clean-ing needs to be activated to achieve the maximumaccuracy (see 5.2 for a discussion of the divergencecleaning parameters). To test the ambipolar term, we now set f a to aconstant value, and f d = f h = 0. We use a particu-lar analytical solution which is axially symmetric incylindrical coordinates with only one component ofthe magnetic field depending on the radial distanceonly ( B z ( (cid:36) ), where (cid:36) = x + y ), so that initiallythe currents are perpendicular to the magnetic field,so that ( (cid:126)j × (cid:126)B ) × (cid:126)B = − B (cid:126)j , and its evolution canbe expressed as10 igure 5: Purely resistive self-similar solutions with WENO3 ina central section plane. Evolution of the solution for N = 64 (top), and solution at t = τ d for different resolutions (bottom).Lines in the top panel represent the analytical solutions whilesymbols represent the numerical solutions. ∂B z ∂t = − f a { (cid:126) ∇ × [ B z ( (cid:126) ∇ × B z ˆ z )] } · ˆ z = f a (cid:36) ∂∂(cid:36) (cid:20) (cid:36) (cid:18) B z ∂B z ∂(cid:36) (cid:19)(cid:21) (31)This form is analogous to a non-linear diffusion equa-tion u t = (cid:126) ∇ · ( mu m − (cid:126) ∇ u ) , (32)where m is a power index. The analytical, 2D so-lutions proposed by Barenblatt and Pattle [60, 61]consist in a delta function of integral Γ at the origin,which diffuses outwards with finite velocity (the diffu-sion front is clearly defined, contrarily to the infinite-speed of the parabolic standard diffusion equation). Figure 6: Barenplatt-Pattle solution for different times for N =256 (top), and for different resolutions at t = t + 4 (bottom).The apparent tail for the low-resolution cases is a visual effect,due to the connection of the points across the front. Linescorrespond to the analytical solutions and symbols representthe numerical solutions. The analytic solution can be explicitly written as fol-lows: u ( (cid:36), t ) = max (cid:40) , t − α (cid:20) Γ − α ( m − dm (cid:36) t αm (cid:21) m − (cid:41) (33)where d is the dimension of the problem, and α =( m − /d ) − . The initial pulse spreads with afront located at (cid:36) f ( t ) = (cid:18) dmα ( m − (cid:19) t α/d (34)In our case, we test a solution for B z in the xy -plane(an expanding cylindrical flux tube), for which d = 2,11 = 3, α = 1 / f a = 3 so that the evolution isanalytically given by: B z ( x, y, t ) = t − (1 / (cid:20) Γ −
118 ( x + y ) t / (cid:21) / . (35)We set t = 1 as the beginning of our simulation,and Γ = 1 /
18, so that the front propagates accordingto ( x + y ) = t / (being x + y = 1 at t = 1).In Fig. 6 we show that our results correctly repro-duce the analytical values, both in the shape of thespreading pulse and in the correct propagation speedof the front. The sharp discontinuity in the slope of B z is well-reproduced even for low resolutions. Theconvergence order in this case is also around 2, as forthe resistive case.
5. Star-like dynamics
We now consider geometric configurations qualita-tively similar to a neutron star crust. We work in acubic domain [ − L, L ] and we set up a shell, definedby R c ≤ r ≤ R (cid:63) , where r is the distance from the cen-ter of the dominion. We set R c = 8 (the core/crust in-terface) and R (cid:63) = 10 (the star’s surface). This setupwill allow to qualitatively compare the results pro-vided by our 3D Cartesian code with the well knownHall dynamics for existing axisymmetric models in 2Dspherical coordinates [19, 20, 21, 62, 22, 24, 25, 26].Note that in a realistic set up, the values of R c and R (cid:63) depend on both the equation of state and the mass ofthe star (the larger the mass, the thinner the shell).The thickness of the crust can vary within the range ∼ − R (cid:63) ∼ −
14 km for a M = 1 . M (cid:12) star.In the shell, we assume f a = f d = 0, and f h ( r ) = [1 − ( r − δr ) /R (cid:63) ] − , (36)where δr > df h /dr in the outer crust. Thesmaller the value of δr , the larger the gradient. Weset δr = 0 .
05, which gives values of f h ∈ [1 , − ,while in the magnetosphere the dynamical timescaleis given by the speed of light (10-15 orders of magni-tude faster than the interior). This transition occursin a thin (100 m) liquid layer called the envelope. This is the reason why, in realistic simulations, onecannot numerically simulate both the short-time dy-namics and the long-term evolution. When the latteris the focus of the study, one has to place the outerboundary at densities ρ ∼ g cm − , and provideappropriate boundary conditions connecting the inte-rior field to some equilibrium state outside (usually,a vacuum or a force-free magnetosphere).We explore two simple cases for the initial mag-netic field: either a purely toroidal field, or a purelypoloidal field, both confined in the spherical shelland axially symmetric around the z -axis. These con-figurations are not realistic and are prone to dif-ferent instabilities. As a matter of fact, the mag-netic field configuration, likely relaxed to a MHD-equilibrium state, has to be constituted by an ax-isymmetric mixed poloidal-toroidal configuration, or,more realistically, by a full 3D solution, possibly con-taining small-scale structures. However, in this paperwe study them separately, since each one is a usefultest that allows us to compare them with known pre-vious simulations and to underline some importantfeatures.We test both cases with a uniform grid and witha simple mesh refinement (MR2 hereafter), set up asone additional layer with a refinement factor of two(defined as the ratio of the coarse grid sizes to the re-fined grid size). In a realistic case, we will use manylayers, in order to cover a domain with L (cid:29) R (cid:63) , andto have a fair resolution in the most dynamical re-gion. Note also that below, for these simple tests, weemploy a resolution of about ∆ x = L/N (cid:39) .
15 km.Note that the typical resolutions of the 2D sphericalcoordinates production runs of [21, 22] were given by ∼
60 angular points (arcs of ∼ . ∼ ρ ),meaning that the radial size of the cells varies inthe range ∼ −
100 m (and larger in the core).Our 3D Cartesian code can reach much better res-olutions in all directions, given the adaptability of-fered by the AMR and the efficiency of the code. In-creasing the resolution is challenging due to the non-linear Hall dynamics, which triggers short-scale, veryfast whistler waves, as we discuss below. This is thereason why we typically have to use, for these par-ticular tests, a time-step of ∆ t (cid:46) − ∆ x c , where∆ x c = 2 L/N is the spatial size of the coarse grid.12 igure 7: Evolution of the toroidal field (color scale) in the meridional slice (i.e., B y in the plane x − z ) at different times t = { , . , , . } (from top left to bottom right). The black lines represent the boundaries of the crust. The evolution of a purely toroidal field subject tothe Hall drift has been studied previously, like for ex-ample in Section 5 of [21]. We successfully reproduceda few standard cases. For conciseness, here we onlyshow the case of an initial quadrupolar field, with theopposite polarity of the example in [21], in order toexplore the behavior of the Cartesian grid when thefield becomes stronger near the magnetic axis.We set L = 15 and an initial toroidal magnetic fieldgiven by: B ϕ = B ( r − R c ) ( r − R (cid:63) ) r cos θ sin θ , (37)with the usual spherical to Cartesian transformation,so that B x = − B ϕ sin ϕ , B y = B ϕ cos ϕ , and B z = 0.Outside the shell, we impose (cid:126)B = 0. We use a number of points N = 200 in each direction, corresponding to ∼
15 points covering the thickness of the shell, andwe set c h = 4 (see next test for a discussion aboutthis parameter).In Fig. 7 we display the initial magnetic configura-tion and its evolution at different times, showing thetoroidal component in a meridional slice (i.e., B y inthe x − z plane). The expected evolution consists ina vertical displacement of both toroidal rings towardthe poles, in opposite directions. If f h is constant, thedrift is purely vertical and the toroidal field would bedisplaced until encountering the star’s surface, wherethe current forms a screening sheet at the surface orit propagates outside, depending on the treatment atsurface (see next section). However, in presence of aradial gradient of f h (as in this case), the drift veloc-ity acquires an additional negative radial component13oward the star center. Therefore, the evolution pro-ceeds such that the toroidal rings initially drift towardthe poles up to a latitude for which the different termsin the velocity drift balance, reaching a sort of equilib-rium solution. Since in this case in the core there is noevolution of the field ( f a = f h = f d = 0 for r < R c ),a radial discontinuity of the toroidal field forms atthe boundary between the crust and the core nearthe poles, sustained by a screening current flowingin the meridional direction over the r = R c surface.Most magnetic energy is confined at low magnetic co-latitude, where a stronger ring-shaped magnetic fieldforms. The magnetic field is zero in the axis, where aradial current forms. We have checked that, through-out the evolution, the integrated constraint C is keptat a level below 10 − of the integrated magnetic en-ergy E m .Last, we remind that the evolution of a purelytoroidal, axisymmetric magnetic field only under theHall term (i.e., with f a = f d = 0), notoriously impliesthat the poloidal field is kept identically zero (the op-posite is not true). In this respect, note that the useof the Cartesian grid implies that the code does notdeal with spherical components and the star’s surfaceis not exactly spherical. These features introduce aspurious poloidal field when one recovers the sphericalcomponent from the computed Cartesian ones. Wehave checked that, with this resolution, the poloidalfield is kept at a level (cid:46) (cid:126) ∇ · (cid:126)B shows its maximum. However, these deviationsdo not cause any numerical or physical instability andthey are cured by increasing the resolution, which im-proves the accuracy of the Cartesian discretization ofthe spherical surfaces. We consider now a solution with an initially purelypoloidal field, given by the following expressions for r ≤ R (cid:63) : B r = B π cos θ (cid:18) sin ξξ − cos ξξ (cid:19) ,B θ = B π sin θ (cid:18) sin ξξ − cos ξξ − sin ξξ (cid:19) , (38) B ϕ = 0 . where ξ = πr/R (cid:63) . In the center of the star, r →
0, the solution corresponding to a homogeneous field
Figure 8: Purely poloidal field: initial data in the plane perpen-dicular to the y -axis. Black lines indicate the poloidal compo-nents of the field, B x − B z . The position of the shell is indicatedby the thick line. The numerical mesh is also shown in gray,with the MR2 in the crustal region. aligned with the magnetic axis (cid:126)B → B π ˆ z . Note thatthe poloidal components correspond to the solution interms of spherical Bessel functions of test § r = R (cid:63) with a vacuum dipolar field, considered for r > R (cid:63) : B r = B cos θ R (cid:63) r ,B θ = B sin θ R (cid:63) r , (39) B ϕ = 0 . As usual, we obtain the Cartesian components by us-ing the spherical-to-Cartesian transformation shownin eqs. (30). Such initial configurations imply thepresence of toroidal currents which circulate mainlyin the core and support the dipolar magnetic field.Note that the currents are initially completely per-pendicular to the magnetic field, i.e., far from a force-free configuration.As before, we set f d = 0, keeping the same profileof f h , and we set f a > f a = 10. The higher f a , the faster the core dynamics, but the qualitative14 igure 9: Evolution at t = 0 . t = 0 . t = 1 (right) of the initially purely poloidal field case, for case (i) (top panels)and case (ii) (bottom). We show the plane perpendicular to the y -axis, the poloidal field, the shell boundaries and the numericalmesh are indicated as in Fig.8. The color scale indicates the toroidal field ( B y component). Note the kinks in the magnetic fieldlines at the star’s surface in case (i). behavior is the same. In this case, the treatment ofthe region r > R (cid:63) is important, since our poloidalfield extends outside the star. In order to comparethe dynamics, we show two representative cases of thedifferent treatments for the outer region: (i) screenedevolution , with f a ( r > R (cid:63) ) = 0, like in the toroidaltest above, and (ii) magneto-frictional coupling , set-ting f a ( r > R (cid:63) ) = 100 and a smooth transition be-tween the interior and the exterior regions given by f a ( R m < r ≤ R (cid:63) ) = 100( r − R m ) / ( R (cid:63) − R m ), with R m = 9.In both cases, the boundary of our computationaldomain is far from the star’s surface, at L = 30, whereopen boundary conditions are applied just by extrap-olating the interior solution. We have checked that,at this location, the boundary does not influence sig-nificantly the evolution within the domain (i.e., thereare no spurious waves, line distortions or additionaldissipation). In Fig. 8 we show the initial configura-tion in the meridional plane, together with the baseand the refined mesh (MR2) covering the shell.The evolution at different times (from left to right) is displayed in Fig. 9 for the two cases (case (i) inthe top row, case (ii) in the bottom row). As before,the simulation is 3D, but here we show the magneticfield lines and the toroidal field in a meridional cut,i.e. the B y component in the x − z plane.As expected, the initially poloidal dipolar field(supported by purely toroidal currents) develops asystem of poloidal currents, which gives rise to aquadrupolar toroidal component in the crust, due tothe Hall term. The corresponding toroidal field whichrises has two opposite polarities in the core/innercrust and in the outer crust. The Hall drift in thecrust drags the toroidal field close to the equator bymeans of the Burger-like dynamics analyzed above( § j ⊥ = | (cid:126)j − ( (cid:126)j · (cid:126)B ) (cid:126)B/B | . Theevolution of the j ⊥ distribution in a meridional cutis displayed in color scales in Fig. 10. At the be-15 igure 10: The same as Fig. 9, but with he white-blue-black scale indicating the squared perpendicular current, j ⊥ = | (cid:126)j − ( (cid:126)j · (cid:126)B ) (cid:126)B/B | , at t = 0 .
05 (left), t = 0 . t = 0 . t = 0 . ginning, these currents are large in the core, sinceinitially the magnetic field is purely poloidal and thecurrent field is purely toroidal (initially, j ⊥ = j ).Therefore, the ambipolar term tends to straightenthe poloidal magnetic field lines, thus reducing theunderlying toroidal current. This process (i.e., trans-form both the initially purely poloidal magnetic fieldand purely toroidal current field into mixed poloidal-toroidal configuration) is the natural way to alignthe magnetic and the current fields, to gradually ap-proach a force-free topology. Note that in the crust,where f a = 0, the perpendicular currents are not dis-sipated, consistently with what is expected.The system of currents that is gradually developedis shown in Fig. 11, for the case (ii) at t = 0 .
5. Cur-rents in the poloidal projection are symmetric withrespect to the equator, consistently with an anti-symmetric (mostly quadrupolar) toroidal magneticfield. In each quarter of the meridional plane, wecan distinguish two vortices of currents (hereafter,loops for brevity) pointing to opposite directions, cor-responding to toroidal fields of opposite polarities(blue/red in figure). Loops are separated by a layer16 igure 12: Evolution for the initially purely poloidal field case: 3d view of the magnetic field lines (blue). The orange and yellowsurfaces indicate the interfaces r = R c and r = R (cid:63) , respectively, while the purple color represents j . Top panels represent thecase without tilt of the magnetic axis and bottom panels represent the case with a tilt given by the angles θ y = θ z = π/
4. Notethat, due to technical reasons inherent to the visualization software
VisIt used here, the poloidal field lines in the two cases arenot necessarily the same. in the inner crust with no toroidal field (white). Allloops close at the equatorial plane and the magneticaxis. The current field for the case (i) is the sameas case (ii), with the only difference that the externalloops close with a very strong current sheet at thesurface, instead of extending outside the star.However, note that both cases shown here are verylikely unphysical, since a smooth matching with aforce-free solution outside the star is expected (i.e.,currents and magnetic fields are parallel to each othereverywhere to ensure (cid:126)j × (cid:126)B = 0, not like in our case(ii)), and no current sheet can be sustained at thestar’s surface (where actually the resistivity reachesits highest values, [63]). Enforcing such matchingis not trivial even in 2D polar coordinates (see thefirst simulations of this kind in [63]), unless a poten-tial (i.e., current-free) magnetospheric configurationis assumed, as usually done in 2D and 3D [21, 28].Therefore, in this test the system of currents andthe magnetic dynamics are different from those pre-vious works dealing with realistic cases [21, 28]. In those cases, the boundary condition prevent currentsfrom circulating in the outer layers of the star, givingrise to a different dynamics of the magnetic geom-etry. We defer for future works to explore furthertechniques and configurations to deal with a realisticstar/magnetosphere matching.Generally speaking, the interplay between the mag-netic fields in the different regions depends on the dif-ference in their timescales, which are set by the co-efficients of the induction equation. If the timescalesare very different at the two sides of a given interface(in contrast with what is shown in our toy model,where coefficients are of the same order across the in-terfaces), then discontinuities of the tangential mag-netic field components can rise, producing a kink inthe magnetic lines supported by current sheets flow-ing along the interface.For the same reason, at the surface, in case (i),for which f a = 0 outside the star, kinks appear in themagnetic field lines (top panels of Fig. 9) and the cur-rents close at the surface (see also the growing value of17 igure 13: Evolution of E m (normalized by its initial value)and C / E m for different choices of c h , for a non-dissipative casewith only f h in the crust. Notice that, due to the non-linearcharacter of the divergence cleaning equation, the constraintdoes not decay monotonically for c h (cid:38)
4, affecting also theevolution of E m . j ⊥ at the surface in the top panels of Fig. 10), while incase (ii) the magneto-frictional treatment employedavoids the appearance of the current sheet and allowsthe propagation of current and magnetic helicity out-side the star (colors in the bottom panels of Fig. 9),accompanied by an inflation of the poloidal field lines(note that, in this particular toy model, the config-uration outside is not force-free). We remark thatwe have explored the behavior for different profilesand magnitudes of f a and f h , but here we only showthe qualitative difference between the two mentionedcases. We refer to the next papers for a detailed dis-cussion of the magneto-frictional coupling, the inter-nal dynamics, and their physical meaning.We also performed several numerical checks forthis case. First, we run the same setups describedabove with a uniform grid with N = 400, for whichthe evolved solutions are basically indistinguishable.This ensures that, for a given desired resolution ofthe shell, we can extend further away the domain byadding layers of refined meshes, decreasing the CPUtime of the simulation by orders of magnitude com-pared to the uniform resolution case. The possibilityof extending the domain very far away while keepinga high-resolution within the star will be particularlyimportant in realistic cases, where a typical resolu-tion of O (10) m is required in the crust, followed bya smooth connection with a force-free magnetospherethat extends up to hundreds of stellar radii. The maindifferences obtained by changing the resolutions arenegligible and result mainly from the possible discon-tinuities and surface geometry at the shell boundary.Second, we check the conservation of energy andthe constraint (cid:126) ∇· (cid:126)B by evolving the same case studiedabove but only with the Hall term activated, whichis the only one strictly not dissipative (i.e., we set f a = f d = 0 everywhere, thus allowing the appear-ance of screening currents along both boundaries).The evolution of the integrated magnetic energy andthe constraint deviation C are displayed as a functionof time in Fig. 13, where we have considered five dif-ferent choices of the c h parameter (always maintain-ing c h = κ ) in the divergence cleaning equation. Themagnetic energy is maintained at a level of (cid:46) − t = 2, and the missing energy is numerically dis-sipated in the discontinuities at the equator and atthe surfaces. On the other hand, the evolution of thesolution, which is initially divergence-less, producesa rise of C . After a short transient, it approaches asaturation level of C (cid:28) − − − E m , depending18n the divergence damping timescales. As a matterof fact, the asymptotic value of C is inversely pro-portional to the value of c h , if the latter has valuesO(1). As a compromise, we set c h = 4 in our simula-tions, which is a suitable value high enough to satisfythe constraints to a good accuracy but low enoughnot to restrain further the already small time-step.Note that, for this resolution and with this choice, C (cid:28) − E m during all the simulations. The opti-mal values are problem-dependent, but we noted thatwith the chosen values of c h we obtain small devia-tions which have no impact on the evolution.Last, we have also checked that the axial symmetryis preserved in our 3D Cartesian setup: no apprecia-ble differences are seen among different meridionalcuts. For the same reason, we also consider the samesimulation with N = 200, but tilting the magneticaxis of the initial condition by different angles aroundthe z and y -axis (by appropriately applying the trans-formation to coordinates and initial magnetic fieldcomponents). Fig. 12 shows the 3D view of magneticfield lines and current intensity squared, j , at differ-ent times, for two cases: with the magnetic axis alongthe z -axis (top panels, the standard initial data de-scribed above) and with tilts around x and y axes of45 ◦ and 45 ◦ , respectively (bottom panels). The axialsymmetry is maintained in both cases. We found nosignificant differences between the two cases in thegraphic visualizations and in the integrated quanti-ties. The latter show negligible variations ( (cid:46)
2% forthe chosen resolution) for different tilt angles at agiven time. Note that such small differences can bealso attributable to the numerical integration over thedomain, which can give such differences due to thedifferent orientation of the numerical cells.This last set of tests shows that the Cartesian gridworks in the same way, with or without symmetriesaround a given direction in the initial conditions.
6. Conclusions
We have introduced a new finite-difference code tostudy a generalized 3D induction equation. Our codewas automatically generated by using
Simflowny , anopen platform able to combine AMR optimal man-agement, high-order accuracy both in space and time,modularity, scalability and user-friendliness.We have explored three terms of a general induc-tion equation: Ohmic, Hall and ambipolar, which areexpected to play a role in the magnetic and thermal evolution of different regions and epochs of a neu-tron star. We considered HRSC methods, suitablefor MHD, and explored several spatial reconstructionmethods, combined with a fourth order RK tempo-ral discretization scheme. These schemes are able tocapture both smooth and non-smooth solutions, thusbeing really suitable for the Hall dynamics.We tested the code through a set of benchmarkinitial data, showing that different versions of well-known third-order and fifth-order schemes, belong-ing to the WENO and FDOC families, reproduce theHall waves smooth analytical solutions with the ex-pected accuracy. In other tests regarding the resistiveand ambipolar terms, the convergence is limited toan order of about two, probably due to the numericaldiscretization of the curl operator used for the calcu-lation of the current. Third-order schemes seem ingeneral much more stable than fifth-order ones.The novelty of our code lies in: (i) the use ofa Cartesian 3D grid to solve at the same time allthe regions of the star, including the magnetosphere;(ii) the efficient implementation, using high-ordernumerical schemes, on the infrastructure SAMRAIwhich allows for high-scalability parallelization andAMR [38, 55]; (iii) the divergence cleaning method toensure the constraint (cid:126) ∇ · (cid:126)B = 0; (iv) the generaliza-tion of the induction equation, including all possibleterms that can act, either independently or combined.We have checked that the 3D Cartesian grid workswell in axisymmetric cases, and we performed anextended battery of checks to ensure that the well-known dynamics are correctly recovered. Also, al-though we only showed a simple illustrative case forsimplicity, we underline that our code is able to man-age an arbitrary number of dynamically refined lay-ers, with an arbitrary refinement factor. These fea-tures might allow us to extend the external domainfar away from the surface, and still resolve very accu-rately the most dynamical regions. This work pavesthe way for future 3D studies of the evolution of mag-netic fields in neutron stars. We note that, in mosttests presented here, we set the resistivity to zero,to challenge the numerical platform with an extremecase. Furthermore, we have not explored the param-eter space of initial configurations of magnetic field,which could certainly be very different from the oneshere assumed.Similarly, the crust-magnetosphere coupling hasbeen briefly faced for one case, but it will be the mainfocus of a follow-up paper. Following the ideas pro-19osed by [64], we will explore the applications of themagneto-frictional method, in order to have a smoothconnection with a magnetospheric configuration, and,notably, to study quantities like the rate of magneticenergy transferred to the magnetosphere, possibly re-lated to the observed magnetar outburst activity.In subsequent works, we aim at studying the long-term evolution for realistic profiles of the coefficients f d , f h , f a . In particular, it will be interesting toexplore the interior evolution under 3D ambipolar-dominated dynamics and the appearance of crustalsmall-scale magnetic structures, which could possiblybe associated with the hotspots commonly inferredfrom the X-ray observations of isolated neutron stars.In order to do this, we will need to couple the temper-ature evolution, i.e., solving a parabolic PDE, takinginto account the feedback between magnetic field andthermal evolution.The platform Simflowny is easily adaptable to in-clude the coupled fluid equations if needed, thus thiscode could actually be applied to other scenarioswhere Hall and ambipolar dynamics are relevant, in-cluding proto-planetary formation, molecular clouds,or the long-term magnetic evolution of other sourceslike exoplanets and white dwarfs.
Acknowledgments
We acknowledge support from the Spanish Min-istry of Economy, Industry and Competitivenessgrants AYA2016-80289-P and AYA2017-82089-ERC(AEI/FEDER, UE). CP also acknowledges supportfrom the Spanish Ministry of Education and Sci-ence through a Ramon y Cajal grant. JAP ac-knowledges the Spanish MINECO/FEDER grantAYA2015-66899-C2-2-P, the grant of Generalitat Va-lenciana PROMETEOII-2014-069. DM acknowl-edges support from Vicepresid`encia i Conselleriad’Innovaci´o, Recerca i Turisme del Govern de lesIlles Balears. The work has been done within thePHAROS COST action CA16214.
References [1] R. W. Schunk, Mathematical structure of transportequations for multispecies flows, Reviews of Geophysicsand Space Physics 15 (1977) 429–445. doi:10.1029/RG015i004p00429 .[2] E. A. Witalis, Hall magnetohydrodynamics and its appli-cations to laboratory and cosmic plasma, IEEE Transac-tions on Plasma Science 14 (1986) 842–848. doi:10.1109/TPS.1986.4316632 . [3] X. H. Deng, H. Matsumoto, Rapid magnetic reconnectionin the Earth’s magnetosphere mediated by whistler waves,Nature 410 (2001) 557–560. doi:10.1038/410557A0 .[4] F. S. Mozer, S. D. Bale, T. D. Phan, Evidence of Diffu-sion Regions at a Subsolar Magnetopause Crossing, Phys-ical Review Letters 89 (1) (2002) 015002. doi:10.1103/PhysRevLett.89.015002 .[5] C. Bard, J. C. Dorelli, On the role of system size inHall MHD magnetic reconnection, Physics of Plasmas25 (2) (2018) 022103. arXiv:1710.03612 , doi:10.1063/1.5010785 .[6] M. W. Kunz, S. A. Balbus, Ambipolar diffusion inthe magnetorotational instability, MNRAS348 (2004)355–360. arXiv:astro-ph/0309707 , doi:10.1111/j.1365-2966.2004.07383.x .[7] B. P. Pandey, M. Wardle, Hall magnetohydrodynamicsof partially ionized plasmas, MNRAS385 (2008) 2269–2278. arXiv:0707.2688 , doi:10.1111/j.1365-2966.2008.12998.x .[8] W. B´ethune, G. Lesur, J. Ferreira, Self-organisation inprotoplanetary discs. Global, non-stratified Hall-MHDsimulations, AAP589 (2016) A87. arXiv:1603.02475 , doi:10.1051/0004-6361/201527874 .[9] J. D. Huba, Theory and simulation of a high-frequencymagnetic drift wave, Physics of Fluids B 3 (1991) 3217–3225. doi:10.1063/1.859752 .[10] J. D. Huba, Hall Magnetohydrodynamics - A Tutorial, in:B¨uchner J., Dum C. & Scholer M. (Ed.), Space PlasmaSimulation, Vol. 615 of Lecture Notes in Physics, BerlinSpringer Verlag, 2003, pp. 166–192.[11] E. Khomenko, M. Collados, Heating of the MagnetizedSolar Chromosphere by Partial Ionization Effects, ApJ747(2012) 87. arXiv:1112.3374 , doi:10.1088/0004-637X/747/2/87 .[12] R. Soler, M. Carbonell, J. L. Ballester, On the Spa-tial Scales of Wave Heating in the Solar Chromosphere,ApJ810 (2015) 146. arXiv:1508.01497 , doi:10.1088/0004-637X/810/2/146 .[13] R. A. Fiedler, T. C. Mouschovias, Ambipolar diffusion andstar formation: Formation and contraction of axisymmet-ric cloud cores. I - Formulation of the problem and methodof solution, ApJ391 (1992) 199–219. doi:10.1086/171336 .[14] A. C. Jones, T. P. Downes, The Kelvin-Helmholtz insta-bility in weakly ionized plasmas: ambipolar-dominatedand Hall-dominated flows, MNRAS418 (2011) 390–400. arXiv:1107.4241 , doi:10.1111/j.1365-2966.2011.19491.x .[15] O. Gressel, N. J. Turner, R. P. Nelson, C. P. McNally,Global Simulations of Protoplanetary Disks With OhmicResistivity and Ambipolar Diffusion, ApJ801 (2015) 84. arXiv:1501.05431 , doi:10.1088/0004-637X/801/2/84 .[16] P. Goldreich, A. Reisenegger, Magnetic field decay inisolated neutron stars, ApJ395 (1992) 250–258. doi:10.1086/171646 .[17] R. Hollerbach, G. R¨udiger, The influence of Hall drift onthe magnetic fields of neutron stars, MNRAS337 (2002)216–224. arXiv:arXiv:astro-ph/0208312 , doi:10.1046/j.1365-8711.2002.05905.x .[18] R. Hollerbach, G. R¨udiger, Hall drift in the stratifiedcrusts of neutron stars, MNRAS347 (2004) 1273–1278. doi:10.1111/j.1365-2966.2004.07307.x .[19] J. A. Pons, U. Geppert, Magnetic field dissipation in eutron star crusts: from magnetars to isolated neutronstars, AAP470 (2007) 303–315. arXiv:arXiv:astro-ph/0703267 , doi:10.1051/0004-6361:20077456 .[20] J. A. Pons, J. A. Miralles, U. Geppert, Magneto-thermalevolution of neutron stars, AAP496 (2009) 207–216. arXiv:0812.3018 , doi:10.1051/0004-6361:200811229 .[21] D. Vigan`o, J. A. Pons, J. A. Miralles, A new codefor the Hall-driven magnetic evolution of neutron stars,CoPhC 183 (2012) 2042–2053. arXiv:arXiv:astro-ph/1204.4707 , doi:10.1016/j.cpc.2012.04.029 .[22] D. Vigan`o, N. Rea, J. A. Pons, R. Perna, D. N. Aguil-era, J. A. Miralles, Unifying the observational diversityof isolated neutron stars via magneto-thermal evolutionmodels, MNRAS434 (2013) 123–141. arXiv:1306.2156 , doi:10.1093/mnras/stt1008 .[23] U. Geppert, D. Vigan`o, Creation of magnetic spots atthe neutron star surface, MNRAS444 (2014) 3198–3208. arXiv:1408.3833 , doi:10.1093/mnras/stu1675 .[24] K. N. Gourgouliatos, A. Cumming, A. Reisenegger, C. Ar-maza, M. Lyutikov, J. A. Valdivia, Hall equilibria withtoroidal and poloidal fields: application to neutron stars,MNRAS434 (2013) 2480–2490. arXiv:1305.6269 , doi:10.1093/mnras/stt1195 .[25] K. N. Gourgouliatos, A. Cumming, Hall effect in neutronstar crusts: evolution, endpoint and dependence on initialconditions, MNRAS438 (2014) 1618–1629. arXiv:1311.7004 , doi:10.1093/mnras/stt2300 .[26] K. N. Gourgouliatos, A. Cumming, Hall Attractor in Ax-ially Symmetric Magnetic Fields in Neutron Star Crusts,Physical Review Letters 112 (17) (2014) 171101. arXiv:1311.7345 , doi:10.1103/PhysRevLett.112.171101 .[27] T. S. Wood, R. Hollerbach, Three Dimensional Simu-lation of the Magnetic Stress in a Neutron Star Crust,Physical Review Letters 114 (19) (2015) 191101. arXiv:1501.05149 , doi:10.1103/PhysRevLett.114.191101 .[28] K. N. Gourgouliatos, R. Hollerbach, Magnetic Axis Driftand Magnetic Spot Formation in Neutron Stars withToroidal Fields, ApJ852 (2018) 21. arXiv:1710.01338 , doi:10.3847/1538-4357/aa9d93 .[29] J. G. Elfritz, J. A. Pons, N. Rea, K. Glampedakis, D. Vi-gan`o, Simulated magnetic field expulsion in neutron starcores, MNRAS456 (2016) 4461–4474. arXiv:1512.07151 , doi:10.1093/mnras/stv2963 .[30] F. Castillo, A. Reisenegger, J. A. Valdivia, Magneticfield evolution and equilibrium configurations in neutronstar cores: the effect of ambipolar diffusion, MNRAS471(2017) 507–522. arXiv:1705.10020 , doi:10.1093/mnras/stx1604 .[31] M. E. Gusakov, E. M. Kantor, D. D. Ofengeim, Evolu-tion of the magnetic field in neutron stars, Phys. Rev.D96 (10) (2017) 103012. arXiv:1705.00508 , doi:10.1103/PhysRevD.96.103012 .[32] A. Passamonti, T. Akg¨un, J. A. Pons, J. A. Miralles, Onthe magnetic field evolution time-scale in superconductingneutron star cores, MNRAS469 (2017) 4979–4984. arXiv:1704.02016 , doi:10.1093/mnras/stx1192 .[33] A. Passamonti, T. Akg¨un, J. A. Pons, J. A. Miralles, Therelevance of ambipolar diffusion for neutron star evolution,MNRAS465 (2017) 3416–3428. arXiv:1608.00001 , doi:10.1093/mnras/stw2936 .[34] E. M. Kantor, M. E. Gusakov, A note on the am-bipolar diffusion in superfluid neutron stars, MNRAS473 (2018) 4272–4277. arXiv:1703.09216 , doi:10.1093/mnras/stx2682 .[35] A. Bransgrove, Y. Levin, A. Beloborodov, Magnetic fieldevolution of neutron stars - I. Basic formalism, numeri-cal techniques and first results, MNRAS473 (2018) 2771–2790. arXiv:1709.09167 , doi:10.1093/mnras/stx2508 .[36] D. D. Ofengeim, M. E. Gusakov, Fast magnetic field evo-lution in neutron stars: The key role of magnetically in-duced fluid motions in the core, Phys. Rev. D98 (4) (2018)043007. arXiv:1805.03956 , doi:10.1103/PhysRevD.98.043007 .[37] A. Arbona, A. Artigues, C. Bona-Casas, J. Mass´o,B. Mi˜nano, A. Rigo, M. Trias, C. Bona, Simflowny: Ageneral-purpose platform for the management of physicalmodels and simulation problems, Computer Physics Com-munications 184 (2013) 2321–2331. doi:10.1016/j.cpc.2013.04.012 .[38] A. Arbona, B. Minano, A. Rigo, B. C., C. Palenzuela,A. Artigues, C. Bona-Casas, J. Mass, Simflowny 2: Anupgraded platform for scientific modelling and simulation,CoPhC 231. doi:10.1016/j.cpc.2018.03.015 .[39] R. D. Hornung, S. R. Kohn, Managing application com-plexity in the samrai object-oriented framework, Concur-rency and Computation: Practice and Experience 14 (5)(2002) 347–368. doi:10.1002/cpe.652 .URL http://dx.doi.org/10.1002/cpe.652 [40] B. T. Gunney, R. W. Anderson, Advances in patch-based adaptive mesh refinement scalability, Journal ofParallel and Distributed Computing 89 (2016) 65 – 84. doi:https://doi.org/10.1016/j.jpdc.2015.11.005 .URL [41] U. Geppert, H.-J. Wiebicke, Amplification of neutron starmagnetic fields by thermoelectric effects. I - General for-malism, AAPS 87 (1991) 217–228.[42] S. I. Vainshtein, S. M. Chitre, A. V. Olinto, Rapid dissipa-tion of magnetic fields due to the Hall current, Phys. Rev.E61 (2000) 4422–4430. arXiv:arXiv:astro-ph/9911386 , doi:10.1103/PhysRevE.61.4422 .[43] G. Roumeliotis, P. A. Sturrock, S. K. Antiochos, A Nu-merical Study of the Sudden Eruption of Sheared Mag-netic Fields, ApJ423 (1994) 847–+. doi:10.1086/173862 .[44] D. Vigan`o, J. A. Pons, J. A. Miralles, Force-free twistedmagnetospheres of neutron stars, AAP533 (2011) A125. arXiv:1106.5934 , doi:10.1051/0004-6361/201117105 .[45] R. A. Helliwell, Whistlers and Related Ionospheric Phe-nomena, 1965.[46] D. Nunn, A self-consistent theory of triggered VLFemissions, Planss 22 (1974) 349–378. doi:10.1016/0032-0633(74)90070-1 .[47] G. T´oth, Journal of Computational Physics 161 (2000)605–652. doi:10.1006/jcph.2000.6519 .[48] A. Dedner, F. Kemm, D. Kr¨oner, C.-D. Munz,T. Schnitzer, M. Wesenberg, Hyperbolic DivergenceCleaning for the MHD Equations, Journal of Computa-tional Physics 175 (2002) 645–673. doi:10.1006/jcph.2001.6961 .[49] E. Toro, Riemann Solvers and Numerical Methods forFluid Dynamics: A Practical Introduction, Springer,1997.URL https://books.google.es/books?id=6QFAAQAAIAAJ
50] P. Colella, P. R. Woodward, The Piecewise ParabolicMethod (PPM) for Gas-Dynamical Simulations, Jour-nal of Computational Physics 54 (1984) 174–201. doi:10.1016/0021-9991(84)90143-8 .[51] A. Suresh, H. Huynh, Accurate monotonicity-preservingschemes with rungekutta time stepping, Journalof Computational Physics 136 (1) (1997) 83 – 99. doi:https://doi.org/10.1006/jcph.1997.5745 .URL [52] C. Bona, C. Bona-Casas, J. Terradas, Linear high-resolution schemes for hyperbolic conservation laws: TVBnumerical evidence, Journal of Computational Physics228 (2009) 2266–2281. arXiv:0810.2185 , doi:10.1016/j.jcp.2008.12.010 .[53] G.-S. Jiang, C.-W. Shu, Efficient implementa-tion of weighted eno schemes, Journal of Com-putational Physics 126 (1) (1996) 202 – 228. doi:https://doi.org/10.1006/jcph.1996.0130 .URL [54] C.-W. Shu, Essentially non-oscillatory and weighted essen-tially non-oscillatory schemes for hyperbolic conservationlaws, Springer Berlin Heidelberg, Berlin, Heidelberg, 1998,pp. 325–432. doi:10.1007/BFb0096355 .URL https://doi.org/10.1007/BFb0096355 [55] C. Palenzuela, B. Mi˜nano, D. Vigan`o, A. Arbona,C. Bona-Casas, A. Rigo, M. Bezares, C. Bona,J. Mass´o, A Simflowny-based finite-difference code forhigh-performance computing in numerical relativity, Clas-sical and Quantum Gravity 35 (18) (2018) 185007. arXiv:1806.04182 , doi:10.1088/1361-6382/aad7f6 .[56] P. A. Gonz´alez-Morales, E. Khomenko, T. P. Downes,A. de Vicente, MHDSTS: a new explicit numericalscheme for simulations of partially ionised solar plasma,AAP615 (2018) A67. arXiv:1803.04891 , doi:10.1051/0004-6361/201731916 .[57] S. O’Sullivan, T. P. Downes, An explicit scheme formultifluid magnetohydrodynamics, MNRAS366 (2006)1329–1336. arXiv:astro-ph/0511478 , doi:10.1111/j.1365-2966.2005.09898.x .[58] G. Calabrese, L. Lehner, O. Reula, O. Sarbach, M. Tiglio,Summation by parts and dissipation for domains withexcised regions, Classical and Quantum Gravity 21(2004) 5735–5757. arXiv:gr-qc/0308007 , doi:10.1088/0264-9381/21/24/004 .[59] N. K. Yamaleev, M. H. Carpenter, Third-orderenergy stable weno scheme, Journal of Compu-tational Physics 228 (8) (2009) 3025 – 3047. doi:https://doi.org/10.1016/j.jcp.2009.01.011 .URL [60] G. I. Barenblatt, On some unsteady fluid and gas mo-tions in a porous medium., Prikladnaya Matematika iMekhanika 16 (1).[61] R. E. Pattle, DIFFUSION FROM AN IN-STANTANEOUS POINT SOURCE WITH ACONCENTRATION-DEPENDENT COEFFICIENT,The Quarterly Journal of Mechanics and AppliedMathematics 12 (4).[62] D. Vigan`o, J. A. Pons, Central compact objects and thehidden magnetic field scenario, MNRAS425 (2012) 2487– 2492. arXiv:1206.2014 , doi:10.1111/j.1365-2966.2012.21679.x .[63] T. Akg¨un, P. Cerd´a-Dur´an, J. A. Miralles, J. A. Pons,Crust-magnetosphere coupling during magnetar evolu-tion and implications for the surface temperature, MN-RAS481 (2018) 5331–5338. arXiv:1807.09021 , doi:10.1093/mnras/sty2669 .[64] T. Akg¨un, P. Cerd´a-Dur´an, J. A. Miralles, J. A. Pons, Theforce-free twisted magnetosphere of a neutron star - II. De-generacies of the Grad-Shafranov equation, MNRAS474(2018) 625–635. arXiv:1709.10044 , doi:10.1093/mnras/stx2814 . Appendix A. Discretized current
We tried with different prescription of the current.Note that its definition, eq. (3), can be discretized intwo ways: either by a direct finite difference in thetwo direction, or applying the Stokes theorem andgiving some prescription for the circulation C kl of themagnetic field (where ( k, l ) is the plane perpendicularto the current component J m that we are calculating,and i, j identify their correspondent discretized posi-tion along the two directions). The circulation canbe chosen by considering the line integral along thesquare centered in the point i, j , with sides given by2∆ x i , 2∆ x j , i.e., the square delimited by the 8 clos-est points of the grid in the plane: the values at thecenter and corners of the four sides. In general, theline integral can be evaluated as: J m = (cid:15) klm C kl (A.1) C kl = { w s [ B l ( i + 1 , j − − B l ( i − , j − w c [ B l ( i + 1 , j ) − B l ( i − , j )]+ w s [ B l ( i + 1 , j + 1) − B l ( i + 1 , j − }−{ w s [ B k ( i − , j + 1) − B k ( i − , j − w c [ B k ( i, j + 1) − B k ( i, j − w s [ B k ( i + 1 , j + 1) − B k ( i + 1 , j − } w s and w c the non-negative weights to the side andcentral neighbors, respectively, so that 2 w s + w c = 1.We tested different combination of weights within therange w c ∈ [0 , w c = 1 , w s = 0, which is actually the purely finitedifference scheme with no correction coming from thecorner neighbors. This can be explained by the factthat the inclusion of the values of the corner means22hat we are including a correction in one direction(let us say, k ), when evaluating the finite differencealong the perpendicular direction ( l ). In this sense,the form of the curl operator is intrinsically differentfrom a n th -order derivative operator acting in one di-rection only: for the latter, extending the stencil inthe direction in a proper way improves the accuracy.We also note that enlarging the number of stencilpoints for the finite difference version ( w c = 1) makesthe code very unstable.A potentially interesting alternative is to evolve thecurrent density (cid:126)J independently, instead of definingit by eq. (3). By taking the curl to eq. 6, we obtainthe following equation in the non-conservative form ∂ (cid:126)J∂t − (cid:126) ∇ ( (cid:126) ∇ · (cid:126)E ) = −∇ (cid:126)E (A.2)where (cid:126)E is given by eq. 5. The flux term con-tains first-order derivatives, and the source term hassecond-order derivatives of the variables. This makesthe resulting code very unstable with our HRSCmethods, thus making this formulation not apt forthe numerical techniques here employed. Appendix B. Some extensions to the charac-teristic analysis of the system
Appendix B.1. Particular solutions in Cartesian co-ordinates
If we particularize the general characteristic solu-tions (10)-(11) to a Cartesian coordinate system, byassuming B ok (cid:54) = 0, we may chose: (cid:126)k = k √ x + ˆ y ) , (cid:126)B o = B o ˆ x (B.1)(i.e. an angle 45 ◦ between (cid:126)k and (cid:126)B o ) . A straightfor-ward calculation then yields, (cid:126)B ± = 2 f h (ˆ x − ˆ y ) + (cid:20) f a B o ∓ (cid:113) f a B o − f h (cid:21) ˆ ziω ± = k (cid:26) f d + 34 f a B o ± B o (cid:113) f a B o − f h (cid:27) which in the ideal case, f d = f a = 0 (selecting the“+” wave-mode), we can re-write as: (cid:126)B = B o ˆ x + B (cid:104) ˆ x − ˆ y + i √ z (cid:105) e i (cid:16) k √ ( x + y ) − ωt (cid:17) Notice there is no loss of generality by adopting this par-ticular choice and it could be generalized easily, as long as (cid:126)k · (cid:126)B o (cid:54) = 0. where ω = f h k B o / √
2. This is equivalent to theinitial data of the whistler test eq. (24), by setting t = 0, k x ≡ k/ √ v ω ≡ ω/k x , and taking just the realpart. The solution consists of a circularly polarizedwave in the plane perpendicular to ˆ x + ˆ y .Another interesting solution can be found by look-ing at the case (cid:126)k (cid:107) (cid:126)B o (i.e. (cid:126)B op = 0). Thus, (cid:126)k = k ˆ x , (cid:126)B o = B o ˆ x (B.2)leading to the damped whistler wave described by, (cid:126)B = B o ˆ x + B [ˆ y + i ˆ z ] e i ( kx − ω h t ) e − t/τ d with ω h = f h k B ok being the whistler frequency ofthe propagation and 1 /τ d = k (cid:0) f d + f a B o (cid:1) control-ling the parabolic decay. Appendix B.2. Non-homogeneous coefficients
If we relax the assumption of high-frequency limit,we can consider more general solutions to the lin-earized problem by studying equations with non-negligible gradients of the coefficients: f d = f d ( (cid:126)x ), f a = f a ( (cid:126)x ) or f h = f h ( (cid:126)x ).Let us define the gradient as, (cid:126)H := (cid:126) ∇ f h . Then, thegeneral equation we need to solve is: iω (cid:126)B = ω d (cid:126)B p − ω h (cid:126)B q + f a k ( (cid:126)B o · (cid:126)B q ) (cid:126)B oq + ik (cid:104) ( (cid:126)H · (cid:126)B o ) (cid:126)B q − ( (cid:126)H · (cid:126)B q ) (cid:126)B o (cid:105) (B.3)where we have used again that ω h = f h k B ok and ω d = k (cid:0) f d + f a B o (cid:1) . This leads to a complicateddispersion relation for the physical modes of the sys-tem, namely:( X + ikH oq ) (cid:0) X + f a k B op (cid:1) + Y ( Y − ikH op ) = 0where X := iω − ω d , Y := ω h − ikH k B ok and H op/q ≡ (cid:126)H · (cid:126)B op/q . Note it is a quadratic system for X, butwith complex coefficients. To our knowledge, there isno trivial generic solution.We shall then particularize the above system tosome interesting cases within the Hall MHD approx-imation ( f d = f a = 0), setting the propagation di-rection orthogonal to both (cid:126)B o and (cid:126)H . In this casethe dispersion relation reduces to ω ( ω + kH oq ) = 0,thus recovering the so-called Hall drift mode (see [9]),which can be written, (cid:126)B = (cid:126)H q , ω = − kH oq (B.4)It might be relevant to note that within Hall MHDthere are no allowed solutions that propagate alongthe f h gradient, i.e. (cid:126)k (cid:107) (cid:126)H(cid:126)H