Magnetic, thermal and rotational evolution of isolated neutron stars

The strong magnetic field of neutron stars is intimately coupled to the observed temperature and spectral properties, as well as to the observed timing properties (distribution of spin periods and period derivatives). Thus, a proper theoretical and numerical study of the magnetic field evolution equations, supplemented with detailed calculations of microphysical properties (heat and electrical conductivity, neutrino emission rates) is crucial to understand how the strength and topology of the magnetic field vary as a function of age, which in turn is the key to decipher the physical processes behind the varied neutron star phenomenology. In this review, we go through the basic theory describing the magneto-thermal evolution models of neutron stars, focusing on numerical techniques, and providing a battery of benchmark tests to be used as a reference for present and future code developments. We summarize well-known results from axisymmetric cases, give a new look at the latest 3D advances, and present an overview of the expectations for the field in the coming years.


Introduction
Neutron stars (NSs), the endpoints of the evolution of massive stars, are fascinating astrophysical sources that display a bewildering variety of manifestations. They are arguably the only stable environment in the present Universe where extreme physical conditions of density, temperature, gravity, and magnetic fields, are realized simultaneously. Thus, they are ideal laboratories to study the properties of matter and the surrounding plasma under such extreme limits. NSs were first discovered as rotationpowered radio pulsars (standing for pulsating stars, due to their periodic signal), sometimes called standard pulsars, the most numerous class with about three thousand identified members. 1 The number is continuously increasing thanks to new extended surveys and the use of high-sensitivity instruments like LOFAR (van Haarlem 2013), and a few more thousand sources are expected to be observed by the soonavailable Square Kilometre Array. To a lesser extent, NSs have also been observed in X rays (about one hundred NSs so far), as persistent or transient sources, and/or as γ-ray pulsars (over two hundred and fifty so far). In most cases, this high-energy radiation is non-thermal, originated by particle acceleration (synchro-curvature emission, Zhang and Cheng 1997;Viganò et al 2015) or Compton up-scattering of lower-energy photons by the particles composing the magnetospheric plasma (Lyutikov and Gavriil 2006).
A particularly intriguing class of isolated NSs are the magnetars (Mereghetti et al 2015;Turolla et al 2015;Kaspi and Beloborodov 2017), relatively slow rotators with typical spin periods of several seconds and ultra-strong magnetic fields (10 13 -10 15 G). In most cases, they show a relatively high persistent (i.e., constant over many years) X-ray luminosity (L x ≈ 10 33 -10 35 erg/s), well exceeding their rotational energy losses, in contrast with radio (standard) and γ-ray pulsars. This leads to the conclusion that the main source of energy is provided by the strong magnetic field, instead of rotational energy, in agreement with the high values of the surface dipolar magnetic field inferred from the timing properties. Magnetars are also identified for their complex transient phenomenology in high energy X-rays and γ-rays, including short (tenths of a second) bursts, occasional energetic outbursts with months-long afterglows (Rea and Esposito 2011;Coti Zelati et al 2018) and, much more rarely (only three observed so far), giant flares (Hurley et al 1999;Palmer et al 2005). During giant flares, the energy release is as large as 10 46 erg in less than a second. The source of energy of such transient, violent behavior is also generally agreed to be of magnetic origin, as proposed in Thompson andDuncan 1995, 1996. Alternative or complementary power sources, such as accretion, nuclear reactions, or residual cooling from the interior, are less effective to account for the transient activity.
Although isolated NSs have been historically differentiated in sub-classes, mostly based on observational grounds (detectability in X and/or radio, transient vs. persistent properties, and presence/absence of pulsations), there is no sharp boundary between classes, and the distributions of their physical properties, such as the inferred magnetic field, partially overlap. Indeed, the evidence accumulated in the last decade has shown that the presence of a strong dipolar field is not a sufficient condition to trigger observable magnetar-like events and, conversely, there has been an increasing number of low-magnetic-field magnetars discovered in the recent past (Rea et al 2010(Rea et al , 2012(Rea et al , 2014. They are NSs with relatively low values of the inferred surface dipolar magnetic fields, showing nevertheless magnetar-like activity. Similar activity has been displayed by a couple of high-magnetic-field radio pulsars with inferred B ∼ 10 13 G (Gavriil et al 2008;Gögüs et al 2016), and by a puzzling young, extremely slowly spinning NS , belonging to the so-called sub-class of central compact objects, a handful of young NSs surrounded by a supernova rem-nant, detectable due to a persistent, mostly non-pulsating X-ray emission (De Luca 2017).
It seems now clear that the non-linear, dynamical interplay between the internal and external magnetic field evolution plays a key role to understand the observed phenomenology, and their study requires numerical simulations. Particularly important issues are the transfer of energy between toroidal and poloidal components and between different scales, the location and distribution of long-lived electrical currents within the star, how magnetic helicity can be generated and transferred to the exterior to sustain magnetospheric currents (i.e., how to twist the magnetic field lines), and how instabilities leading to outbursts and flares are triggered. In order to answer all these questions, 2D and 3D numerical simulations are required. The problem is similar to other scenarios in plasma physics or solar physics, but with extreme conditions and additional ingredients (strong gravity and possibly superconductivity). The goal of this paper is to provide an overview of the subject of modeling NS evolution accessible not only to specialists on the subject, but to a wider community including astrophysicists in general, and particularly students. For this purpose, we will review the basic equations and the numerical techniques applied to each part of the problem, with a special focus on the distinctive features of NSs, compared to other stellar sources.
This work is organized as follows. In Sect. 2 the theory of the cooling of NSs is reviewed; the magnetic field evolution is described in detail in Sect. 3, where we discuss the physical processes in different parts of the star. In Sect. 4 we review the specific numerical methods and techniques used to model the magnetic evolution. They can be implemented and tested with the benchmark cases presented in Sect. 5. In Sect. 6 we discuss the challenging coupling between the slowly evolving interior and the force-free magnetosphere, and how it determines the evolution of the spin period. Some examples of realistic evolution models from the recent literature are presented in Sect. 7. Finally, in Sect. 8 we comment on future developments and open issues.

Neutron star cooling
For a few tens of isolated NSs, the detected X-ray spectra show a clear thermal contribution directly originated from a relatively large fraction of the star surface. For the cases in which an independent estimate of the star age is also available, one can study how temperatures correlate with age, which turns out to be an indirect method to test the physics of the NS interior. The evolution of the temperature in a NS was theoretically explored even before the first detections, in the 1960s (Tsuruta 1964). Today, NS cooling is the most widely accepted terminology for the research area studying how temperature evolves as NSs age and their observable effects. We refer the interested reader to the introduction in a recent review (Potekhin et al 2015b) for a thorough historical overview of the foundations of the NS cooling theory.
According to the standard theory, a proto-NS is born as extremely hot and liquid, with T 10 10 K, and a relatively large radius, ∼ 100 km. Within a minute, it becomes transparent to neutrinos and shrinks to its final size, R ∼ 12 km (Burrows and Lattimer 1986;Keil and Janka 1995;Pons et al 1999). Neutrino transparency marks the starting point of the long-term cooling. At the initially high temperatures, there is a copious production of thermal neutrinos that abandon the NS core draining energy from the interior. In a few minutes, the temperature drops by another order of magnitude to T ∼ 10 9 K, below the melting point of a layer where matter begins to crystallize, forming the crust. Since the melting temperature depends on the local value of density, the gradual growth of the crust takes place from hours to months after birth. The outermost layer (the envelope, sometimes called the ocean) with a typical thickness O(10 2 m), remains liquid and possibly wrapped by a very thin O(cm) gaseous atmosphere. In the inner core, a mix of neutrons, electrons, protons and plausibly more exotic particles (muons, hyperons, or even deconfined quark matter), the thermal conductivity is so large that the dense core quickly becomes isothermal.
The central idea of NS cooling studies is to produce realistic evolution models that, when confronted with observations of the thermal emission of NSs with different ages Yakovlev and Pethick 2004;Yakovlev et al 2008;Page 2009;Tsuruta 2009;Potekhin et al 2015a), provide useful information about the chemical composition, the magnetic field strength and topology of the regions where this radiation is produced, or even the properties of matter at higher densities deeper inside the star. Two interesting examples are the low temperature (and thermal luminosity) shown by the Vela pulsar, arguably a piece of evidence for fast neutrino emission associated to higher central densities or exotic matter, or the controversial observational evidence for fast cooling of the supernova remnant in Cassiopeia A (Heinke and Ho 2010;Posselt and Pavlov 2018), proposed to be a signature of the core undergoing a superfluid transition (Page et al 2011;Shternin et al 2011;Ho et al 2015;Wijngaarden et al 2019). We now review the theory of NS cooling, beginning with a brief revision of the stellar structure equations and by introducing notation.

Neutron star structure
The first NS cooling studies (and most of the recent works too) considered a spherically symmetric 1D background star, in part for simplicity, and in part motivated by the small deviations expected. The matter distribution can be assumed to be spherically symmetric to a very good approximation, except for the extreme (unobserved) cases of structural deformations due to spin values close to the breakup values (P 1 ms) or ultra-strong magnetic fields (B 10 18 G, unlikely to be realized in nature). Therefore, using spherical coordinates (r, θ , ϕ), the space-time structure is accurately described by the Schwarzschild metric ds 2 = −e 2ν(r) c 2 dt 2 + e 2λ (r) dr 2 + r 2 (dθ 2 + sin 2 θ dϕ 2 ), where λ (r) = − 1 2 ln 1 − 2G c 2 m(r) r 2 accounts for the space-time curvature, m(r) = 4π r 0 ρ(r)r 2 dr is the gravitational mass inside a sphere of radius r, ρ is the mass-energy density, G is the gravitational constant, and c is the speed of light. The lapse function e 2ν(r) is determined by the equation with the boundary condition e 2ν(R) = 1 − 2GM/c 2 R at the stellar radius r = R. Here, M ≡ m(R) is the total gravitational mass of the star. The pressure profile, P(r), is determined by the Tolman-Oppenheimer-Volkoff equation Throughout the text, we will keep track of the metric factors for consistency, unless indicated. The Newtonian limit can easily be recovered by setting e ν = e λ = 1 in all equations. To close the system of equations, one must provide the equation of state (EoS), i.e., the dependence of the pressure on the other variables P = P(ρ, T,Y i ) (Y i indicating the particle fraction of each species). Since the Fermi energy of all particles is much higher than the thermal energy (except in the outermost layers) the dominant contribution is given by degeneracy pressure. The thermal and magnetic contributions to the pressure, for typical conditions, are negligible in most of the star volume. Besides, the assumptions of charge neutrality and beta-equilibrium uniquely determine the composition at a given density. Thus, one can assume an effective barotropic EoS, P = P(ρ), to calculate the background mechanical structure. Therefore, the radial profiles describing the energy-mass density and chemical composition can be calculated once and kept fixed as a background star model for the thermal evolution simulations.
In Fig. 1 we show a typical profile of a NS, obtained with the EoS SLy4 (Douchin and Haensel 2001), which is among the realistic EoS supporting a maximum mass compatible with the observations, M max ∼ 2.0-2. . We show the enclosed radius and mass, and the fractions of the different components, as a function of density, from the outer crust to the core. For densities ρ 4 × 10 11 g cm −3 , neutrons drip out the nuclei and, for low enough temperatures, they would become superfluid. Note that the core contains about 99% of the mass and comprises 70-90% of the star volume (depending on the total mass and EoS). Envelope and atmosphere are not represented here. For a more detailed discussion we refer to, e.g., Haensel et al (2007); Potekhin et al (2015b).

Heat transfer equation
Spherical symmetry was also assumed in most NS cooling studies during the 1980s and 1990s. However, in the 21st century, the unprecedented amount of data collected by soft X-ray observatories such as Chandra and XMM-Newton, provided evidence that most nearby NSs whose thermal emission is visible in the X-ray band of the electromagnetic spectrum show some anisotropic temperature distribution (Haberl 2007;Fig. 1 Structure and composition of a 1.4 M NS, with SLy EoS. The plot shows, as a function of density from the outer crust to the core, the following quantities: mass fraction in the form of nuclei X h (blue dot-dashed line), the fraction of electrons per baryon Y e (black dashes), the fraction of free neutrons per baryon Y n (red dashes), the atomic number Z (dark green triple dot-dashed), the mass number A (cyan long dashes), radius normalized to R (pink solid), and the corresponding enclosed mass normalized to the star mass (green solid). Posselt et al 2007;Kaplan et al 2011) . This observational evidence made clear the need to build multi-dimensional models and gave a new impulse to the development of the cooling theory including 2D effects (Geppert et al 2004(Geppert et al , 2006Page et al 2007;Aguilera et al 2008b,a;Viganò et al 2013). The cooling theory builds upon the heat transfer equation, which includes both flux transport and source/sink terms.
The equation governing the temperature evolution at each point of the star's interior reads: where c v is specific heat, and the source term is given by the neutrino emissivity Q (accounting for energy losses by neutrino emission), and the heating power per unit volume H, both functions of temperature, in general. The latter can include contributions from accretion and, more relevant for this paper, Joule heating by magnetic field dissipation. All these quantities (including the temperature) vary in space and are measured in the local frame, with the metric (redshift) corrections accounting for the change to the observer's frame at infinity. 2 The heat flux density F is given by withκ being the thermal conductivity tensor. In Fig. 2 we show the different contributions to the specific heat by ions, electrons, protons and neutrons, for T = {10, 5, 1, 0.5}× 10 8 K, respectively, computed again with SLy EoS. For the superfluid/superconducting 2 Throughout the text, we will use the ∇ operator for conciseness, but we note that it must include the metric factors, e.g., using the metric (1), the gradient would be ∇ ≡ e −λ ∂ ∂ r , 1 r ∂ ∂ θ , 1 r sin θ ∂ ∂ ϕ .

Fig. 2
Contributions to the specific heat from neutrons (red dashes), protons (green dot-dashed), electrons (blue dots), and ions (black solid line) as a function of density, from the outer crust to the core, and for different temperatures in each panel (as indicated). The superfluid gaps employed are the same as in (Ho et al 2012).
gaps we use the phenomenological formula for the momentum dependence of the energy gap at zero temperature employed in Ho et al (2012), in particular their deep neutron triplet model. The bulk of the total heat capacity of a NS is given by the core, where most of the mass is contained. The regions with superfluid nucleons are visible as deep drops of the specific heat. The proton contribution is always negligible. Neutrons in the outer core are not superfluid, thus their contribution is dominant. The crustal specific heat is given by the dripped neutrons, the degenerate electron gas and the nuclear lattice (van Riper 1991). The specific heat of the lattice is generally the main contribution, except in parts of the inner crust where neutrons are not superfluid, or for temperatures 10 8 K, when the electron contribution becomes dominant. In any case, the small volume of the crust implies that its heat capacity is small in comparison to the core contribution. For a detailed computation of the specific heat and other transport properties, we recommend the codes publicly available at http://www.ioffe.ru/astro/EIP/, describing the EoS for a strongly magnetized, fully ionized electron-ion plasma (Potekhin and Chabrier 2010).
The second ingredient needed to solve the heat transfer equation is the thermal conductivity (dominated by electrons, due to their larger mobility). For weak magnetic fields, the conductivity is isotropic: the tensor becomes a scalar quantity times the identity matrix. Since the background is spherically symmetric, at first approxi-mation, the temperature gradients are essentially radial throughout most of the star. In this limit, 1D models are accurately representing reality, at least in the core and inner crust. However, for strong magnetic fields (needed to model magnetars), the electron thermal conductivity tensor becomes anisotropic also in the crust: in the direction perpendicular to the magnetic field the conductivity is strongly suppressed, which reduces the heat flow orthogonal to the magnetic field lines.
In the relaxation time approximation, the ratio of conductivities parallel (κ ) and orthogonal (κ ⊥ ) to the magnetic field is Here we have introduced the so-called magnetization parameter (Urpin and Yakovlev 1980), ω B τ e , where τ e is the electron relaxation time and ω B = eB/m * e c is the gyrofrequency of electrons with charge −e and effective mass m * e moving in a magnetic field with intensity B. Equation (6) is only strictly valid in the classical approximation (see Potekhin and Chabrier 2018 for a recent discussion of quantizing effects), but this dimensionless quantity is always a good indicator of the suppression of the thermal conductivity in the transverse direction. We will see later that this is also the relevant parameter to discriminate between different regimes for the magnetic field evolution. Figure 3 shows the thermal conductivity including the contributions of all relevant carriers, for two different combinations of temperatures and magnetic field, roughly corresponding to a recently born magnetar (T = 10 9 K, B = 10 15 G), or after ∼ 10 4 yr (T = 10 8 K, B = 10 14 G). Note that the thermal conductivity of the core is several orders of magnitude higher than in the crust, which results in a nearly isothermal core. Thus, the precise value of the core thermal conductivity becomes unimportant, and thermal gradients can only be developed and maintained in the crust and the envelope. In the crust, the dissipative processes responsible for the finite thermal conductivity include all the mutual interactions between electrons, lattice phonons (collective motion of ions in the solid phase), impurities (defects in the lattice), superfluid phonons (collective motion of superfluid neutrons) or normal neutrons. The mean free path of free neutrons, which is limited by the interactions with the lattice, is expected to be much shorter than for the electrons, but a fully consistent calculation is yet to be done (Chamel 2008). Quantizing effects due to the presence of a strong magnetic field become important only in the envelope, or in the outer crust for very large magnetic fields (B 10 15 G). For comparison, we also plot the B = 0 values. The quantizing effects are visible as oscillations around the classical (non-magnetic) values, corresponding to the gradual filling of Landau levels. More details about the calculation of the microphysics input (κ, c v , Q) can be found in Sect. 2 of Potekhin et al (2015b).
We can understand how and where anisotropy becomes relevant by considering electron conductivity in the presence of a strong magnetic field (and for now, ignoring quantizing effects). The heat flux is then reduced to the compact form (Pérez-Azorín et al 2006): where b ≡ B/B is the unit vector in the local direction of the magnetic field. The heat flux is thus explicitly decomposed in three parts: heat flowing in the direction of the redshifted temperature gradient, ∇(e ν T ), heat flowing along magnetic field lines (direction of b), and heat flowing in the direction perpendicular to both.
In the low-density region (envelope and atmosphere), radiative equilibrium will be established much faster than the interior evolves. The difference by many orders of magnitude of the thermal relaxation timescales between the envelope and the interior (crust and core) makes computationally unpractical to perform cooling simulations in a numerical grid including all layers up to the star surface. Therefore, the outer layer is effectively treated as a boundary condition. It relies on a separate calculation of stationary envelope models to obtain a functional fit giving a relation between the surface temperature T s , which determines the radiation flux, and the temperature T b at the crust/envelope boundary. This T s − T b relation provides the outer boundary condition to the heat transfer equation. The radiation from the surface is usually assumed to be blackbody radiation, although the alternative possibility of more elaborated atmosphere models, or anisotropic radiation from a condensed surface, have also been studied (Turolla et al 2004;van Adelsberg et al 2005;Pérez-Azorín et al 2005;Potekhin et al 2012). A historical review and modern examples of such envelope models are discussed in Sect. 5 of Potekhin et al (2015b). Models include different values for the curst/envelope boundary density, magnetic field intensity and geometry, and chemical composition (which is uncertain).
The first 2D models of the stationary thermal structure in a realistic context (including the comparison to observational data) were obtained by Geppert et al (2004Geppert et al ( , 2006 and Pérez-Azorín et al (2006), paving the road for subsequent 2D simulations of the time evolution of temperature in strongly magnetized NS (Aguilera et al 2008b,a;Kaminker et al 2014). In all these works, the magnetic field was held fixed, as a background, exploring different possibilities, including superstrong (B ∼ 10 15 -10 16 G) toroidal magnetic fields in the crust to explain the strongly non-uniform distribution of the surface temperature. Only recently , the fully coupled evolution of temperature and magnetic field has been studied with detailed numerical simulations. In the remaining of this section, we focus on the main aspects of the numerical methods employed to solve Eq. (4) alone, and we will return to the specific problems originated by the coupling with the magnetic evolution in the following sections.

Numerical methods for 2D cooling
There are two general strategies to solve the heat equation: spectral methods and finite-difference schemes. Spectral methods are well known to be elegant, accurate and efficient for solving partial differential equations with parabolic and elliptic terms, where Laplacian (or similar) operators are present. However, they are much more tedious to implement and to be modified, and usually require some strong previous mathematical understanding. On the contrary, finite-difference schemes are very easy to implement and do not require any complex theoretical background before they can be applied. On the negative side, finite-difference schemes are less efficient and accurate, when compared to spectral methods using the same amount of computational resources. The choice of one over the other is mostly a matter of taste. However, in realistic problems with "dirty" microphysics (irregular or discontinuous coefficients, stiff source-terms, quantities varying many orders of magnitude, etc), simpler finite-difference schemes are usually more robust and more flexible than the heavy mathematical machinery normally carried along with spectral methods, which are often derived for constant microphysical parameters. For this last reason, here we will discuss the use of finite-difference methods to solve our particular problem.
Let us consider the energy balance equation (4), with the flux given by Eq. (7). We first note that, in axial symmetry, the ϕ−component of the flux is generally non-zero but need not to be evaluated since it is independent of ϕ, so that its contribution to the flux divergence vanishes. For example, in the case of a purely poloidal field (only r, θ components), we can ignore the last term in Eq. (7) because it does not result in the time variation of the temperature. However, in the presence of a significant toroidal component B ϕ , the last term gives a non-negligible contribution to the heat flux in the direction perpendicular to ∇(e ν T ) (it acts as a Hall-like term).
In Aguilera et al (2008b,a); Viganò et al (2013) and related works, they assume axial symmetry and adopt a finite-differences numerical scheme. Values of temperature are defined at the center of each cell, where also the heating rate and the neutrino losses are evaluated, while fluxes are calculated at each cell-edge, as illustrated in Fig. 4. The boundary conditions at the center (r = 0) are simply F = 0, while on the axis the non-radial components of the flux must vanish. As an outer boundary, they consider the crust/envelope interface, r = R b , where the outgoing radial flux, F out , is given by a formula depending on the values of T b and B in the last numerical cell. For example, assuming blackbody emission from the surface, for each outermost numerical cell, characterized by an outer surface Σ r and a given value of T b and B, one has F out = σ B Σ r T 4 s where σ B is the Stefan-Boltzmann constant, and T s is given by the T s − T b relation (dependent on B), as discussed in the previous subsection.
To overcome the strong limitation on the time step in the heat equation, ∆t ∝ (∆ x) 2 , the diffusion equation can be discretized in time in a semi-implicit or fully implicit way, which results in a linear system of equations described by a block tridiagonal matrix (Richtmyer and Morton 1967). The "unknowns" vector, formed by the temperatures in each cell, is advanced by inverting the matrix with standard numerical techniques for linear algebra problems, like the lower-upper (LU) decomposition, a common Gauss elimination based method for general matrices, available in open source packages like LAPACK. However, this is not the most efficient method for large matrices. A particular adaptation of the Gauss elimination to the blocktridiagonal systems, known as Thomas algorithm Thomas (1949) or matrix-sweeping algorithm, is much more efficient, but its parallelization is limited to the operations within each of the block matrices. A new idea that has been proposed to overcome parallelization restrictions is to combine the Thomas method with a different decomposition of the block tridiagonal matrix (Belov et al 2017). A word of caution is in order regarding the treatment of the source term. The thermal evolution during the first Myr is strongly dominated by neutrino emission processes, which enter the evolution equation through a very stiff source term, typically a power-law of the temperature with a high index (T 8 for modified URCA processes, T 6 for direct URCA processes). These source terms cannot be handled explicitly without reducing the time step to unacceptable small values but, since they are local rates, linearization followed by a fully implicit discretization is straightforward and results in the redefinition of the source vector and the diagonal terms of the matrix. A very basic description to deal with stiff source terms can be found in Sect. 17.5 of Press et al (2007). This procedure is stable, at the cost of losing some precision, but it can be improved by using more elaborated implicit-explicit Runge-Kutta algorithms (Koto 2008).

Temperature anisotropy in a magnetized neutron star
An analytical solution that can be used to test numerical codes in multi-dimensions is the evolution of a thermal pulse in an infinite medium, embedded in a homogeneous magnetic field oriented along the z-axis, which causes the anisotropic diffusion of heat. Assuming constant conductivities, and neglecting relativistic effects, the following analytical solution for the temperature profile can be obtained for t > t 0 : where T 0 is the central temperature at the initial time t 0 . In Fig. 5 we show the comparison between the analytical (solid) and numerical (stars) solution for a model with t 0 = 10 −4 , T 0 = 1, κ ⊥ = 10 2 and ω B τ e = 3. The boundary conditions employed are F = 0 at the center and the temperature corresponding to the analytical solution at the surface (r = 1). Pérez-Azorín et al (2006) found deviations from the analytical solution to be less than 0.1% in any particular cell within the entire domain, even with a relatively low grid resolution of 100 radial zones and 40 angular zones.
To conclude this section, the induced anisotropy in a realistic NS reported by Pérez-Azorín et al (2006) is shown in Fig. 6. The figure shows equilibrium thermal solutions, in the absence of heat sources and sinks. The core temperature is kept at 5 × 10 7 K, and the surface boundary condition is given by the T s − T b relation, assuming blackbody radiation. The poloidal component is the same in all models (B p = 10 13 G). The effect of the magnetic field on the temperature distribution can be easily understood by examining the expression of the heat flux (7). When ω B τ e 1, the dominant contribution to the flux is parallel to the magnetic field and proportional to b · ∇(e ν T ). Thus, in the stationary regime (i.e., ∇ · (e 2ν F) = 0 if no sources are present), the temperature distribution must be such that b ⊥ ∇(e ν T ): magnetic field lines are tangent to surfaces of constant temperature. This is explicitly visible in the left panel, which corresponds to the stationary solution for a purely poloidal configuration with a core temperature of 5 × 10 7 K. Only near the surface, the large temperature gradient can result in a significant heat flux across the magnetic field lines. When we add a strong toroidal component, the Hall term (proportional to ω B τ e ) in Eq. (7), activates meridional heat fluxes which lead to a nearly isothermal crust. The central panel shows the temperature distribution for a force-free magnetic field with a global toroidal component, present in both the crust and the envelope. The right panel shows a third model with a strong toroidal component confined to a thin crustal region (dashed lines). It acts as an insulator maintaining a temperature gradient between both sides of the toroidal field.
3 Magnetic field evolution in the interior of neutron stars: theory review The interior of a NS is a complex multifluid system, where different species coexist and may have different average hydrodynamical velocities. In most of the crust, for instance, nuclei have very restricted mobility and form a solid lattice. Only the "electron fluid" can flow, providing the currents that sustain the magnetic field. In the inner crust superfluid neutrons are partially decoupled from the heavy nuclei, providing a third neutral component. In the core, the coexistence of superfluid neutrons and superconducting protons makes the situation even less clear. Since a full multifluid, reactive MHD-like description of the system is far from being affordable, one must rely on different levels of approximation that gradually incorporate the rele-vant physics. In this section we give an overview of the theory, trying to capture the most relevant processes governing the magnetic field evolution in a relatively simple mathematical form. For consistency with the previous section, we assume the same spherically symmetric background metric and we keep track of the most important relativistic corrections.
The evolution of the magnetic field is given by Faraday's induction law: which needs to be closed by the prescription of the electric field E in terms of the other variables (constituent component velocities and the magnetic field itself), either using simplifying assumptions (e.g., Ohm's law) or solving additional equations. Very often, this prescription involves the electrical current density, which in many MHD variations can be obtained from Ampére's law, neglecting the displacement currents In a complete multi-fluid description of plasmas, the set of hydrodynamic equations complements Faraday's law. From the multi-fluid hydrodynamics equations, a generalized Ohm's law -in which the electrical conductivity is a tensor -can be derived (Yakovlev and Shalybkov 1990;Shalybkov and Urpin 1995) j =σ E.
Expressing the tensor components in a basis referred to the magnetic field orientation, one can identify longitudinal, perpendicular and Hall components, that give rise to a complex structure when the equation is inverted to express E as a function of j, and B. However, in some regimes, one can make simplifications to make the problem affordable (Urpin and Yakovlev 1980;Jones 1988;Goldreich and Reisenegger 1992). The three main processes are Ohmic dissipation, Hall drift (only relevant in the crust) and ambipolar diffusion (only relevant in the core) (Goldreich and Reisenegger 1992;Shalybkov and Urpin 1995;Cumming et al 2004), although additional terms could in principle be also included in the induction equation. For instance, there are theoretical arguments proposing additional slow-motion dynamical terms, such as plastic flow (Beloborodov and Levin 2014;Lander 2016;Lander and Gourgouliatos 2019), magnetically induced superfluid flows (Ofengeim and Gusakov 2018) or vortex buoyancy (Muslimov and Tsygan 1985;Konenkov and Geppert 2000;Elfritz et al 2016;Dommes and Gusakov 2017). Typically, all these effects are introduced as advective terms, of the type E = −v × B, with v being some effective velocity. Thermoelectric effects have also been proposed to become significant in regions with large temperature gradients Wiebicke and Geppert 1991, 1995Geppert and Wiebicke 1995;Wiebicke and Geppert 1996); These additional terms are not included in most of the existing literature, and no detailed numerical simulations are known so far. However, some of them may play a more important role than expected and should be carefully revisited. Here, we review the principal characteristics of the most standard and better understood physical processes.

Ohmic dissipation
In the simplest case, the electric field in the reference frame comoving with matter is simply related to the electrical current density, j, by: where the conductivity σ , dominated by electrons, must take into account all the (usually temperature-dependent) collision processes of the charge carriers. Here, σ actually represents the longitudinal (to the magnetic field) component of the general conductivity tensorσ . In the weak field limit, the tensor becomes a scalar (σ ≡ σ ) times the identity, and possible anisotropic effects are absent. The induction equation, when we have only Ohmic dissipation, conforms a vector diffusion equation: where we have defined the magnetic diffusivity η ≡ c 2 4πσ . In the relaxation time approximation, the electrical conductivity parallel to the magnetic field, σ = e 2 n e τ e /m * e , with n e being the electron number density. Typical values of the electrical conductivity in the crust are σ ∼ 10 22 -10 25 s −1 , several orders of magnitude larger than in the most conductive terrestrial metals described by the band theory in solid state physics. In the core, the even larger electrical conductivity (σ ∼ 10 26 -10 29 s −1 ) results in much longer Ohmic timescales, thus potentially affecting the magnetic field evolution only at a very late stage (t 10 8 yr), when isolated NSs are too cold to be observed. In Fig. 7 we show typical profiles of the electrical conductivity, for the same combinations of T and B shown for the thermal conductivity in Fig. 3. Since, neglecting inelastic scattering, both thermal and electrical conductivities are proportional to the collision time τ e , they share some trends: the suppression of the conduction in the direction orthogonal to a strong magnetic field, and the quantizing effects visible as oscillations around the classical value (Potekhin et al 2015b;Potekhin and Chabrier 2018). We note that, if inelastic scattering contributes significantly, τ e can be different for thermal and electrical conductivities.

The Hall drift
At the next level of approximation, one must consider not only Ohmic dissipation but also advection of the magnetic field lines by the charged component of the fluid, say the electrons, with velocity v e . The electric field has the following form In the crust, the electron velocity is simply proportional to the electric current Here, the first term on the right-hand side is the same as in Eq. (12) and accounts for Ohmic dissipation, while the second term is the nonlinear Hall term. Note that the latter does not depend on the temperature, but it varies by orders of magnitude in the crust due to the inverse dependence with density. We can factor out the magnetic diffusivity and express the Hall induction equation in the form This form of the induction equation makes explicit that the magnetization parameter ω B τ e , which also determined the degree of anisotropy in the heat transfer, Eq. (6), plays the role of the magnetic Reynolds number: it gives the relative weight of the Hall and Ohmic dissipation terms. Generally speaking, as we approach the surface from the interior, ω B τ e increases. We note that, given these considerations, one has to be careful interpreting analytical estimates of the Ohmic or Hall timescales, since both vary by many orders of magnitude depending on the local conditions. The vast majority of the existing studies of magnetic field evolution in NS crust (Hollerbach andRüdiger 2002, 2004;Pons and Geppert 2007;Reisenegger et al 2007;Pons et al 2009a;Kondić et al 2011;Viganò et al 2012Viganò et al , 2013Gourgouliatos et al 2013;Marchant et al 2014;Cumming 2014b, 2015;Wood and Hollerbach 2015) are restricted to 2D simulations, but the few recent 3D simulations suggest that the main aspects of 2D results partially hold: although the Hall term itself conserves energy, the creation of small-scale structures results in an enhanced Ohmic dissipation. Some distinctive 3D features are the Hall-induced, small scale, azimuthal magnetic structures that seem to persist on long timescales (see Sect. 7).

Plasticity and crustal failures
The main idea for the Hall-MHD description of the crust is that ions are locked in the crustal lattice and only electrons are mobile. However, molecular dynamics simulations (Horowitz and Kadau 2009) show that the matter has an elastic behavior until certain maximum stress. Above it, the magnetic stresses, quantified by the Maxwell tensor M ≡ B i B j /4π, cannot be compensated by the elastic response (a more rigorous global condition is the von Mises criterion applied in Lander et al 2015). Crustal failures are treated in the most simplified manner as star-quakes. By evaluating the accumulated stress, Pons and Perna (2011); Perna and Pons (2011) simulated the frequency and energetics of the internal magnetic rearrangements, which was proposed to be at the origin of magnetar outbursts. This model mimics earthquakes since, under terrestrial conditions, the low densities of the material allow for propagation of sudden fractures: the Earth mantle in this sense can be thought as brittle. However, materials subject to very slow shearing forces could behave differently and enter a plastic regime where, instead of sudden crustal failures, a slow plastic flow takes place. Despite the different dynamics, the energetic arguments relating the release of energy due to the accumulation of magnetic stresses are similar. Recent simulations (Lander and Gourgouliatos 2019) show the features of such plastic flow under the assumption of Stokes flow, where a viscous term balances magnetic and elastic stresses. They compare the crustal response under Ohmic and Hall evolution and find that there can be significant plastic-like motions in the external layers of the star. Similar arguments have also been proposed to account for the deposition of heat by the visco-plastic flow and the propagation of thermo-plastic waves (Beloborodov and Levin 2014). Depending on which hypotheses we make, the interpretation of the velocities in the advective term (v × B) of the induction equation requires a proper physical and mathematical approach.

Ambipolar diffusion in neutron star cores
The number of works concerning mechanisms operating in NS cores is sensibly smaller, and most contain far less detail than the studies of the crust. Owing to its cubic dependence on B, ambipolar diffusion could be the dominant process driving the evolution of magnetars during the first 10 3 − 10 5 yr, although there is some controversy. In particular, we refer the reader interested in the role of chemical potential gradients, which is out of the scope of this review, to the literature. For example, Goldreich and Reisenegger (1992) or Passamonti et al (2017b) derived an elliptic equation from the continuity and momentum equations to determine the small deviations from beta equilibrium. However, Gusakov et al (2017) question the validity of that approach in stratified matter, and obtain a different equation from the momentum equation (implicitly assuming magnetostatic equilibrium), in which the small deviations of the chemical potentials from their equilibrium values do not depend on temperature and are determined by the Lorentz force. With the same methodology, Ofengeim and Gusakov (2018) calculate the instantaneous particle velocities and other parameters of interest, determined by specifying the magnetic field configuration, and found that the evolution timescales could be shorter than expected.
The short way to incorporate ambipolar diffusion is to generalize the form of the electric field by introducing the "ambipolar velocity" v a : The simplest case is realized in the regime where the system attains β −equilibrium faster than it evolves, and the ambipolar velocity is proportional to the Lorentz force where f a is a positive-defined drag coefficient. For simplicity we only consider this case in the next sections. We also note that, alternatively, the ambipolar term can be written as: where it explicitly takes the form of a resistive-like term, with a B 2 -dependent coefficient, only acting on the currents perpendicular to the magnetic field (j ⊥ ) aligning the magnetic field with the current and bringing the system into a force-free configuration, characterized by definition by j × B = 0. It is important to remark that the effect of this term is very sensible to the magnetic geometry, besides its strength: it has no consequences on the current flowing along magnetic field lines. This property has been used to introduce a formally similar term (differing only by a re-normalization factor ∝ 1/B 2 ) in the so-called magneto-frictional method, used to obtain configurations of twisted force-free solar (Roumeliotis et al 1994) and NS (Viganò et al 2011) magnetospheres (see also Sect. 6.3).
Most previous works studying ambipolar diffusion rely on timescale estimates, with few exceptions. Simulations are only available in a simplified 1D approach (Hoyos et al 2008(Hoyos et al , 2010 and very recently in 2D (Castillo et al 2017;Passamonti et al 2017b;Bransgrove et al 2018), usually for constant coefficients. However, in a realistic scenario, there is a further complication. The NS core cools down below the neutron-superfluid and proton-superconducting critical temperatures very fast, which has important implications, sometimes controversial. Goldreich and Reisenegger (1992) argued that ambipolar diffusion would still be a significant process, but Glampedakis et al (2011) studied in detail the ambipolar diffusion in superfluid and superconducting stars and concluded that its role on the magnetic field evolution would be negligible. Other recent works (Graber et al 2015;Elfritz et al 2016) have also shown that, without considering ambipolar diffusion, the magnetic flux expulsion from the NS core with superconducting protons is very slow. In Passamonti et al (2017a) the various approximations employed to study the long-term evolution of the magnetic field in NS cores were revisited, solving a recent controversy (Graber et al 2015;Dommes and Gusakov 2017) on the correct form of the induction equation and the relevant evolution timescale in superconducting NS cores.

Mathematical structure of the generalized induction equation
In order to understand the dynamical evolution of the system and to design a successful numerical algorithm, it is important to identify the mathematical character of the equations and the wave modes. The magnitude of ω B τ e defines the transition from a purely parabolic equation (ω B τ e 1) to a hyperbolic regime (ω B τ e 1). The Hall term introduces two wave modes into the system. Huba (2003) has shown that, in a constant density medium, the only modes of the Hall-MHD equation are the whistler or helicon waves. They are transverse field perturbations propagating along the field lines. In presence of a charge density gradient, additional Hall drift waves appear. These are transverse modes that propagate in the B × ∇n e direction. We also note that the presence of charge density gradients results in a Burgers-like term (Vainshtein et al 2000). Furthermore, even in the constant density case but without planar symmetry, the evolution of the toroidal component also contains a quadratic term that resembles the Burgers equation (Pons and Geppert 2007) with a coefficient dependent on the distance to the axis. This term leads to the formation of discontinuous solutions (current sheets) that require proper treatment. It is fundamental for a numerical Hall-MHD code to reproduce these modes and features, which are easily testable, as illustrated in Sect. 5.
In Viganò et al (2019) they give a complete description of the characteristic structure of the induction equation, including the Ohmic, Hall and ambipolar terms, in a flat spacetime, e ν = e λ = 1. By assuming a generic perturbation over a fixed background field B o : with a wavelength k much shorter than any other typical length of the system (typical variation scales of the Ohmic, ambipolar and Hall pre-coefficients), the eigenvalues are given by where B ok =k · B o , and B op = |B o − B okk |. This relation explicitly confirms that the Hall term is the only one that could be associated with waves (take the limit η = f a = 0), while the Ohmic and ambipolar terms are intrinsically dissipative.

Magnetic field evolution in the interior of neutron stars: numerical methods
In this section, we go through the most relevant aspects of numerical methods. The first important choice is the formalism to be adopted. There are two options: i) to work directly with the magnetic field components, which does not require any further mathematical manipulation but implies to care about how to preserve the divergence-free condition, and ii) exploiting the solenoidal constraint to work with only two functions representing the two true degrees of freedom instead of three components: the socalled poloidal-toroidal decomposition (see Appendix A). Finite-difference schemes have been developed for both formalisms, while spectral methods more often built on the poloidal-toroidal decomposition. We begin with an overview of spectral methods, before turning into some key aspects of finite-difference schemes.

Spectral methods with the toroidal-poloidal decomposition
Using the notation of Geppert and Wiebicke (1991), the basic idea is to expand the poloidal (Φ) and toroidal (Ψ ) scalar functions in a series of spherical harmonics where n = 1, . . . , n max and m = −n, . . . , +n.
Assuming a radial dependent diffusivity, η = η(r), it can be shown that the Ohmic term for each multipole effectively decouples, and the set of coupled evolution equations for the radial parts (Φ nm and Ψ nm ) can be readily obtained : where we use D nm and C nm as a shorthand for the nonlinear Hall terms (the full expressions can also be found in Geppert and Wiebicke 1991). These include sums over running indices and coupling constants related to Clebsch-Gordan coefficients (the sum rules to combine angular momentum operators are used to determine which multipoles are coupled to each other). All these coefficients can be evaluated once at the beginning of the evolution and stored in a memory-saving form since only specific combinations of indices are non-zero.
In the most general case, however, the magnetic diffusivity also depends on the angular coordinates, for example through the temperature dependence of η when the temperature is non-uniform. In this case we can also expand the magnetic diffusivity in spherical harmonics where the sum must include the monopole term, n = 0, . . . , n max . These new terms couple different multipoles of the same component (poloidal or toroidal). The inclusion of additional terms in the electric field (e.g. ambipolar diffusion) would introduce even more complicated non-linear couplings (the theory has not yet been developed).
In general, we end up with a system of the order of ≈ 2n 2 max , strongly coupled, differential equations. The choice now is whether using a different spectral decomposition in the radial direction (usually Chebyshev polynomials) or employing a hybrid method, applying standard finite-difference techniques in the radial direction to solve the system of equations.
The first multi-dimensional (2D) simulations of the evolution of the crustal magnetic field assumed a constant density shell (Hollerbach and Rüdiger 2002) and were later extended to include density gradients (Hollerbach and Rüdiger 2004). They used an adapted version of the spherical harmonic code described in Hollerbach (2000), including modes up to l = 100, and 25 Chebychev polynomials in the radial direction, but they were restricted to ω B τ e < 200 by numerical issues. In Pons and Geppert (2007); Pons et al (2009b), they used a hybrid code (spectral in angles but finitedifferences in the radial direction) to perform 2D simulations in realistic profiles of NSs over relevant timescales (typically, Myr). This approach allowed us to reach higher values of the magnetization parameter (ω B τ e ≈ 10 3 ), and to study the Hall instability (Pons and Geppert 2010). The same approach is used in the 3D simulations of Wood and Hollerbach (2015); Gourgouliatos et al (2016), which were limited to magnetization parameters of the order of 100. The main problem arises from the presence of non-linear Burgers-like terms, which naturally lead to discontinuities (see § 3.5), which are notoriously poorly handled by spectral codes. For this reason, subsequent works aiming at extending the simulations to more general cases have been gradually shifting towards the use of finite-difference schemes.

Finite-difference and finite-volume schemes
To study the interesting magnetar scenario in detail, the numerical codes must be able to go a bit further. In Viganò et al (2012), a novel approach making use of the wellknow High-Resolution Shock-Capturing (HRSC) techniques (Toro 1997), designed to handle shocks in hydrodynamics and MHD, was proposed. These techniques have been successfully applied to a range of problems, from a simple 1D Burgers equation to complex ideal MHD problems (Antón et al 2006;Giacomazzo and Rezzolla 2007;Cerdá-Durán et al 2008), avoiding the appearance of spurious oscillations near discontinuities. We refer to Martí and Müller (2015) for a general review on grid-based methods and to Balsara (2017) for a review on finite-volume methods, applied to other astrophysical scenarios. Let us review some of the main characteristics of these methods, of particular interest in our problem.

Conservation form and staggered grids
In hydrodynamics and MHD, the system of partial differential equations (PDEs) involve the divergence operator acting on vector or tensor fields. Thus, Gauss' theorem is usually employed in the design of the algorithms, exploiting the formulation of the equations in conservation form. Analogously, for problems involving the induction equation, the presence of the curl operator makes it natural to apply Stokes' theorem to the equation. Considering a numerical cell and its surface Σ α normal to the α direction, delimited by the curve C Σ , we have a discretized version of eq. (9): The space-discretized evolution equation for the average of the magnetic field component normal to the surface over the cell surface is then Here, the circulation of the electric field is approximated by the sum ∑ k E k l k , where E k is the average value of the electric field over each cell of length l k , and k identifies each of the four edges of the face. For clarity, in this section, we omit relativistic metric factors that must be consistently incorporated in the definitions of lengths, areas, and volumes. The problem is then reduced to design an accurate and stable discretization method to calculate the E k components at each edge. A natural choice is to use staggered grids, for which in each numerical cell the locations of the different field components are conveniently displaced, instead of being all located at the same position (typically, the center), as in standard centered schemes. In our case, we allocate the normal magnetic field components at each face center and electric field components along cell edges. Fig. 8 shows an example of the location of the variables in a numerical cell in spherical coordinates (r, θ , ϕ), considering axial symmetry (in the general 3D case, there would be a displacement of B ϕ , E θ , E r in the direction orthogonal to the plane of the figure). Making use of Gauss' theorem, the numerical divergence can be evaluated, for each cell with volume ∆V , as follows: Fig. 9 Illustration of the procedure to calculate the electric field in a staggered grid: location of the components of velocity (red arrows) and magnetic field (blue) involved in the definition of contribution to E ϕ (black dot) from the Hall term.
With this definition, the divergence-preserving character of the methods using the conservation form and advancing in time B α components, becomes evident: taking the time derivative of eq. (27), and using eq. (26), every edge contributes twice with a different sign and cancels out. By construction, the divergence condition is preserved to machine error for any divergence-free initial data. Examples of applications of such methods can be found, among many others, in Tóth (2000); Viganò et al (2012); Balsara and Dumbser (2015).

Evaluation of the current and the electric field
Let us consider a general electric field of the form: where nonlinear (Hall and/or ambipolar) dependences on the magnetic field are implicitly contained in the expression of v.
By considering the allocation of the components in the staggered grid (Fig. 8), the components of the current density can be naturally defined along the edges of the cells, in the same positions as the electric field components, exploiting the discretized version of the Stokes' theorem applied to j ∝ ∇ × (e ν B). Therefore, the ohmic term in the electric field can be directly evaluated, but the other terms involving vector products require special care since they involve products of field components that are not defined at the same place as the desired electric field component. The simplest option is a direct interpolation of both v and B using the first neighbors, but this often results in numerical instabilities.
In the spirit of HRSC methods, we can instead think of the interpolated value of v as the advective velocity acting at that point (although it depends on B itself), and consistently take the upwind components B w α of the magnetic field at each interface. For example, in the axisymmetric case and considering the evolution of the poloidal components (B r , B θ ), the contributions of E r and E θ to the circulation cancel out and we only need to evaluate the contribution of E ϕ , which is given by In Fig. 9 we explicitly show the location of E ϕ (black point) and the location on the staggered grid of the quantities needed for its evaluation. First, v r and v θ are calculated taking the average of the two closest neighbors; in the example, they point outward and to the right, respectively. Second, one considers the upwind values of B w r and B w θ ; in the example, they are taken from the bottom and left sides.

Divergence cleaning methods in finite-difference schemes
An algorithm built on a staggered grid can be designed to preserve the divergence constraint by construction, but the different allocation of variables makes its implementation relatively complex, particularly in 3D problems and with the inclusion of quadratic and cubic terms in the electric field. Among alternative formulations that have recently gained popularity, and can also handle many MHD-like problems, a relatively simple option is the family of divergence-cleaning schemes built on standard grids (all components of the fields are defined and evolved at every grid node). A popular divergence-cleaning method (Dedner et al 2002), extensively used in MHD, consists in the extension of the system of equations as follows: where χ is a scalar field that allows the propagation and damping of divergence errors, and c h and γ are two parameters to be tuned: c h is the propagation speed of the constraint-violating modes, which decay exponentially on a timescale 1/γ. In principle, a large value of γ will damp and reduce divergence errors very quickly, but in practice the optimal cleaning is reached for c h ≈ γ ∼ O(1) because, if γ is too large, the source term becomes stiff and more difficult to handle with explicit numerical schemes.

Cell reconstruction and high-order accuracy
The original upwind (Godunov's) method is well known for its ability to capture discontinuous solutions, but it is only first-order accurate: the variables are assumed to be constant on each cell. This method can be easily extended to give second-order spatial accuracy on smooth solutions, but still avoiding non-physical oscillations near discontinuities, by using a reconstruction procedure that improves the piecewise constant approximation. A very popular choice for the slopes of the linear reconstructed function is the monotonized central-difference limiter, proposed by van Leer (1977). Given three consecutive points x i−1 , x i , x i+1 on a numerical grid, and the numerical values of the function f i−1 , f i , f i+1 , the reconstructed function within the cell i is given by The minmod function of three arguments is defined by Other popular higher order reconstructions, are PPM (Colella and Woodward 1984), PHM (Donat and Marquina 1996), MP5 (Suresh and Huynh 1997), the FDOC families (Bona et al 2009), or the Weighted-Essentially-Non-Oscillatory (WENO) reconstructions (Jiang and Shu 1996;Shu 1998;Yamaleev and Carpenter 2009;Balsara 2017). In Viganò et al (2019) they presented and thoroughly tested a two-step method consisting of the reconstruction with WENO methods of a combination of fluxes and fields at each node, known as flux-splitting (Shu 1998). This reconstruction scheme does not require the characteristic decomposition of the system of equations (i.e., the full spectrum of characteristic velocities) and, at the lowest order of reconstruction, their flux formula reduces to the popular and robust Local-Lax-Friedrichs flux (Toro 1997).

Courant condition and time advance
In explicit algorithms to solve PDEs involving propagating waves, the time step is limited by the Courant condition, which essentially states that waves cannot travel more than one cell length on each time step, avoiding numerical instabilities. Since we want to evolve our system on long (Ohmic) timescales, the Courant condition makes the simulation computationally expensive for Hall-dominated regimes, ω e B τ e 1. For each cell, we can estimate the Courant time related to the Hall term by where L is a typical distance in which the magnetic field varies (e.g., the curvature radius of the lines), ∆ l is the minimum length of the cell edges in any direction. In the case of a spectral code, ∆ l L dom /n max , i.e., the ratio between the length of the dominion and the maximum number of multipoles calculated. The Courant condition related to the ambipolar diffusion term is which becomes more restrictive than the Hall term when en c f a B 1. The Courant condition is then where k c is a factor < 1 and the minimum is calculated among all the numerical cells.
For test-bed problems in Cartesian coordinates, taking k c = 0.1 − 0.3 is usually sufficient. In realistic models, however, numerical instabilities caused by the quadratic dispersion relation of the whistler waves arise. It becomes particularly problematic with spherical coordinates unless we use a very restrictive k c ≈ 10 −3 . Recent work (González-Morales et al 2018) includes other stabilizing techniques introduced in O' Sullivan and Downes (2006) for the time advance of the non-linear terms. These techniques, namely the Super Time-Stepping and the Hall Diffusion Schemes, allow us to maintain stability and efficiently speed up the time evolution when the ambipolar or the Hall term dominates. Another common technique is the use of high-order dissipation (also called hyper-resistivity; Huba 2003), or a predictorcorrector step advancing alternatively different field components.
Viganò et al (2012) used a particularly simple method that significantly improves the stability of the scheme in spherical coordinates. Their procedure to advance the solution from t n to t n+1 = t n + ∆t can be summarized as follows: • starting from B n , all currents and electric field components are calculated t are used to calculate the modified current components and the toroidal part of the electric field E t : B n+1 t → J p → E t ; • finally, we use the values of E t to update the poloidal components E t → B n+1 p . In Tóth et al (2008), the authors discussed that such a two-stage formulation is equivalent to introduce a fourth-order hyper-resistivity. Since the toroidal component is advanced first, it follows that the hyper-resistive correction only acts on the evolution of the poloidal components. In Viganò et al (2012) it was also shown that the additional correction given by E t contains higher-order spatial derivatives and scales with (∆t) 2 , which is characteristic of hyper-resistive terms. They found a significant improvement in the stability of the method when comparing a fully explicit algorithm with the two-steps method, allowing to work with k c ≈ 10 −2 − 10 −1 .
In the finite-difference schemes of Viganò et al (2019), the authors used a fourthorder Runge-Kutta scheme and found that the instabilities are especially significant when using fifth-order-accurate methods for the flux reconstruction (i.e. WENO5), which needed to be combined with the application of artificial Kreiss-Oliger dissipation along each coordinate direction (Calabrese et al 2004). A sixth-order derivative dissipation operator has a similar stabilizing effect, filtering the high-frequency modes which can not be accurately resolved by the numerical grid, at the cost of a potential loss of accuracy  . For this reason, they recommend using third-order schemes, that do not require any additional artificial Kreiss-Oliger dissipation. The typical Courant factors used were again quite low, k c ≈ 10 −2 − 10 −1 .
The most advanced 3D code currently available (Wood and Hollerbach 2015;Gourgouliatos and Cumming 2015;Gourgouliatos et al , 2016Gourgouliatos and Hollerbach 2018) uses spherical harmonic expansions of the magnetic potential functions for the angular directions (see Appendix A), and a discretized grid in the radial one. The linear Ohmic terms are evaluated using a Crank-Nicolson scheme, while for the non-linear Hall terms an Adams-Bashforth scheme is used. The code is parallelized by considering spherical shells and uses the infrastructure of the PAR-ODY code (Dormy et al 1998;Aubert et al 2008). Further details are available in Gourgouliatos et al (2016).

Numerical tests
In order to calibrate the performance of numerical methods or algorithms, it is crucial to provide analytical solutions against which the numerical results can be confronted. Unfortunately, there are not many such solutions in the 3D case with arbitrary coefficients in the generalized Ohm's law. For reference, we collect in this section a number of testbed cases with analytical solutions (most of them used in previous works Viganò et al 2012Viganò et al , 2019, which probe different terms the induction equation. The successful completion of this battery of tests should be a good indicator of the performance of the codes. For the smooth tests below, § 5.1,5.2,5.4,5.5, one can also check the convergence order of the numerical scheme, by computing the dependence of the relative errors (assessed for instance by a L2-norm) on the resolution used. The remaining two tests, where discontinuities form, are instead useful to test the robustness of the code, because near discontinuous solutions the convergence reduces to first order, regardless of the scheme. In all the following tests, we work in the Newtonian limit, e ν = e λ = 1.

Whistler waves
We begin by considering the case when only the Hall term is present in the induction equation. In a constant density medium, the only modes of the Hall-MHD equation are the whistler or helicon waves (see § 3.5), which consist in transverse field perturbations propagating along the magnetic field lines (notably known also in the terrestrial ionosphere Helliwell 1965;Nunn 1974). The first test we discuss is to follow the correct propagation of whistler waves. Consider a two-dimensional slab, extending from z = −L to z = +L in the vertical direction, with periodic boundary conditions in the x-direction, and assume that all variables are independent of the y-coordinate. For the following initial magnetic field: where k x = nπ/L, n = 1, 2, ..., and B 1 B 0 , the linear regime admits a pure wave solution confined in the vertical direction and traveling in the x-direction, that is, the same eq. (34) replacing Here we have defined the reference Hall timescale as As an example, in Fig. 10 and attached movie -online version only -, we report the evolution of this initial configuration with B 0 = 10 3 B 1 , and k x L = π, from t = 0 to t = 2 (in units of τ 0 ), in a 200 × 50 Cartesian grid. The perturbations travel through the horizontal domain twice, with negligible dissipation or dispersion. Viganò et al (2012) ran the test for hundreds of Hall timescales without any indication of instabilities, even though electrical resistivity is set to zero. By varying the values k x and B 0 , one can confirm that the velocity of the perturbations in the simulation scales linearly with both parameters. An additional twist is to consider the same problem in a 2D or 3D box, but with an arbitrary rotation of the coordinates, in order to test the correct propagation in a more general direction, not aligned to any axis.

Hall drift waves
In the second test of the Hall term, we remove the assumption of a constant charge density background. In presence of a charge density gradient, additional transverse modes appear, the so-called Hall drift waves, which propagate in the B × ∇n e direction. Let us consider the same domain as in the previous test, but with a stratified background in the y-direction with n e given by n e (y) = n 0 1 + β L y , where n 0 is a reference density, with an associated Hall timescale τ 0 defined in eq. (36), and β L is a parameter with dimensions of inverse length. We apply periodic boundary conditions in the x-direction, while in the y-direction, an infinite domain can be simulated by copying the values of the magnetic field in the uppermost and lowermost cells (y = ±L) into their first neighbor ghost cells.
For the following initial configuration: and small perturbations (B 1 B 0 ), the solution at early times consists in pure Hall drift waves traveling in the x-direction with speed The solution in the linear regime can be obtained by replacing x by (x − v hd t) in eq. (38). For the particular model shown in Fig. 11, with B 0 = 10 3 B 1 , k x = π/2, L = 1 and β L = 0.2, we have a horizontal drift velocity of v hd = 0.2L/τ 0 , corresponding to a crossing time of 20 τ 0 . The figure shows the initial configuration of B z (top left) and the evolution of the perturbation after 0.5, 1.25 and 2 crossing times, respectively (top right, bottom left and bottom right). The shadow increases with the value of 1/n(y). For the Hall drift modes, the propagation velocity scales linearly with both B 0 and the gradient of n −1 e , but it is independent of the wavenumber of the perturbation. All these properties are correctly reproduced. After many cycles (the number depending on the B 1 /B 0 ratio), deviations from the purely advected, smooth solution begin to be visible. This is an expected non-linear effect that we discuss next.

The nonlinear regime and Burgers flows
With the two previous tests, we can check if a numerical code can reproduce the propagation of the fundamental modes at the correct speeds. However, these are valid solutions only in the linear regime. Let us consider more carefully the evolution of the B y component in a medium stratified in the z-direction. Assuming that B x = B z = 0, the governing equation reduces to: This is a version of the Burgers equation (which solution is well known) in the xdirection with a coefficient that depends on the z coordinate: If we consider the following initial configuration: Fig. 12 Horizontal section of the evolution of the initial configuration defined by eq. (43) with B 0 = 10 3 and k x L = π at t = 0 (crosses), t = 2τ 0 (triangles) and t = 4τ 0 (diamonds). The shock forms at t = 2τ 0 . The classical sawtooth shape developed during the evolution of the Burgers equation is evident. on a stratified background with n e (z) = n 0 1+β L z , we have g(z) = −β L L 2 /τ 0 B 0 and we can directly compare to the solution of the Burgers equation in one dimension, which evolves to form discontinuities from smooth initial data.
To handle this problem, we can make use of well-known HRSC numerical techniques to design a particular treatment to the quadratic term in B y 3 . A key issue is to consider the Burgers-like term in conservation form: whereF = g(z)B 2 y /2, which can now be treated with an upwind conservative method (Viganò et al 2012). In this case, the wave velocity determining the upwind direction is given by g(z)B y .
Expressing the evolution equations in conservative form is crucial when solving problems with shocks or other discontinuities, since non-conservative methods may result in the incorrect propagation speed of discontinuous solutions (Toro 2009). In Fig. 12 we show snapshots of the evolution of the initial conditions (43) with k x L = π, B 0 = 10 3 , and β L L = 0.2, taken from Viganò et al (2012). It follows the typical Burgers evolution. The wave breaking and the formation of a shock at t = 2τ 0 is clearly captured. We remark again that this test is done with zero physical resistivity, i.e., in the limit ω B τ → ∞, which is not reachable by spectral methods or centered- Fig. 13 Evolution of the purely Ohmic modes, eqs. (46)-(48), with α = 1, at t = 0, 1, 2, and 3 diffusion times (τ d ). The simulation has been run in a [−10, 10] 3 cubic domain, with a resolution of 128 3 equallyspaced Cartesian grid points with the Simflowny-based code . In the figure we compare the analytical (black lines) and numerical (color symbols) profiles of B z (x = 0, y = 0, z) (i.e., B r (r, θ = 0), top left), B z (x, y = 0, z = 0) (i.e., −B θ (r, θ = π/2), top right), and B x (x = 0, y, z = 0) (i.e., B ϕ (r, θ = π/2, ϕ = π/2), bottom). difference schemes in non-conservative form. In Viganò et al (2019), the reader can find more details about the solutions obtained with different reconstruction schemes.

Ohmic dissipation: self-similar axisymmetric force-free solutions
In spherical geometry, one of the few existing analytical solutions is the evolution of pure Ohmic dissipation modes. Considering the limit ω B τ → 0, and a constant η, the induction equation reads: It is straightforward to show that a force free magnetic field satisfying ∇ × B = αB, with constant α, is an Ohmic eigenmode, since the induction equation is reduced to Therefore, each component of the magnetic field decays exponentially with the diffusion timescale τ d = (ηα 2 ) −1 . We note that the evolution of each component is completely decoupled in this case. In spherical coordinates, the solutions of eq. (45) are described by factorized functions, which radial parts involve the spherical Bessel functions. The regularity condition at the center selects only one branch of the spherical Bessel functions (of the first kind), which, for the (l, m) = (1, 0) mode are where x = αr and k = ±1. With this initial condition, we follow the evolution of the modes during several τ d , until the magnetic field is almost completely dissipated.
As boundary conditions, we impose the analytical solutions for B θ and B ϕ . Fig. 13 compares the numerical (crosses) and analytical (solid lines) solutions of B r and B ϕ at different times, for a model with α = 1, and a [−10, 10] 3 cubic domain, run with the Simflowny-based code .

Ambipolar diffusion: the Barenblatt-Pattle solution
To test the ambipolar term, we now consider the case of a constant f a and set to zero the Hall and Ohmic coefficients. In axial symmetry and cylindrical coordinates, there exists an analytical solution corresponding to the diffusion of an infinitely long magnetic flux. Let us consider the evolution of the only component of the magnetic field B = B zẑ , in the direction of the fluxtube, which only depends on the cylindrical radial coordinate (ϖ). The currents are perpendicular to the magnetic field so that (j × B) × B = −B 2 j, and the induction equation with the ambipolar term is reduced to This form is analogous to the non-linear diffusion equation where m is a power index. The analytical 2D solutions proposed by Barenblatt and Pattle (Barenblatt 1952;Pattle 1959) consist of a delta function of integral Γ at the origin, which diffuses outwards with finite velocity. We note that the diffusion front is clearly defined, contrarily to the infinite front speed of a linear diffusion problem. This analytic solution can be explicitly written as follows: where d is the dimension of the problem, and α = (m − 1 + 2/d) −1 . The initial pulse spreads with a front located at a distance ϖ f from the origin, given by In Viganò et al (2019), they studied the evolution of the model with d = 2, m = 3, α = 1/3, and f a = 3, which gives the explicit solution In Fig. 14 we show three snapshots of the evolution, starting with t 0 = 1, and Γ = 1/18. The front propagates according to ϖ f (t) = t 1/6 . The numerical results correctly reproduce the expected shape of the expanding flux tube and the propagation speed of the front. The sharp discontinuity in the slope of B z near the front end was found to be well-reproduced even for low resolutions.

Evolution of a purely toroidal magnetic field
Finally, to conclude our proposed series of tests and examples, we consider the evolution of a pure toroidal magnetic field confined into a spherical shell, R core < r < R, under the combined action of both Ohmic dissipation and the Hall term. This case does not have an analytical solution, but we believe it is an important (yet relatively simple) test that can highlight some relevant issues. For simplicity, we impose as boundary conditions that all components of the magnetic field vanish at both boundaries. We consider the realistic NS background profile of Fig. 1, with R core = 10.8 km, R = 11.6 km, and we set a constant temperature of 10 8 K, which corresponds to a density-dependent magnetic diffusivity in the range η ∼ 0.01 − 10 km 2 /Myr. Our initial magnetic field is given by the following expression: where B 0 is a normalization factor adjusted to fix the initial maximum value that the toroidal magnetic field reaches across the star (denoted by B 0 t ). According to the Hall induction equation, any initial toroidal configuration must remain purely toroidal during the evolution, but its shape and location (and the associated currents) vary with time. As discussed in the literature (Hollerbach and Rüdiger 2002;Pons and Geppert 2007;Viganò et al 2012), the evolution has two characteristics: (i) a vertical drift, northward or southward depending on the sign of B ϕ ; (ii) a drift towards the interior of the star due to the existence of a charge density gradient. In the top panels of Figure 15 we show three snapshots of the evolution of an ultra-strong toroidal field (B 0 t = 3 × 10 15 G), such that the first effect (drift towards the equator of both rings) occurs faster. In this model the maximum value of ω B τ e is ≈ 300, although it varies throughout the crust. After 1000 yr, a radial current sheet (i.e., a sharp discontinuity of the toroidal magnetic field in the meridional direction) is created in the equator and Ohmic dissipation is locally enhanced. We also notice a global drift towards the interior: compare the distance to the surface of the models at t = 1000 and t = 3000 yr. The bottom panels show the evolution of the initial model with the reverse sign. We observe how the drift proceeds in the opposite direction, creating a strong toroidal ring around the axes, near the poles. This simple model is useful to understand the evolution when the initial toroidal field is dominant, even in the presence of a weaker poloidal field. As Geppert and Viganò (2014) have shown, starting with a very large fraction ( 99%) of magnetic energy stored in the toroidal component is a potential way to create local magnetic spots near the poles, where the lines are more concentrated and the magnetic field intensity can be one or two orders of magnitude higher than the average value.
In addition, because of the formation of the current sheet in the equator or the localized rings in the poles, this model is also useful to check several issues concerning energy conservation, numerical viscosity, and current sheet formation, amply discussed in Sect. 5.1 of Viganò et al (2012). A necessary test for any numerical code is to check the instantaneous (local and global) energy balance. To remark one of the Fig. 15 Evolution of a quadrupolar toroidal magnetic field, with a maximum initial value of 3 × 10 15 G confined into the NS crust. We show snapshots t = 0, 1000, and 3000 yr. Color contours show the toroidal magnetic field strength, the reddish corresponding to negative B ϕ and the yellowish to positive B ϕ . The only difference between the top and bottom panels is the reverse sign of the initial field. In the figure, the size of the crust has been amplified by a factor of 4 for better visualization. key points, let us recall the magnetic energy balance equation: where Q j = 4πη j 2 /c 2 is the Joule dissipation rate and S = cE×B/4π is the Poynting flux. During the evolution, the magnetic energy in a cell can only vary due to local Ohmic dissipation and by the interchange between neighbor cells (Poynting flux). Integrating eq. (55) over the volume of the numerical domain, we obtain the following energy balance equation: where E b = V (e ν B 2 /8π)dV is the total magnetic energy , Q tot = V e 2ν Q j dV the total Joule dissipation rate, and S tot = ∂V e 2ν S · dΣ the Poynting flux through the boundaries. Numerical instabilities usually show up as a strong violation of energy conservation and careful monitoring of the energy balance is a powerful diagnostic. This was one of the simplest possible initial configurations, yet capturing interesting physics. When we introduce an initial poloidal component, higher multipoles, stratified microphysical (Ohmic/ambipolar) coefficients, etc., it becomes non-trivial to design benchmark tests. In § 7, we give a summary of previous attempts to gradually approach realistic scenarios.

Magnetosphere-interior coupling and rotational evolution
An open issue in realistic simulations of the magnetic field evolution of NSs concerns the correct implementation of boundary conditions at the star surface. In the external region of a NS, the mass density is over twenty orders of magnitude smaller than in the outer crust, where numerical grids usually end. One needs to match two regions, with radically different physical conditions and timescales, through the thin (≈ 100 m) layer between them. The usual procedure is to assume that, on the slow secular evolution timescales, the exterior is immediately readjusted (light crossing timescale) to the stationary solution imposed by the surface values of magnetic fields and currents. In other words, for long timescales, the magnetosphere can be seen as a perfect conductor where currents quickly respond to cancel electromagnetic forces out. Thus, the interior evolution provides the surface values of the field that determine the external configuration. However, in a numerical code the interior also needs, at each time step, a recipe for the outer boundary condition to proceed with the evolution, so both problems are interlinked and must be consistently treated.
Under the assumption that the dynamics of the magnetosphere is dominated by the electro-magnetic field, and the plasma pressure as well as its inertia are negligible, a reasonable approximation is to consider that the large-scale structure of the magnetosphere is given by force-free configurations, in which the electric and magnetic forces on the plasma balance each other. For magnetar conditions, one can safely neglect the effects of rotation in the magnetospheric region near the star. Under this approximation, which we follow hereafter, the electric force is neglected (E = −v × B, with v c). Thus, the force-free condition reduces to j × B = 0: the electric currents flow parallel to the magnetic field lines that they sustain (since j ∝ ∇ × B, a force-free magnetic field is a Beltrami vector field). Within the family of possible solutions, the most trivial (and popular) one is the current-free, or potential solution, j = 0, which also holds in vacuum. Matching the interior magnetic field to a magnetospheric potential field is equivalent to physically avoid that the current escapes (enters) from (into) the star, although the non-vanishing Poynting flux across the boundary allows the two regions to interchange magnetic energy (but not magnetic helicity).
While the potential solution is acceptable as a first approximation, to advance toward more realistic models, we need more general solutions. As a matter of fact, electrical currents can stably flow in the closed magnetic field line region, similar to the Solar coronal loops. The current system lasts on relatively long timescales, from months to decades (Beloborodov 2009), presumably sustained by the interior dynamics. There is indirect observational evidence of such currents in some magnetars, where the presence of a plasma much denser than the Goldreich-Julian value has been inferred. Soft X-ray photons emitted from the star surface are up-scattered to higher energy (Lyutikov and Gavriil 2006;Rea et al 2008;Beloborodov 2013) through resonant Compton processes, resulting in the observed spectra. Equilibrium solutions of force-free twisted magnetospheres in the magnetar context were considered by several recent works (Fujisawa and Kisaka 2014;Glampedakis et al 2014;Pili et al 2015;Akgün et al 2016;Kojima 2017). However, the evolution of the interior sometimes leads to solutions that cannot be smoothly connected to a force-free solution. This implies discontinuities in the tangential components at the surface, corresponding to current sheets, that may cause numerical instabilities.
While rotation has negligible effects on the magnetic evolution, the opposite is not true: the spin period evolves due to electromagnetic torques determined by the magnetospheric configuration. Compared to the magnetic and thermal evolution, the equations describing the rotational evolution are simpler, but they predict the observable timing properties of isolated NSs. In the remaining of this section we review the methodology to prescribe boundary conditions to the magnetic field when one solves the induction equation with different types of code, commenting on some problems that arise at the practical level, and we provide the recipe for the rotational evolution.

Potential boundary conditions
Spectral methods. Using the same notation as in Sect. 4.1 for the poloidal/toroidal decomposition, the requirement that all components of the magnetic field be continuous (no current sheets at the surface) implies that the scalar potentials Φ nm and Ψ nm , and their derivatives ∂ Φ nm ∂ r , are continuous through the outer boundary. Therefore, the ∇ × B = 0 condition translates into and the following differential equation for each radial function Φ nm (r) where we assume the metric (1), and z ≡ 2GM c 2 r . We note that there is no m−dependence in the equation, so that the solution depends only on n and we will omit the m subindex hereafter.
In general, the family of solutions of Eq. (58) for any value of n can be expressed in terms of generalized hypergeometric functions (F([], [], z)), also known as Barnes' extended hypergeometric functions, as follows: where C n and D n are arbitrary integration constants that correspond to the weight of each magnetic multipole n. Note that regularity at r = ∞ requires D n = 0 for each n. For any given value of n, one can also express the solution in closed analytical form. The explicit expressions for n = 1 and n = 2 are If we consider the Newtonian limit (z → 0), Eq. (58) simplifies to: The only physical solution (regular at infinity) of this equation is Φ n = C n r −n . Therefore, the requirement of continuity across the surface results in In the relativistic case, we can implement Eq. (59) directly, or the most practical form, analogous to the Newtonian case: where the f n 's are relativistic corrections that only depend on the value of z at the star surface, z(r = R) (in the Newtonian limit all f n = 1), and can be evaluated numerically only once with the help of any algebraic manipulator and stored 4 .
Finite-difference schemes. If we do not use a spectral method, we must apply boundary conditions to magnetic field components instead of the individual multipoles. In general, we need to provide the field components, in one or more ghost cells outside the physical grid, in terms of the components at the last grid point. However, we still can make use of the previous form of the boundary conditions, as a relation between the poloidal radial function and its derivative for each multipole,. Let us explain an accurate and elegant procedure to impose the current-free constraint in the axisymmetric case. From Eq. (64) and the expression of the field components in terms of the poloidal and toroidal functions, one can easily show that the potential solution expansion in terms of Legendre polynomials (P l ) reads: where we denote by b l the weights of the multipoles At a practical level, one can proceed as follows: • First, at each time step, obtain the b l coefficients from the Legendre decomposition of the radial component of the magnetic field over the surface B r (r = R, θ ).
In a discretised scheme, values of b l can be calculated up to a maximum multipole l max = n θ /2, where n θ is the number of angular points of the grid. • Second, from the b l 's, reconstruct the values of B r and B θ in the external ghost cells, as required by the method, by using Eq. (66). • Finally, simply set B ϕ = 0 for any cell r ≥ R.
This method is very accurate for smooth functions B r . In the case of sharp features in B r , which may be created by the Hall term, the largest multipoles acquire a nonnegligible weight, and, since l max is limited, fake oscillations in the reconstructed B θ may appear (Gibbs phenomenon). An alternative method to impose potential boundary conditions is based on the Green's representation formula, a formalism often used in electrostatic problems able to correctly handle the angular discontinuities in the normal components. Details about the derivation of the Green's integral relation between B r and B θ at the surface are given in Appendix B.
Note that, in 3D, applying the potential boundary conditions is a challenge for parallelization. The easiest solution is to parallelize the dominion by spherical shells, so that the integration over the star's surface, needed by either the spherical harmonic expansion or the Green's method, is in charge of a single processor (as in Wood and Hollerbach 2015 and following works). If, on the other hand, the parallelization is done by geometrically optimized patches (cubic in the simplest case, Viganò et al 2019), then the star's surface would be covered by different processors. In this case, the calculation at each point depends on calculations done by other processors, thus enlarging the needed stencil. This results in an excessive intercommunication load and prevents optimal scaling.

Force-free boundary conditions
The construction of relativistic, axisymmetric, force-free magnetospheres for (nonrotating) magnetars is a well studied problem (see, e.g., Kojima 2017 and references therein). In Akgün et al (2018) the authors explored a method to impose such boundary conditions by solving the Grad-Shafranov equation, at each time step, to match the internal evolution of the star. Let us review their approach. Considering axial symmetry, the magnetic field can be written as follows: where P and T are functions defining the poloidal and toroidal components, respectively (see more details in Appendix A).
The force-free condition (j × B = 0) implies that the electrical currents flow along magnetic surfaces, which are defined by constant P. Thus, the mathematical requirement of a vanishing azimuthal component of the local Lorentz force implies that the poloidal and toroidal functions must be functions of one another, say T = T (P), that is, the poloidal and toroidal functions P and T are constant on the same magnetic surfaces 5 .
From the definition of the current, one can arrive at the so-called Grad-Shafranov equation where T (P) = dT /dP. The current-free limit (potential solution) is simply recovered by taking the right hand side equal to zero. In principle, there is an infinite family of external force-free solutions for a given radial magnetic field at the surface, because of the freedom to choose the functional form of T (P). The main problem of this approach is how to continuously match the arbitrary field configuration, resulting from the evolution in the crust, while enforcing the force-free solution outside. In the crust, any line bundle marked by a given magnetic flux P has in general different values of T because, internally, the force-free condition does not hold. As discussed in Akgün et al (2018), there is an intrinsic inconsistency in the possibly multi-valued function T (P), if we strictly take it from the values at the surface (r = R). They address this problem by symmetrizing the numerical function T (P), which is physically equivalent to allow to propagate through the surface only the modes compatible with solutions of the Grad-Shafranov equation. This is motivated by the results from MHD simulations of the propagation of internal torsional oscillations (Gabler et al 2014), who found that antisymmetric modes cannot propagate into the magnetosphere and are reflected back into the interior.
In Fig. 16, we show the evolution of a magnetospheric configuration physically connected to the interior. The initial model consists of both poloidal and toroidal dipolar components, with the latter extending beyond the surface. As the internal magnetic field evolves, the external magnetic field is consistently twisted, by the injection of magnetic helicity (i.e., currents) in the magnetosphere. Force-free solutions are calculated at each time step until a critical point, where numerical solutions cannot be found anymore. At this point, the magnetosphere is expected to become unstable, resulting in a global reconfiguration with the opening of the twisted field lines and magnetic reconnection. This mechanism was studied in more detail in force-free electrodynamics simulations in 2D and 3D, and both the Newtonian and general relativistic cases (Parfrey et al 2013;Carrasco et al 2019).

Extended domains
An alternative to imposing a precise mathematical boundary condition at the surface is to consider an extended domain, where we evolve at the same time all compo- nents of the field, but with physical coefficients that enforce the solution to meet the required conditions. Instead of imposing a boundary condition at the last numerical cell, this approach considers a generalized induction equation, where, at the surface, there is a sharp transition in the values of the pre-coefficients describing the physics (η, 1/n e , f a ). In the numerical GRMHD context, this approach has been successfully used to describe at the same time the resistive and ideal MHD inside and outside a NS (Palenzuela 2013).
The idea is that, since the magnetospheric timescales are many orders of magnitude shorter than the interior, the long-term evolution of the magnetosphere can be seen as a series of equilibrium states, attained immediately after every time step of the interior. Therefore, one can activate an artificial term that dynamically leads to the force-free solution. This approach is similar to the magneto-frictional method (Yang et al 1986;Roumeliotis et al 1994), as known in solar physics. The modified induction equation employed in the exterior of the star has a mathematical structure equivalent to an ambipolar term, which forces currents to gradually align to magnetic field lines, without having to solve the elliptical Grad-Shafranov equation at every time step (which is numerically expensive). This also allows us to account for the transfer of helicity and provides a mechanism to continuously feed currents that twist the magnetosphere. The caveat is that the ambipolar coefficient must be fine-tuned to prevent the exterior dynamics from being neither too fast (it would excessively limit the time step), nor too slow (it would not manage to relax to a force-free configuration and would cause non-negligible, unphysical feedback on the interior).
In our context, such a strategy has only been explored preliminarily in the 3D Cartesian parallelized code used in Viganò et al (2019). This code was built by using Simflowny (Arbona et al 2013(Arbona et al , 2018, a versatile platform able to automatically generate parallelized codes for partial differential equations. It employs the adaptive mesh refinement libraries of SAMRAI (Hornung and Kohn 2002), and a graphical user interface that easily allows us to implement equations and to choose among different time and space discretization schemes. The code has not been applied yet to realistic simulations and is very different from previous codes. The Cartesian grid, with parallelization by regular cubic patches, imposes numerical challenges. One problem associated to Cartesian grids is that the geometry is not adapted to the physically preferred radial direction, along which gradients are usually much larger. Thus, one cannot improve the resolution in the radial direction alone, causing a rapid increase in the computational cost, compared to spherical coordinates-based codes (∝ N 3 instead of N r ). Furthermore, the Cartesian discretization implies the appearance of numerical noise at the (physically spherical) surface and crust/core interfaces (where the precoefficients in the induction equation show a sharp transition as mentioned above). Finally, one has also to take care that the outer domain (placed far enough from the surface by using different mesh refinements) does not introduce noise and allow a regular solution at infinity. These challenges need to be tackled soon to make this alternative method, still at its infancy, numerically feasible.

Evolution of spin period and obliquity
As NSs age, they spin down because of angular momentum losses due to magnetospheric torques (Spitkovsky 2006;Beskin et al 2013;Philippov et al 2014). This mechanism is effectively ruled by the dipolar component since other multipoles decay faster with the distance to the star and are negligible. Thus, the magnetic field evolution not only affects the surface temperature but also determines the rotational properties of the star. In the general case, the equations describing the coupled evolution of the spin period P and the angle between the magnetic dipolar moment and the rotation axis, χ, are (Philippov et al 2014): where we have defined the auxiliary quantity I is the moment of inertia of the part of the star co-rotating with the magnetosphere, and B p is the value of the dipolar component of the magnetic field at the magnetic pole. The latter is a function of time and can be provided by the simulations of the internal magnetic field. The coefficients κ 0 , κ 1 , κ 2 depend on the magnetosphere geometry and its physical conditions and determine the magnetospheric torque. For the classical vacuum dipole formula κ 0 = 0, κ 1 = κ 2 = 2/3, while for a realistic plasma-filled magnetosphere, κ 0 ≈ κ 1 ≈ 1, with the last coefficient varying between 0 and 1, depending on the assumptions. This coefficient can be fitted from results from 3D simulations for force-free and resistive magnetospheres (Philippov et al 2014), who found that the alignment of the rotation and magnetic axis in pulsars with vacuum magnetospheres proceeds much faster (exponential, with characteristic time τ 0 = P 0 2β B 0 ) than for realistic plasma-filled magnetospheres (a power-law).
A very important remark is that using the classical dipole formula, with κ 0 = 0, may mislead to the wrong conclusion that an aligned rotator does not exert any torque, thus stopping the spin-down of the star (Johnston and Karastergiou 2017). This is not physically correct, and realistic models predict κ 0 ≈ 1, which at most results in about a factor two correction (Spitkovsky 2006;Philippov et al 2014). Therefore, alignment cannot completely stop the period evolution, which can only happen if the magnetic field becomes negligible.
Besides the magnetic field strength and the inclination angle, I could also change with time. If there is a superfluid component (e.g., neutrons in the core or the inner crust), it is generally rotationally decoupled from the rest of the star. Therefore, it does not contribute to I, which only accounts for the matter rigidly co-rotating with the magnetosphere. Then, there are two possible effects. First, the volume of the superfluid component can change with time, since the phase transition depends on density and temperature (as the star cools down superfluid components occupy a larger volume, thus β slowly increases). Second, normal and superfluid components can suddenly and temporarily couple during glitches, modifyingf I. This can be formally considered with a two-fluid description or neglected, by assuming a constant I corresponding to the rigid co-rotation of the whole star. For realistic stars, the moment of inertia is I ∼ 1.5 × 10 45 g cm 2 , with a 50% uncertainty. This gives β ∼ 6 × 10 −40 s G −2 . We note that either angle variations (alignment) or moment of inertia variations, result in corrections to the torque by a factor 2, while magnetic field decay can result in torque variations of several orders of magnitude, and consequently in large and relatively fast variations of P andṖ.
Finally, we would like to address why the effects of rotation have negligible feedback on the magneto-thermal evolution. Magnetospheres of spinning isolated NSs have been studied in the last 50 years, analytically (starting from the seminal work of Goldreich and Julian 1969) and numerically for some geometries (mostly inclined rotating dipoles). Examples of 2D and 3D numerical simulations can be found in Spitkovsky (2006); Contopoulos et al (1999). First of all, the rotationallyinduced electric field, E = −(Ω r sin θ /c)φ × B, is negligible in the interior regions of the magnetosphere (r sin θ c/Ω ), thus justifying the non-rotating force-free approximation above, as a boundary condition at the star surface. Secondly, rotation opens up a bundle of lines close to the magnetic poles, which are stretched and twisted and are supposedly responsible for the radiated emission. The surface polar cap containing the footprints of the open lines has an opening angle of only θ ∼ arcsin[(RΩ /c) 1/2 ] ∼ 0.8 • P[s] −1/2 : a negligible fraction of the NS surface, especially for magnetars (P ∼ 1 − 10 s). A recent work (Karageorgopoulos et al 2019), based on the minimization of the Joule dissipation rate, has proposed how the rotationallyinduced polar currents close within the crust. The dissipated power by the Joule effect was estimated to be more than 10 orders of magnitude smaller than the rotational energy losses, which does not affect the cooling history.

Magneto-thermal evolution of neutron stars
The first question to answer when one plans to simulate the evolution of magnetic fields in NSs is the choice of the initial model. Since we are mostly interested in understanding the evolution of highly magnetized NSs, we should use a physically motivated initial model. A first approach is to consider that the hot, liquid initial phase lasts long enough to establish an MHD equilibrium, and reduce the pool of possible initial models to perfect MHD equilibria solutions, which have been calculated for different geometries (Colaiuda et al 2008;Ciolfi and Rezzolla 2013). Unfortunately, the formation of a NS during a supernova explosion is a very complex process, and the origin of strong magnetic fields in NSs and their presumed topology remains unclear. One of the most promising mechanisms at work is the magnetorotational instability, which, in the presence of differential rotation, can amplify exponentially fast a weak initial magnetic field in a proto-NS to a dynamically relevant strength (Guilet et al 2017). A recent simulation showed an effective dynamo (Mösta et al 2015), able to create an amplified large-scale toroidal magnetic field. A different mechanism is based on the interplay between compression and convection in the hot-bubble region between the proto-NS and the stalled shock (Obergaulinger et al 2015). Another relatively recent idea gaining popularity is the formation during a NS-NS merger (Ciolfi et al 2019). In any case, all the viable mechanisms involve some degree of turbulence, so that the outcome is plausibly different from a perfectly ordered dipole. With all these caveats in mind, one must choose to start the simulations, preferably with a few free parameters, to establish some qualitative trends in the long term evolution. For this reason, most of the existing works simply take a dipolar configuration, or at most some combination of dipolar and quadrupolar poloidal and toroidal components.
In the literature, there is a clear distinction between crustal-confined and corethreading fields. The former has been studied in-depth, with special focus on the Hall term. It assumes a type I superconducting core (not realistic) or, equivalently, that some other mechanism acts on very short timescales to expel most of the magnetic flux from the core. Core-threading magnetic fields are less studied, due to uncertain core physics. A popular configuration is the twisted-torus (e.g., Ciolfi and Rezzolla 2013): an MHD equilibrium solution with a large magnetic helicity where a dipole threads the core and the closed field lines within the star contain a toroidal field. We note the two fundamental differences between crustal-confined and core-threading configurations: first, the field curvature changes by one order of magnitude (roughly the size of the star versus the size of the crust); second, the location of most of the currents (core or crust) determines where Ohmic dissipation occurs, and the two regions have very different conductivities (see Fig. 7).
Regarding the temperature evolution, the initial conditions are much easier: it is well known that a few hours or days after birth, most of the star is nearly isothermal, so it is a good approximation to assume a constant temperature, between 10 9 to 10 10 K. Moreover, the particular choice of the initial temperature only affects the evolution in the first few days, which is completely irrelevant for following the NS evolution for thousands or millions of years. In any case, since the majority of existing works do not couple the magnetic field evolution to the temperature (most of them assume a constant temperature, or in some cases an independently prescribed function of time), we will begin by reviewing the main results of models with only magnetic evolution (2D and 3D), to conclude the section revising the only fully consistent magnetothermal simulations available (in 2D). 7.1 Magnetic field evolution in the neutron star crust 2D simulations. Many works (Pons and Geppert 2007;Viganò et al 2012Viganò et al , 2013Gourgouliatos et al 2013;Gourgouliatos and Cumming 2014b,a) have agreed on the general picture of the Hall-driven dynamics of a crustal-confined field in axial symmetry. For typical field strengths of 10 14 G, and starting from a predominantly poloidal dipolar field, we observe a stage dominated by the Hall drift (readjusting from initial conditions), which creates higher-order multipoles, followed by a quasi-stationary Ohmic stage. This structure, which has been called the Hall attractor (Gourgouliatos and Cumming 2014a), is characterized by a nearly constant angular velocity of the electron fluid (Ω ≈ j/en e r) along each poloidal field line, and proportional to the magnetic flux. This result holds even if the initial state is a high multipole, say l, with the system relaxing to a mixture of modes dominated by the l and l + 2, but again with the electron angular velocity linearly related to the flux (Gourgouliatos and Cumming 2014b). It is also relevant to remark that the Hall drift may noticeably accelerate the dissipation of magnetic fields, by continuously redistributing magnetic field energy towards smaller scales, where Ohmic dissipation is more effective. In the supplementary material, we provide the animations of two models with an initial dipolar poloidal magnetic field with surface polar intensity B p = 10 14 G plus a toroidal field with a maximum intensity B tor = 10 15 G. The models differ in the initial multipole of the toroidal field (l = 1 or l = 2).
3D simulations. Using a mixed spectral/finite-difference code, Wood and Hollerbach (2015); Gourgouliatos et al (2016) presented the first 3D simulations of crustalconfined fields, with an exterior boundary condition consisting of a general potential solution. The temperature was not included in the simulations, and the resistivity and density profiles were some radial dependent analytical functions, fitted to mimic a realistic model at T = 10 8 K (Cumming et al 2004). These 3D studies show new dynamics and the creation of km-size magnetic structures persistent over long timescales. Even using initial axisymmetric conditions, the Hall instability breaks the symmetry and new 3D modes quickly grow. These have lengthscales of the order of the crust thickness. A typical model is shown in Fig. 17. The surface field is highly irregular, with small regions in which the magnetic energy density exceeds by at least an order of magnitude the average surface value. By exploring many different initial models, Gourgouliatos et al (2016) found that magnetic instabilities can efficiently transfer energy to small scales, which in turn enhances Ohmic heating and powers the star persistent emission, confirming the 2D results.
More recently, Gourgouliatos and Hollerbach (2018) explored magnetic field configurations that lead to the formation of magnetic spots on the surface of NSs, extending previous 2D works (Geppert and Viganò 2014), as described in the final part of Fig. 17 Left: Magnetic field lines and magnetic energy density maps on the star surface (in colors), at t = 15 kyr, for an initial model consisting of an l = 1 poloidal field, and l = 2 toroidal field, plus a small non-axisymmetric perturbation. Right: Contour plot of the azimuthal component of the magnetic field at r = 0.995R , with R being the star radius, for the same model. Figures courtesy of Gourgouliatos et al (2016). §5.6. They show how an ultra-strong toroidal component is essential for the generation of a single spot, possibly displaced from the dipole axis, which can survive on very long timescales. We must note that boundary conditions arguably play a very important role to determine the scale of the, initially unstable, dominant modes since the thickness of the crust sets a preferred scale in crustal-confined models.

Coupled magneto-thermal simulations
We now turn to the complete problem: solving the temperature evolution coupled to the induction equation with realistic microphysics. To our knowledge, the only existing work studying the fully coupled magneto-thermal evolution of a realistic NS was , where they presented the results of 2D simulations. This work also re-analysed in a consistent way the available data on isolated, thermally emitting NSs (a sample of 40 sources), and compared the theoretical models to the data, concluding that the evolutionary models can explain the phenomenological diversity of isolated NSs by only varying their initial magnetic field, NS mass, and envelope composition.
As an example, in Fig. 18 we show three snapshots of the evolution of a crustal confined model, initially an l = 1 poloidal field with B p = 10 14 G (labelled as model A14 in Viganò et al 2013). Many of the general features described in previous more simple cases are also visible in this realistic model. Let us recap the most important details: • The first effect of the Hall term in the induction equation is to couple the poloidal and toroidal components so that, even if the latter is zero at the beginning, it is quickly created. After ∼ 10 3 yr, a quadrupolar toroidal magnetic field with Fig. 18 Snapshots of the magneto-thermal evolution of a NS model at 10 3 , 10 4 , 10 5 yr, from left to right. Top panels: the left hemisphere shows in color scale the surface temperature, while the right hemisphere displays the magnetic configuration in the crust. Black lines are the projections of the poloidal field lines and the color scale indicates the toroidal magnetic field intensity (yellow: positive, red: negative). Middle panels: intensity of currents; the color scale indicates J 2 /c 2 , in units of (G/km) 2 . Bottom panels: temperature map inside the star. In all panels, the thickness of the crust has been enlarged by a factor of 4 for visualization purposes. a maximum strength of the same order of the poloidal magnetic field has been created, with B ϕ being negative in the northern hemisphere and positive in the southern hemisphere. • Thereafter, under the effect of the Hall drift, the toroidal magnetic field rules the evolution, dragging the currents into the inner crust (see middle panels), and compressing the magnetic field lines. The Hall term is thus responsible for the energy redistribution from the large scale dipole to small scales (higher-order multipoles are locally very strong), possibly creating current sheets in some situations (here, in the equator).
• Where sufficiently small-scale components are present, the locally enhanced ohmic dissipation balances the effect of the Hall drift and a quasi-stationary state (resembling the Hall attractor) is reached. After ∼ 10 5 yr, the toroidal magnetic field is mostly contained in the inner crust. • We note that, at this point, most of the current circulates close to the crust/core interface. Therefore, the dissipation of magnetic energy is regulated by the resis- tivity in this precise region. In the model, there was a highly resistive layer in the nuclear pasta region leading to a rapid decay of the magnetic field, which has a direct imprint on the observable rotational properties of X-ray pulsars . • Joule heating modifies the map of the internal temperature. We can observe in the bottom panels of Fig. 18 how, at t = 10 3 yr, the equator is hotter than the poles by a factor of 3. This is caused by the insulating effect of the strong magnetic field discussed in §2.4. The presence of strong tangential components (B θ and B ϕ ) insulates the surface against the interior. In a dipolar geometry, the magnetic field is nearly radial at the poles, which remain thermally connected with the interior, while the equatorial region is insulated by tangential magnetic field lines. This has a two-fold effect: if the core is warmer than the crust, the polar regions will be warmer than the equator; however, if ohmic dissipation heats the equatorial regions, the situation is reverted. The temperature reflects the geometry of the poloidal magnetic field lines, which channel the heat flow.
In order to show more clearly the enhanced dissipation caused by the combined action of Hall and Ohmic terms, in Fig. 19 we show the evolution of the total magnetic energy stored in each component, comparing the evolution of the previous model with another model with the same initial data but switching off the Hall term (purely resistive case). In this case, there is no creation of a toroidal magnetic field or smaller scales. When the Hall term is included, ∼ 99% of the initial magnetic energy is dissi- Fig. 20 Evolutionary tracks in the P −Ṗ diagram of a typical NS with mass of 1.4M and radius of 11.6 km, with different initial magnetic field strengths: B 0 p = 3 × 10 12 , 10 13 , 3 × 10 13 , 10 14 , 3 × 10 14 , 10 15 G, evolving under the action of the Hall drift and Ohmic dissipation. Asterisks indicate the points when the star reaches the age of t = 10 3 , 10 4 , 10 5 , 5 × 10 5 yr. Dashed lines show the tracks followed without considering magnetic field decay. The figure includes the sample of X-ray pulsars with thermal emission analysed in Viganò et al (2013), which contains magnetars (MAG), nearby X-ray isolated NSs (XINS), rotation powered pulsars (RPP), and high magnetic field pulsars (HB). Figure courtesy of Viganò et al (2013). pated in the first ∼ 10 6 yr, compared to only the 60% in the purely resistive case. At the same time, a ∼ 10% of the initial energy is transferred to the toroidal component in 10 5 yr, before it begins to decrease. Note that the poloidal magnetic field, after 10 5 yr, is dissipated faster than the toroidal magnetic field. The poloidal magnetic field is supported by toroidal currents concentrated in the inner, equatorial regions of the crust. Here the resistivity is high for two reasons: the effect of the nuclear pasta phase, and the higher temperature (see right bottom panel of Fig. 18). Conversely, the toroidal magnetic field is supported by larger loops of poloidal currents that circulate in higher latitude and outer regions, where the resistivity is lower. As a result, at late times most of the magnetic energy is stored in the toroidal magnetic field. This example is very illustrative of the importance of knowing in detail the topology of the field and the location of currents at different stages. We refer the interested reader to Viganò et al (2013) for an extended analysis of different models, and how the initial magnetic field configuration affects the evolution. The qualitative behavior is similar to that shown in model A14, but subtle differences can arise when the strength or geometry of the initial field is modified.
Finally, we turn our attention to the rotational evolution of NSs. In Fig. 20 we show evolutionary tracks in the P −Ṗ diagram for a typical NS of 1.4 M with dif-ferent initial values of the initial magnetic field strength. The magnetic field configuration employed is the type A geometry (crustal confined) in Viganò et al (2013). Dashed lines show the results for models assuming a constant magnetic field, which are straight lines in the diagram. The solid lines, which account for realistic field evolution, show significant differences from the constant field models. Initially, the tracks overlap (B p is almost constant during an initial epoch, t 10 3 − 10 5 yr, which depends on the initial B 0 p ), but eventually, the field dissipates faster than the spin period evolution timescale and the lines bend down, at nearly constant P. This effect has been proposed to be the main reason for the observed clustering of periods of isolated X-ray pulsars . The particular value of the limit period mainly depends on the initial magnetic field and the resistivity at the crust-core interface. The large differences between the spin period evolution for models assuming a constant magnetic field and more realistic models make evident that the coupling between temperature, magnetic, and rotational evolution has to be considered.

Future prospects
After having reviewed the status-of-the art of the field of the long-term magnetothermal evolution of NSs, we highlight the three main areas that, in our opinion, need the focus of the researchers for the near-and mid-term future: • Although the numerical solution of the 3D heat equation is a well-studied problem in the literature of numerical methods, in our particular context of NSs, the first paper implementing a full 3D temperature evolution with realistic microphysics is yet to come. The main reason for this lack of models is that the study of the 3D temperature evolution alone does not add much to the problem, and only its coupling with the magnetic field evolution is of great interest. From the observational point of view, the possible existence of small-scale hotspots associated with the properties of the X-ray spectra is plausibly connected to the creation of small magnetic structures, and localized heat deposition. Thus, a necessary future step is to implement consistent 3D temperature evolution in the few existing 3D magnetic field evolution codes. • As discussed in Sect. 6, another open issue is the correct implementation of realistic, more general, boundary conditions at the star surface. Going beyond the popular and simple potential/vacuum solution seems a necessity, that has only begun to be considered. The magnetar observational data favor the presence of twisted magnetospheres, that can influence the interior dynamics in a significant manner. Among the different possible solutions, some mentioned in the text (solving elliptic equations, extended domains), it is unclear which has a better balance between physical motivation and computational cost. This point is of particular relevance when connected to the first one: the creation of localized hot spots may be strongly dependent on the applied boundary condition, because currents passing through the envelope may have the key to understand the very high temperatures of magnetars. • Finally, the core evolution is arguably the less explored part of the problem. Concerning the temperature, the core is almost isothermal to a very good approxi-mation. But the complexities of the interaction between superfluid neutrons and superconducting protons result in uncertainties of many orders of magnitude in the transport coefficients that determine the magnetic field evolution. A full 3D study of the evolution of the field penetrating the core and including all relevant physics does not exist, and it should be a high priority task for the incoming years.
All these efforts, in combination with the continuous upgrades of the microphysics ingredients, and the improving quality of the observational data with the new instruments, will allow us to decipher some of the fascinating physical processes taking place in the interiors of NSs.

A Poloidal-toroidal decomposition of the magnetic field
Any three-dimensional, solenoidal vector field B, can be expressed in terms of its poloidal and toroidal components In the literature, one can find different formalisms and notations to describe the two components. In this appendix we go through some of the ideas of the mathematical formalism and compare the most common notations. Adopting the notation of Geppert and Wiebicke (1991), the magnetic field can be written in terms of two scalar functions Φ(r,t) and Ψ (r,t) (analogous to the stream functions in hydrodynamics) as follows: where k is an arbitrary vector. This decomposition is particularly useful in situations where k is taken to be normal to one of the physical boundaries. Therefore, for a spherical domain, and using spherical coordinates (r, θ , ϕ), a suitable choice is k = r. In this case, ∇ × r = 0, and we can write: B tor = ∇Ψ × r .
Generally speaking, the radial component of the magnetic field is included in the poloidal part, while the θ and ϕ components are shared between poloidal and toroidal components. In axial symmetry, Φ = Φ(r, θ ) and Ψ = Ψ (r, θ ), the expressions are further simplified: the toroidal magnetic field is directed along the azimuthal direction ϕ. In this case the ϕ−component of the potential vector is given by A ϕ = −r × ∇Φ , and the poloidal field can be directly derived from B pol = ∇ × A ϕ .

B Potential solutions with Green's method
For potential configurations, we can express the potential magnetic field in terms of the magnetostatic potential χ m , so that The second Green's identity, applied to a volume enclosed by a surface S, relates the magnetostatic potential χ m with a Green's function G (see Eq. (1.42) of Jackson 1991): wheren is the normal to the surface. Comparing with the electrostatic problem, we see that no volume integral is present, because ∇ · B ≡ ∇ 2 χ m = 0. Note also that the factor 2π appears instead of the canonical 4π, because inside the star Eq. (83) does not hold, thus 2π is the solid angle seen from the surface. The Green's function has to satisfy ∇ 2 G(r, r ) = −2πδ (r − r ). The functional form of G is gauge dependent: given a Green's function G, any function F(r, r ) which satisfied ∇ 2 F = 0 can be used to build a new Green's functionG = G + F. The boundary conditions determine which gauge is more appropriate for a specific problem. In our case the volume is the outer space, S is a spherical boundary of radius R (e.g., the surface of the star), andn = −r . We face a von Neumann boundary condition problem, because we know the form of the radial magnetic field In order to reconstruct the form of we have to solve the following integral equation for χ m : 2π χ m (r) = R 2 π 0 2π 0 ∂ G ∂ r (r, r )χ m (R, θ ) sin θ dϕ dθ + − π 0 2π 0 G(r, r )B r (θ ) sin θ dϕ dθ .
For numerical purposes, we can express Eq. (94) in matrix form, introducing f i j = f (θ i , θ j ) evaluated on two grids with vectors θ i , θ j , with m steps ∆ θ . The coefficients of the matrix f i j are purely geometrical, therefore they are evaluated only once, at the beginning. The grid θ i coincides with the locations of B r (R, θ ), while the resolution of the grid θ j is M times the resolution of the grid θ i (M 5) to improve the accuracy of the integral function f i j near the singularities θ i → θ j . The resolution of the grid of ϕ k barely affects the result, provided that it avoids the singularities ϕ = 0, π/2. We typically use M = 10 and n ϕ = 1000. The calculation of the factors f i j is performed just once and stored. The matrix form is: From this, we obtain B θ by taking the finite difference derivative of χ m (θ ).