3D modelling of magneto-thermal evolution of neutron stars: method and test cases
Davide De Grandis, Roberto Turolla, Toby S. Wood, Silvia Zane, Roberto Taverna, Konstantinos N. Gourgouliatos
DDraft version September 10, 2020
Typeset using L A TEX twocolumn style in AASTeX63
3D modelling of magneto-thermal evolution of neutron stars: method and test cases
Davide De Grandis,
1, 2
Roberto Turolla,
1, 2
Toby S. Wood, Silvia Zane, Roberto Taverna,
4, 1 andKonstantinos N. Gourgouliatos
5, 6 Department of Physics and Astronomy, University of Padova, via Marzolo 8, I-35131 Padova, Italy Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Surrey, RH5 6NT, United Kingdom School of Mathematics and Statistics, Newcastle University, Newcastle upon Tyne, NE1 7RU, United Kingdom Department of Mathematics and Physics, University of Roma Tre, via della Vasca Navale 84, I-00146 Roma, Italy Department of Physics, University of Patras, Patras 26500, Greece Department of Mathematical Sciences, Durham University, Durham DH1 3LE, United Kingdom (Received September 10, 2020; Revised September 10, 2020; Accepted September 10, 2020)
Submitted to ApJABSTRACTNeutron stars harbour extremely strong magnetic fields within their solid outer crust. The topologyof this field strongly influences the surface temperature distribution, and hence the star’s observationalproperties. In this work, we present the first realistic simulations of the coupled crustal magneto-thermal evolution of isolated neutron stars in three dimensions with account for neutrino emission,obtained with the pseudo-spectral code parody . We investigate both the secular evolution, espe-cially in connection with the onset of instabilities during the Hall phase, and the short-term evolutionfollowing episodes of localised energy injection. Simulations show that a resistive tearing instabilitydevelops in about a Hall time if the initial toroidal field exceeds ≈ G. This leads to crustal failuresbecause of the huge magnetic stresses coupled with the local temperature enhancement produced bydissipation. Localised heat deposition in the crust results in the appearance of hot spots on the starsurface which can exhibit a variety of patterns. Since the transport properties are strongly influencedby the magnetic field, the hot regions tend to drift away and get deformed following the magnetic fieldlines while cooling. The shapes obtained with our simulations are reminiscent of those recently derivedfrom NICER X-ray observations of the millisecond pulsar PSR J0030+0451.
Keywords: neutron stars — magnetars — magnetohydrodynamical simulations — pulsars — stellarmagnetic fields INTRODUCTIONNeutron stars (NSs) are unanimously believed topower some of the most violent phenomena observedin the high-energy sky, from the hyper-energetic giantflares of ultra-magnetised NSs ( magnetars ; see e.g. Tur-olla et al. 2015; Kaspi & Beloborodov 2017, for reviews),to the spectacular merging of a binary NS system andthe associated emission of gravitational waves (Abbottet al. 2017). Despite this, many aspects of NS physicsare still poorly understood, mainly—but not only— con-
Corresponding author: Davide De [email protected] cerning their internal structure and composition, as wellas the topology of their magnetic field.Isolated NSs, from which (thermal) emission com-ing directly from star surface is visible in the X-ray-to-optical bands, provide an ideal laboratory to inves-tigate the physics of the interior of these objects, asfirst suggested by Tsuruta & Cameron (1966, see alsoTurolla 2009). NSs cool down as they age and theirthermal evolution is coupled to that of their magneticfield. Knowledge of the secular magneto-thermal evo-lution can discriminate between different cooling sce-narios when compared to observations, thus constrain-ing the equation of state of ultra-dense matter (see e.g.Page et al. 2006; Haensel et al. 2007). Moreover, it pro-vides a self-consistent map of the surface temperature, a r X i v : . [ a s t r o - ph . H E ] S e p De Grandis et al. which is essential in deriving any reliable estimate ofthe star radius from X-ray observations of (passively)cooling NSs (see e.g. Greif et al. 2019, and referencestherein). A detailed model of the short-term evolution,following an impulsive energy release in the NS surfacelayers, is equally desirable since it directly bears to theorigin of magnetar outbursts (see e.g. Rea & Esposito2011; Pons & Rea 2012; Coti Zelati et al. 2018) and ofthermal X-ray emission in radio pulsars (see e.g. Becker2009; Miller et al. 2019).Magneto-thermal evolution of NSs has been the focusof many investigations over the past decades (see Vigan`o2013, for a complete historical outline and further refer-ences). First attempts dealt with cooling in one dimen-sional (i.e., spherically symmetric) models with little orno account for the magnetic field (see e.g. Yakovlev &Urpin 1981; Page & Baron 1990). As a further step,axysimmetic, 2D, calculations were produced, but theseeither assumed a known evolution of the temperaturewhen solving for the magnetic field (Pons & Geppert2007) or the opposite (Aguilera et al. 2008). Moreover,inherent numerical difficulties prevented for a long timefrom including the Hall term in the induction equation,despite its importance in rearranging the magnetic fieldon the smaller spatial scales where dissipation is faster(Pons & Geppert 2007). The first consistent treatmentof the coupled magneto-thermal evolution in 2D waspresented in Vigan`o et al. (2013), who also succeededin coping with the Hall term. Recent efforts were de-voted to investigate the magnetic evolution with a fully3D approach and confirmed the role of the Hall term inshaping the magnetic field in the earlier stages of the NSevolution when a peculiar magnetic structure develops(the
Hall attractor ; Gourgouliatos et al. 2016).According to the commonly accepted picture, the coreof NSs is in a superfluid and superconducting state, forwhich the ground state is magnetic flux-free. Up to nowvery few investigations dealt with the magnetic evolu-tion including the core (see e.g. Ciolfi & Rezzolla 2013)and the structure of the magnetic field in the core of aNS is poorly understood as yet. Most studies of the mag-netic field evolution, both in two and three dimensions,were restricted to the NS crust, relying on the assump-tion that the Meissner effect has been able to expel anyflux from the core in a timescale shorter than those ofmagnetic and thermal evolution (see e.g. Lander 2014,but also Ho et al. 2017 for a different perspective). Inthis work the same approach is followed.Over the last few years X-ray (e.g. Bilous et al.2019) observations provided increasing evidence for thepresence of localised region(s) on the surface of dif-ferent classes of isolated NSs with non trivial ther- mal/magnetic properties and evolution. To explainthese observations, as well as to validate results obtainedin 2D calculations, fully coupled magneto-thermal 3Dsimulations are necessary. In this work we present someof the first simulations of such kind, showing some of thepossible applications in which a 3D treatment is neces-sary to fully tackle the observed phenomenology.The paper is organised as follows. In section 2 wepresent the basic equations and their numerical imple-mentation. In section 3 some study cases for the long-term evolution of NSs are presented; in particular, theonset of eMHD instabilities is discussed in section 3.2.Some examples of the short-term evolution following alocalised crustal heating are illustrated in section 4, witha view to applications to magnetar outbursts (section4.1) and to surface heating in pulsars (section 4.2). Dis-cussion and conclusions follow in section 5. THE MODEL2.1.
Input physics and evolution equations
The NS crust comprises a Coulomb lattice in whichnuclei have negligible motion. Hence, the crustal cur-rents are produced entirely by the flow of electrons,which form a highly relativistic and strongly degener-ate Fermi gas. Still, their mean velocity is typicallyonly a tiny fraction of the speed of light. We can there-fore resort to the (non-relativistic) electron magneto-hydrodynamics (eMHD) approximation in treating thecrustal dynamics. The evolution of the magnetic field B in the crust is described by the induction equation that,taking into account also the effects of thermal coupling,can be written in the form ∂ B ∂t = − c ∇ × (cid:2) σ − · J + G · ∇ T − ∇ µ/ e (cid:3) (1)where the term in square brackets is the electric field E as given by generalised Ohm’s law. Here, c is the speedof light, e the electron charge, σ and G are the electricconductivity and the thermopower tensors, and µ is theelectron chemical potential. The latter, for a degener-ate relativistic Fermi gas, depends only on the density, µ = c (cid:125) (3 π n ) / where (cid:125) is the reduced Planck con-stant. Neglecting the displacement current, the electroncurrent is given by J = c ∇ × B / π .Assuming the temperature of the crust to be well be-low the electron degeneracy temperature, but above theion plasma temperature, scattering of electrons can bedescribed in terms of an energy-dependent relaxationtime τ (e.g. Ziman 1972; Urpin & Yakovlev 1980), andthe electron conductivity can then be approximated as( σ − ) ij = σ − δ ij + (cid:15) ijk B k e cn (2) D magneto-thermal evolution of neutron stars σ = e c nτ ( µ ) µ (3)and the anti-symmetric part represents the Hall effect.The thermopower can be calculated using the Mottformula, and in general has an isotropic part and a partproportional to the conductivity tensor. In our modelwe include only the isotropic part (the so-called Seebeckterm ), which is responsible for the Biermann batteryeffect, G = e ∂ σ − ∂µ (cid:12)(cid:12)(cid:12)(cid:12) T · k (cid:39) − π k B T e µ δ ij (4)where the approximate equality is obtained further as-suming electrons to form perfect Fermi gas. Here, k isthe thermal conductivity tensor, that is taken to be pro-portional to σ , according to the Wiedemann-Franz law k = π k T σ (5)where k B is Boltzmann’s constant.The evolution of temperature is, in turn, governed bythe heat equation C V ∂T∂t = − ∇ · (cid:16) T G · J − k · ∇ T − µ e J (cid:17) + E · J + N ν (6)where C V is the heat capacity (per unit volume) of thecrust, N ν is the neutrino emissivity due to weak pro-cesses, and the term in round brackets is the electronenergy flux.Although general-relativistic effects can be accountedfor with no inherent difficulty in equations (1) and (6),see e.g. Pons et al. (2009) and Vigan`o et al. (2013), theyare not included here. The reason for this is twofold.First, given the small thickness of the crust, they areof limited importance and will not change our resultsqualitatively and, second, a proper general-relativistictreatment impacts on the boundary conditions whichare imposed on the evolution equations (see section 2.2).While this poses no serious problem in 2D, it becomesquite troublesome in 3D. Equations (1) is analogous tothe one solved in Gourgouliatos & Cumming (2014a)(with the addition of thermocoupling terms) and Vigan`oet al. (2013), where (a different version of) equation 6was included as well.In order to simplify the equations, all physical quanti-ties are henceforth expressed in terms of values typical ofthe outer crust in a magnetar. In particular, the temper-ature, magnetic field, relaxation time and chemical po-tential are normalised to T = 10 K, B = 10 G, µ =2 . × − erg, and τ = 9 . × − s. This implies thatthe reference values for the number density n and the conductivity η = c / (4 πσ ) are n (cid:39) . × cm − and η (cid:39) . × − cm s − . It is also useful to introducefour length scales, which are relevant to the electron dy-namics, λ = (cid:18) k B T πn e (cid:19) / Debye length (7) d = (cid:18) µ πn e (cid:19) / skin-depth (8) L = µ e B Larmor radius (9) l = cτ mean free path, (10)supplemented by the star radius R (cid:63) = 10 km as themacroscopic length scale. Furthermore, the ohmic time τ O = R (cid:63) /η (cid:39) × yr is taken as the reference timescale, and C V = k B n as the scale for the heat capacity.The evolution equations to be solved for B and T be-come then ∂ B ∂t = Se ∇ (cid:18) µ (cid:19) × ∇ T + Ha ∇ × (cid:20) µ B × ( ∇ × B ) (cid:21) + − ∇ × (cid:20) τ µ ∇ × B (cid:21) (11)1Ro C V T ∂T ∂t = ∇ · (cid:0) τ µ χ · ∇ T (cid:1) + PeSe | ∇ × B | τ µ ++ Pe µ ( ∇ × B ) · ∇ (cid:18) T µ (cid:19) + 1Ro N ν (12)where we defined χ ij = δ ij + Ha ( τ /µ ) B i B j − Ha( τ /µ ) (cid:15) ijk B k ( τ /µ ) | B | (13)and introduced the adimensional numbersHa = l/L (cid:39) Hall (14)Se = π Llλ d (cid:39) . Seebeck (15)Pe = 3 d Ll (cid:39) × − Peclet (16)1Ro = (cid:114) π PeSe Ha (cid:39) × − Roberts . (17)Performing adimensionalisation with these scales, theneutrino emissivity is expressed in units of N ν =8 π e c R τ /µ k b T = 1 . × erg s − cm − .The large thermal conductivity of the crust, which isreflected in the fact that Ro (cid:29)
1, means that on thetimescale of magnetic evolution the term dependent onthe heat capacity of the crust is subdominant. Ratherthan using a detailed microphysical model for C V , we De Grandis et al. therefore simply take C V = τ µ T , which implies a con-stant effective thermal diffusivity throughout the crust.Under this assumption, equation (12) depends on tem-perature only through T , which proves to be an advan-tageous feature for numerical implementation.Under the eMHD approximation, electrons and pro-tons in the crust have equal, time-independent numberdensity n . We will assume that the crust is spheri-cally symmetric, so that n is a function of the radius r alone. With our definition of the chemical potential(see above), this implies that µ depends on r and wetake µ ( r ) = (cid:18) − r . (cid:19) / , (18)following Gourgouliatos & Cumming (2014a); µ ( r ) in-creases from unity at the outer boundary ( r = 1) to (cid:39) . r = 0 . τ is a function of r only and, in par-ticular, we take τ ≡ Q (cid:39) τ to be independent on temperature is adequate in thelower crust, while it is just an approximation in the up-per crust, where scattering is dominated by phonons(Potekhin et al. 2015). We, nevertheless, note that tak-ing τ ≡
1, the conductivity in the upper crust corre-sponds to the phonon conductivity at a realistic temper-ature, T ≈ K. We assume an Fe-Ni crust, withoutaccounting for chemical composition stratification.The emission of neutrinos in the crust is due to a largevariety of reactions. In this work, the four dominant con-tributions are taken into account, namely phonon decay,neutrino pair production, neutrino bremsstrahlung andneutrino synchrotron emission N ν ( n, T, B ) = N ph ( n, T ) + N pair ( n, T )+ (19) N bre ( n, T ) + N syn ( n, T, B ) . We make reference to Ofengeim et al. (2014) for neutrinopair bremsstrahlung decay and to Kantor & Gusakov(2007) for phonon decay. A complete review can befound in Yakovlev et al. (2001). These papers providefitting formulae for numerical evaluation, that were im-plemented in our code.2.2.
Boundary conditions
Solution of the evolution equations (1) and (6) requiresboundary conditions which reflect a number of physi-cal prescriptions at the core-crust interface and at thestar surface. We assume that all magnetic flux has beenexpelled from the superconducting core. This requires that the normal component of the magnetic field and thetangential component of the electric field must vanish atthe core-crust boundary, r = r c . The latter results in anonlinear boundary condition for the magnetic field dueto the presence of the Hall term. Nevertheless, this con-tribution is negligible near the bottom of the crust dueto the high electron density (see Hollerbach & R¨udiger2004), and can hence be neglected. This allows to writethe boundary conditions in terms of the radial magneticfield and of the tangential component of the current J t , B r ( r c ) = 0 J t ( r c ) = 0 . (20)We assume that the electrical conductivity in the mag-netosphere is negligible in comparison with that of thecrust, and therefore match the field at the crust outerboundary to a potential one. This can be achieved ina very natural way by exploiting the spectral nature ofour code, which decomposes the field using the sphericalharmonics Y m(cid:96) ( θ, φ ) as the basis (see section 2.3), afterintroducing a poloidal-toroidal decomposition B = ∇ × ∇ × ( r B pol ) + ∇ × ( r B tor ) . (21)In such a representation each mode of a potential field ispurely poloidal and such that ( B pol ) m(cid:96) ∝ r − ( (cid:96) +1) . Hence,the boundary condition can be met by requiring that ∂∂r ( B pol ) m(cid:96) + (cid:96) + 1 r ( B pol ) m(cid:96) (cid:12)(cid:12)(cid:12)(cid:12) r = R (cid:63) = 0;( B tor ) m(cid:96) | r = R (cid:63) = 0 . (22)The core conducts heat even more efficiently than thecrust, and is therefore approximately isothermal. Itcools by neutrino emission according to the equation ∂T c ∂t = − N c ( T ) C c (23)where C c is the core specific heat and N c ( T ) = N T k the neutrino emissivity of the core. The star long-termthermal evolution is governed by eq. (23) once the neu-trino emissivity is specified. Our model uses a standardslow-cooling scenario with k = 8, C c = 10 erg / s / K , N = 10 erg / c m / K (Page et al. 2004).The surface temperature is controlled by the proper-ties of the thermal blanketing envelope. This layer isgeometrically very thin, but hosts a large temperaturegradient. Thus, the widespread approach is to treat itseparately, using a plane-parallel approximation to ob-tain a relation between the temperature at the bottomof the envelope T b (that is, the temperature of the top ofthe crust) and the surface temperature T s (Gudmunds-son et al. 1983). Assuming that no energy gains or losses D magneto-thermal evolution of neutron stars − ( k · ∇ T ) · ˆ r = σ sb T s ( T b , B ) , (24)where the left-hand side is evaluated at the top of thecrust and ˆ r is the radial unit vector. We have chosenthe form T s ( T b , g, B ) = T (0) s ( T b , g ) X ( T b , B ) (25)where g is the gravitational acceleration at the surface.We used the expressions for T (0) s as calculated in Gud-mundsson et al. (1983) for an iron envelope neglectingmagnetic fields, and the magnetic correction X ( T b , B )obtained in Potekhin & Yakovlev (2001).2.3. Numerical implementation
Equations (11) and (12) were solved in three dimen-sions using a suitably modified version of the code par-ody , which was originally developed by Dormy et al.(1998) and Aubert et al. (2008). A version of the samecode, which did not include the thermo-magnetic cou-pling, was first used to investigate the magnetic field evo-lution in NSs in Wood & Hollerbach (2015). The code ispseudo-spectral: it uses a finite grid in the radial direc-tion and an expansion in spherical harmonics Y m(cid:96) ( θ, φ )for the angular part. The NS crust is assumed to bea perfect spherical shell. The time-stepping algorithmis Crank-Nicholson for the Ohmic diffusion, backward-Euler for the isotropic part of the thermal diffusion, andAdams-Bashforth for all other terms.We typically use 128 radial grid points, and spheri-cal harmonics up to degree (cid:96) ≈ (cid:46)
100 m on the surface. Parallelisation isimplemented using mpi and the code is run on a clus-ter of CPUs. Work is distributed in such a way thateach thread takes care of a spherical shell containing N r mod N cores points, where N r is the radial grid size. Inorder to compute space derivatives within our finite dif-ference scheme in each thread, a single shell should con-tain at least four grid points, hence to achieve the desiredresolution the code is typically run on (cid:46)
32 cores. STUDY CASESIn order to validate the code and provide comparisonswith previous works, we first address the problem ofthe secular magneto-thermal evolution of highly magne-tised, isolated NSs. The magnetic evolution follows twodifferent timescales, the Hall and Ohm ones (Goldreich& Reisenegger 1992), τ H = 4 πn eR cB ≈ yr (26) τ O = R /η ≈ yr . (27) t (yr)10 U ( c o d e un i t s ) solid: B pol (0) B tor (0)dash.: B pol (0) + B tor (0)dott.: B tor (0) = 0 = 1= 2= 3= 4= 5= 6= 7 Figure 1.
Time evolution of the energy in the first seven (cid:96) modes for three typical initial configurations of the field:a purely poloidal field (solid) and two cases with an addedtoroidal field of opposite polarity (dashed and dotted).
Magnetic field reconfiguration occurs on the Halltimescale, when small scale structures are formed bythe action of the Hall term, while on the Ohm one, dis-sipation takes place. Long-term thermal evolution alsooccurs on a time (cid:46) τ O (see e.g. Potekhin et al. 2015,for a review). A 3D approach is particularly suited, andindeed necessary, to follow the formation and evolutionof small-scale structures in the Hall phase.3.1. Neutron star magneto-thermal evolution
In order to set the initial magnetic configuration forour simulations, we followed the widespread approachof confining the field in simple, large-scale structures(R¨udiger et al. 2013). In particular, we selected a force-free B -field matching our boundary conditions, withboth non-zero poloidal ( (cid:96) = 1, m = 0) and toroidal ( (cid:96) =2, m = 0) components. For such a field, the componentsof equation (21) take the form B m(cid:96) ∝ Y m(cid:96) ( θ, φ ) ζ (cid:96) ( r ) /r ,where ζ (cid:96) is a linear combination of spherical Bessel func-tions of degree ± (cid:96) , constructed in such a way to obeythe boundary conditions (see Chandrasekhar & Kendall1957, for a full derivation). We stress that the evolutionof poloidal/toroidal components is strongly coupled bythe action of the Hall term, that can transfer energyboth ways between them (Pons & Geppert 2007). Theinitial temperature profile is assumed to be a constant, T ( r, t = 0) ≡ K, but we note that the overall evolu-tion is virtually independent on this choice. This is dueto the fact that the term ∝ ∂T /∂t in equation (12) is We remark that throughout the work the initial time is set in cor-respondence to the superfluid transition, which typically occursa few years later than the formation of the proto NS.
De Grandis et al. suppressed by a factor Ro − ≈ − , so that the tem-perature rapidly achieves a quasi-steady state.The evolution of the B-field over a few Hall timescalesis shown in Figure 1 for three different initial magneticconfigurations: a purely dipolar field ( B pol (0) ≈ G, B tor (0) = 0) and a field with poloidal and toroidal com-ponents of the same order but opposite relative po-larity ( B pol (0) ≈ ± B tor (0) ≈ G). Our simula-tions confirm the previous finding that the magneticfield evolves towards the so-called
Hall attractor (Gour-gouliatos & Cumming 2014b), where the magnetic fieldtends to reach a configuration dominated by the modes (cid:96) = 1 , , , , B r , B θ and B φ , at the beginning of the sim-ulation ( t = 0) and at t = 3 × yr (cid:39) τ H are shown infigure 2 for the case with B pol (0) ≈ + B tor (0). A gen-eral feature of the magnetic field is the appearance of anequatorial structure in which the field is stronger (Gour-gouliatos et al. 2016), and of small scale structures dueto the Hall term.As already mentioned, the temperature distributiontends to follow the magnetic field. The structure of theHall attractor, in which an equatorial current ring forms,is reflected in a hotter equatorial region. Moreover, for-mula (5) implies that heat tends to be transported pref-erentially along the field lines. Hence, the equatorialregion is hotter not only because of higher dissipation,but also because heat is trapped by the closed field linesappearing in that region. Figure 3 shows a typical case,that is representative—at least qualitatively—of all ournearly-axisymmetric runs. We note that, owing to thedependence of the properties of the heat blanketing en-velope on the geometry of the magnetic field, the ob-servable surface map can be quite different from the oneon the top of the crust. As an example, the last panelof Figure 3 shows the surface temperature for the verysame case: the overall topology is quite different, asthe equatorial belt is not just hotter, but instead showsa colder ring at the very equator (see e.g. the recentresults in Kondratyev et al. 2019, in which a similarbehaviour is discussed in a 3D stationary framework).Even though the various features are on a large scale,they exhibit a smaller scale –yet well resolved– structure,due to the Hall term. 3.2. Magnetars and eMHD instabilities
As already noted in Gourgouliatos & Pons (2020),the presence of a strong toroidal field can trigger a re-sistive tearing eMHD instability (Wood et al. 2014).This instability, even when starting from an initialcondition that is essentially symmetric, produces non-axisymmetric small-scale magnetic structures, that, dueto Joule dissipation, translate into localised heat depo-sition. A strong toroidal component in the star crustpasses undetected and is invoked to explain the observedactivity in the so-called low-B magnetars, i.e. sourceswith a dipole field comparable to that of the radio-pulsarpopulation (see section 5.1 for further details).To explore better this issue, we ran a simulation as-suming an (cid:96) = 1, m = 0 initial magnetic field witha poloidal field B pol (0) ≈ G and a toroidal one B tor ≈ × G. Given the nature of the solutionwe are looking for, the resolution for this case was im-proved to (cid:96) max = 250, corresponding to cells of a fewtenths of meters on the surface. Indeed, an instabilityis triggered after about a Hall time τ H . The spectrumof all the (cid:96) modes at t (cid:39) yr (Figure 4) exhibits thecharacteristic features of the Hall attractor: even modesare suppressed with respect to the nearby odd ones upto (cid:96) (cid:46)
100 and this produces the typical wavy profile.However, the onset of an instability is marked by theappearance of well-resolved structures that form up to (cid:96) (cid:39) (cid:96) guarantees that the instability is of physicaland not numerical origin. The slight increase at veryhigh (cid:96) is due to numerical aliasing. The spectrum of m modes, on the other hand, is sharply peaked towardszero and hence is not shown.In the small structures, the magnetic field can reachvalues up to ∼ × G and this drives a local tem-perature increase, as shown in figure 5. Such strongfields generate high magnetic stresses in the crust. Asa gauge to determine whether such stresses are strongenough to lead to a crustal failure, we compared them tothe maximum mechanical yield of the crust through thevon Mises criterion (see e.g. Pons & Rea 2012; Lander& Gourgouliatos 2019), (cid:113) ¯ M ij ¯ M ij (cid:38) τ max ( n, T ) (28)where ¯ M ij is the traceless part of the magnetic stresstensor M ij = B i B j / π . Chugunov & Horowitz (2010)derived estimates for τ max by means of molecular dy-namic simulations and elucidated the strong dependenceof the breaking stress on temperature. In our calculationwe used the fit they provide for the maximum crustal D magneto-thermal evolution of neutron stars Figure 2.
A meridional cut of the crust along the prime meridian ( φ = 0) showing the the magnetic field components B r , B θ and B φ (from left to right) at the start (top row) and after 3 × yr ≈ τ H (bottom row) for the run with B pol (0) ∼ + B tor (0).The plots for the φ component also show the field lines of the poloidal field. Here and in all figures where relevant, the crustthickness is enhanced by a factor 4 for better visualisation. Figure 3.
Temperature maps at t = 3 × yr ≈ τ H for acase with B pol0 ≈ + B tor0 ≈ G. Top : meridional cut, withthe field lines of the poloidal component superimposed;
Bot-tom left : temperature at the top of the crust (i.e. under theheat blanketing envelope);
Bottom right : surface tempera-ture according to equation (25), showing how the envelopecan change the very topology of the temperature distribu-tion. yield, τ max = (cid:18) . − . − (cid:19) n i Z e a (29) ( c o d e un i t s ) Figure 4.
Power spectrum (in code units) of the (cid:96) modesat t (cid:39) yr. On top of the wavy profile privileging oddmodes, typical of the Hall attractor, a more complex patternat short wavelengths reflecting the instability is visible. Theslope (cid:96) − , obtained from scaling relations for Hall turbulencein Goldreich & Reisenegger (1992), is shown for reference. where Γ = Z e /ak B T is the classical Coulomb couplingparameter, n i is the ion density, a = (4 πn i / − / is theion sphere radius and, following Horowitz et al. (2007),we took Z = 29 . τ max decreases at higher T , a 3D, coupledmagneto-thermal code provides the most accurate wayto investigate the onset of crustal failures in magnetars.Figure 6 shows the ratio between the magnetic andthe breaking stress in our simulation after a time6 . × yr (cid:46) τ H . The map refers to the region of thecrust where the ratio is maximum, at about half the De Grandis et al.
Figure 5.
Temperature at the top of the crust showing theformation of a hotter equatorial belt with a small-scale, yetnumerically resolved, pattern that reflects the eMHD insta-bility. Note that this is not the surface temperature. crust depth. The magnetic stress reaches values up to ∼
50% of the maximum yield in our simulation so thatvon Mises criterion for crustal yielding is likely to befulfilled. Crustal failures can therefore be triggered inour simulation because of the large magnetic stress cou-pled with the heating produced by magnetic dissipationwhich significantly rises the temperature in the equato-rial region, thus lowering the breaking stress. The insta-bility lasts for some 1000 yr before it is damped by dissi-pation. This directly concerns magnetar activity, sincemagnetically induced crustal failures (“starquakes”) arethought to be responsible for magnetar bursts and out-bursts (Pons & Rea 2012).We conclude this section noticing that the use of thevon Mises criterion as expressed by equation (28), albeitwidespread in the literature, should be taken with somecare. In fact, it does not take into account the effectsof the enormous gravity, that tends to inhibit any radialdisplacement (see e.g. Haskell 2008). However, since theresistive tearing instability arises as a consequence of thepresence of a strong toroidal field, our result is not muchaffected even setting to zero all radial shear terms (themaximum stress-to-yield ratio decreases from ∼
50% to45%). Nevertheless, only a consistent, non-local calcu-lation which takes into account the global hydrostaticstructure of the crust, could unambiguously solve theissue. LOCALISED HEATING IN NS CRUSTIn order to fully exploit the three-dimensionality of ourcode, we investigated models in which a localised heatsource is present in the NS crust. This is accounted forby adding a term ˙ H/ V heated to equation (6), describing Figure 6.
Ratio between magnetic stresses and maximumyield during the instability, shown at the radius where isattains its maximum value, close to half of the crust depth. the heat injection rate per unit volume. In particular, weconsider two cases: (i) localised heat deposition in thedeep crustal layers, and (ii) heating of the star’s externallayers. Although no direct application to real astrophys-ical sources will be attempted, these two models are ofinterest in connection with the evolution of magnetaroutbursts and the X-ray emission from radio-pulsars.4.1.
Heating in the deep crust
As already mentioned, magnetar activity is believedto be associated with crustal failures (see Section 3.2).However, the crustal dynamics in such events is littleexplored as yet, owing to its inherent complexity (e.g.the crust may flow plastically, Lander 2016). Such astudy is beyond the capability of our code, which doesnot incorporate a description of the motion of crustalmatter.As a minimal model to address the physics of crustalfailures, and in particular the way in which heat is trans-ported to the surface, we therefore performed a simula-tion in which energy is injected during a short time in-terval in a localised region of the crust, much in the sameway as in Pons & Rea (2012) but exploiting our fully 3Dapproach. As the background state, we take a NS withan initial field B pol ≈ G and B tor ≈ G thathas been consistently evolved for a Hall time. This hightoroidal field configuration was chosen in the spirit ofthe results of section 3.2 and mimics a low-B magnetar.In our test model heat has been released in the north-ern hemisphere and in the innermost half of the crust,assuming a gaussian profile along the three spatial di-mensions with σ r (cid:39)
100 m, σ θ (cid:39) σ φ (cid:39) π/ H (cid:39) × erg s − , modulated by the gaussian profile.Heating is assumed to be quasi-instantaneous (the ˙ H term is active for ∆ t inj (cid:39) T s from equation(25). The time evolution follows a typical FRED (fast- D magneto-thermal evolution of neutron stars t (yr)10 l o g L t o t ( e r g / s ) Neutrino-onNeutrino-off
Figure 7.
Luminosity evolution after an impulsive heat in-jection in the inner half of the crust. Neutrino emission re-duces the peak luminosity by a factor ∼ . raise-exponential-decay) pattern, as shown in Fig. 7.The two curves in Fig. 7 illustrate the role played byneutrino losses in the crust. The temperature rise pro-duced by heat deposition, in fact, is large enough in thiscase to make neutrino emission sizable (contrary to whatoccurs when the crust is not heated) and this results ina photon luminosity lower by a factor ∼ ∼ erg s − ,with an increase of a factor ≈
10 above the quiescentlevel.The hot structure that develops onto the surface ex-hibits a somehow peculiar evolution. In fact, its shapeis determined by heat diffusion, which is not isotropicbut depends on the magnetic field direction according toequations (2) and (5). Hence, heat tends to flow alongfield lines. Figure 8 shows how during the luminosityrise time heat is not just flowing radially to reach thesurface, but does so following the magnetic field. More-over, once it is formed the hotter region tends to driftas it cools down, both in latitude, towards the equa-tor, and in longitude. Such behaviour is clearly visiblein Figure 9 which shows four snapshots of the heatedsurface patch evolution.The duration of this event is of some thousands ofyears, with a rise time of about a century (however,timescales in this case are affected by code limitations,see section 5.2). Still, on a qualitative basis, it can betaken as a representation of the observed flux variationduring magnetar outbursts, which happens on shortertimescales (Coti Zelati et al. 2018).4.2.
Surface heating
The framework discussed in Section 4.1 can be easilyadapted to study the shape of pulsar hot spots. In fact,
Figure 8.
Meridional cuts (at the same φ ) of the evolutionof the hot spot during the rise phase. The first panel corre-sponds to the initial injection, and the subsequent ones areseparated by ∼
50 yr. Transport of heat to the surface hap-pens preferentially along magnetic field lines, whose planarprojection is superimposed in black. Note that colour barrange decreases between the two rows to improve visualisa-tion.
Figure 9.
Surface thermal evolution of the hot spot pro-ducing the luminosity shown in Fig. 7. Time increases fromleft to right and from top to bottom; snapshots are separatedby ∼
200 yr and the first one corresponds to the peak of theluminosity curve. The magnetic north pole is highlighted forreference. the physical ingredients remain the same as long as itcan be assumed that no other effects apart from heatingcome into play. The major difference with respect tothe model presented in the previous section is that nowenergy is deposited in the outermost crustal layers, as0
De Grandis et al.
Figure 10.
The extrapolated external magnetic field lines,in which the colour bar encodes the strength from black(zero) to copper (maximum value), compared to a purelydipolar field (dashed red lines). Given the high degree ofsymmetry of this case, only a quarter of the star is shown forbetter visualisation. it is the case e.g. for the heat deposited by backflowingcurrents on the surface of radio-pulsars.As a background state, we take a NS with initial field B pol ≈ B tor (cid:39) G and temperature T (cid:39) × K,which was evolved for some Hall times, t ∼ yr. Then,a heating source is placed in a small region of size ∼ . Backflowing currents in pulsars can reach a depthranging from about a tenth to the entire width of thecrust (Karageorgopoulos et al. 2019). In our simula-tions we choose to insert the heat source uniformly fromthe surface down to a quarter of the crust width. Weat any rate checked that different depth values providequite similar results, possibly because this length is any-way much smaller than the other relevant lengths in theproblem.Results show that if heat injection is steady, the hotspot reaches a state of quasi-equilibrium in a few years.Starting with a heated patch of size ∼ Even if our simulations do not explicitly include the dynamics ofthe magnetosphere, its configuration can be extrapolated fromthe boundary condition 22 for the magnetic field at the topof the crust, requiring that for each harmonic B m(cid:96) ( r > R (cid:63) ) = B m(cid:96) ( R (cid:63) ) /r (cid:96) +1 . Figure 11.
Top row: initial (left) and equilibrium (right)stages of the evolution of a hot region (magnetic pole atthe centre, towards the observer) during steady heating fromabove.
Bottom row: cooling of the spot after heating isturned off. The last three snapshots (from top to bottom andfrom left to right) are separated by a time interval ∼
200 yr.Note that the colour bar range decreases between panels tohighlight the effect. heat injection ˙ H = 5 × erg s − . We then followedthe evolution of the same spot after the heating term isturned off. In a time ≈ H =5 × erg s − , we observe a similar phenomenology, D magneto-thermal evolution of neutron stars ≈ Glarge-scale (quasi) dipolar field, where the field strengthcan reach values up to 6 × G. Thus, localised heat-ing may also account for small-scale magnetic struc-tures in the crust, that are therefore not originated bydynamo-like processes.The appearance of these crustal magnetic features re-flects in the creation of local magnetic structures in themagnetosphere, even though the overall B-field remainsvery close to dipolar. Figure 14 shows the external fieldlines for the same case as in Figure 13, as derived by solv-ing the magnetospheric structure (see footnote 2 at page10). After subtracting the contribution of the m = 0modes, which are dominated by the dipolar field, a smallmagnetic field loop is clearly visible above the heatedregion, extending outwards up to a distance (cid:46) R (cid:63) witha typical strength ≈ G. This shows that magneticstructures are not necessarily confined to the crust, butcan extend in the inner magnetosphere.The cooling phase lasts some thousands years. It istherefore possible that the aftermath of powerful heatingevents can produce long-lasting thermal structure on aNS crust, evolving in complex patterns along field lines. DISCUSSION AND CONCLUSIONSIn this paper we presented for the first time 3D numer-ical simulations of the coupled magneto-thermal evolu-tion in isolated neutron stars with full account for neu-trino emission from the crust and a simplified neutrinocore cooling model. While results for the long-term evo-lution show no substantial deviations with respect tothose obtained with 1D and 2D calculations (see e.g.Pons & Vigan 2019, for a review), the capability of a3D approach to consistently deal also with the smallerspatial scales proved essential to highlight the onset ofeMHD instabilities and to follow the evolution of lo-
Figure 12.
Cooling of a spot similar to the one in Figure11 but for a heat flux 10 times higher. Here snapshots areseparated by ∼
300 yr, the first one referring to the time atwhich heating stops. Note that colour bar range decreasesbetween the two rows to highlight the effect.
Figure 13.
Radial component of the magnetic field in thelocalised heating steady state shown on the first panel offigure 12. The crust thickness is enhanced by a factor 4 tohelp visualisation. calised heat injection in the star crust. In particular,out main findings are:(i) the magnetic field evolves towards the so calledHall attractor (Gourgouliatos & Cumming 2014b).In this configuration magnetic energy is storedpreferentially in the odd modes and especially inthe (cid:96) = 1 , , , ≈ G–10 G) can trigger the resistive tearingeMHD instability in less than a Hall time. Oursimulations show that the appearance of small-2
De Grandis et al.
Figure 14.
The extrapolated external magnetic field for thecase in Figure 13. The left half shows the total field and theright one the difference between the total field and its m = 0modes, which are dominated by the dipole. scale high-B structures, mainly along the equator,coupled with a local enhancement of the temper-ature produce the conditions for crustal yieldingaccording to the von Mises criterion;(iii) a localized, impulsive heat injection in the deepcrustal layers results in a cooling hot spot on thestar surface. The emitted luminosity has a sharprise followed by an longer decay;(iv) as a result of anisotropic heat transport in themagnetized crust, the heated region drifts and maychange its shape as it cools;(v) even with an essentially dipolar field, quasi sym-metric hot regions near the poles can cool downassuming a crescent-like shape.Our 3D simulations of the evolution of a locally heatedregion in the star crust revealed a variety of behavioursreflecting the location of the heat source (position anddepth in the crust), the energy injection rate and thecrust magnetic and thermal structure. In particular weconsidered two scenarios, in which heating occurs eitherinside the crust (deep heating) or in the outermost lay-ers (surface heating). Both of them may be relevant formagnetar outbursts, during which a hotter region on thestar surface appears and then progressively cools downand shrinks (see e.g. Coti Zelati et al. 2018). In fact, thishas been explained in terms of dissipated magnetic en-ergy inside the crust (Lyubarsky et al. 2002; Pons & Rea2012) or of Joule heating due to returning currents flow-ing along the field lines of a (locally) twisted magneticfield (Beloborodov 2009, see also Turolla et al. 2015).In the simulations presented in this work, neutrinoemission is relevant only for the case presented in sec-tion 4.1. In fact, neutrinos become important if injectionis fast, so that high local temperatures can be reached. According to 2D simulations Pons & Rea (2012), largeneutrino losses result in an upper limit on the radiativeluminosity released in magnetar outbursts. At presentsuch regimes can not be investigated with our code, dueto numerical hindrances associated with the treatmentin three dimensions of a strongly non linear term (asa rule of thumb, N ν ∝ T . ). In fact, this term cancause the appearance of numerical spurious features inour solutions when temperature gradients become veryhigh. Moreover, if the background state has an ultra-strong magnetic field, turbulent patterns analogous tothose discussed in section 4.2 can be triggered also in acontext of impulsive heat injection; with our present nu-merical set-up, such behaviour has proven to be hard totreat numerically when neutrino losses are important.This prevents a comprehensive treatment of impulsiveheating events. The question if (and how) results ob-tained in a 3D framework are different with respect tothe 2D treatment of Pons & Rea (2012) is a matter thatwill be addressed in a future study.5.1. Ramifications
Comparison with low-B magnetars — The presence ofstrong toroidal fields in magnetars has been invokedsince long to explain their distinguishing activity as com-pared to radio-pulsars with similar spin-down magneticfields, the high-B pulsars with B dip ≈ –10 G (seee.g. Turolla et al. 2015). On the other hand, somesources with B dip as low as ≈ G can show magnetar-like activity (the low-B magnetars; see e.g. Turolla et al.2011, and references therein). According to our sim-ulations, the resistive tearing instability appears on atimescale (cid:46) τ H ≈ yr and lasts for about ≈ (cid:46) yr), which are the vast ma-jority of the magnetar population. Whether such aninstability can be triggered under the conditions typicalof older objects, like the low-B sources SGR 0418+5729and Swift J1822.3-1606 (age ≈ –10 yr; Turolla et al.2011; Rea et al. 2012) or the onset of outbursts is pro-duced by a different, possibly related, mechanism is anopen question. Crescent-shaped features and observations — Non-polar,crescent-like hot spots have been recently detected inNICER X-ray observations of the millisecond pulsarPSR J0030+0451 and interpreted as due to heating frombackflowing currents in a non-dipolar magnetic field ∼ pulsar/magnetar/main.html (Olausen & Kaspi 2014). D magneto-thermal evolution of neutron stars
Battery effects and magnetar magnetospheres — In section4 we showed how the magnetic field created throughbattery effects by an external heating source can reachstrong local values in a turbulent-like pattern. The ex-istence of small-scale magnetic structures, in which thefield strength is orders of magnitude higher than in thesurrounding dipole, has been invoked to explain the (rel-atively) large energy ( ≈ E cp (cid:39) . B/ G)keV. The prototypical source is the low-field ( B dip ∼ × G) magnetar SGR 0418+5729, where a phase-dependent absorption feature at ∼ XMM-Newton data by Tiengo et al. (2013).According to their interpretation, the line arises as ra-diation from a cooling spot on the star surface crossesa baryon-loaded, small ( ≈
100 m) magnetic loop whereabsorption occurs. Albeit associating this kind of mag-netic structure with those produced by the battery effectin our simulations is tempting, we warn that thermocou-pling effects turn out to be less important in the case ofdeep heating (see section 4.1), where the local enhance-ment of the magnetic field is modest.5.2.
Present limitations
In this work, we highlighted the perspectives that anovel three dimensional approach can open in the studyof neutron star magneto-thermal evolution. There are,nevertheless, some limitations that must be taken intoaccount when interpreting our numerical results.In fact, we had to reduce the microphysical input toa realistic yet simplified model for the computing timeto be manageable. This concerns in particular the useof a simplified form for hydrostatic equilibrium densityprofile of equation (18), which was also assumed inde-pendent on temperature and magnetic field, and the useof a constant τ throughout the crust. Moreover, we havechosen some strong prescriptions on thermal conductiv-ity and heat capacity. In particular, the assumption that C V is linearly dependent on the temperature is valid only for the electron contribution, and does not takeinto account the contribution of the lattice. This impliesthat equation (6) depends on T only, which is a keypoint for the efficiency of the numerical scheme. How-ever, such an assumption becomes questionable whenthe term ∝ ∂T /∂t starts to dominate, as in the caseof impulsive heating in section 4.1. In particular, thisaffects our estimates for the duration of thermal relax-ation events. In fact, the heat diffusion timescale acrossa length L can be estimated as (Chaikin et al. 2018) τ diff ∼ (cid:34)(cid:90) L d l (cid:18) C V κ (cid:19) / (cid:35) , (30)hence it is regulated by the specific heat-to-thermal con-ductivity ratio. According to the estimates of Chaikinet al. (2018), for the typical conditions of a NS, thetimescale of heat transport from an internal heater tothe surface is (cid:46) τ diff ≈
50 yr. This may well berelated to our assumptions which make the ratio C V /κ is independent on temperature, while it is expected todepend from temperature as well as from the propertiesof crustal superfluidity (Potekhin et al. 2015, and refer-ences therein). Hence, our results for this case shouldsimply be regarded as indicative of the general evolu-tion of such events. For the model discussed above theevolution timescale is longer than what expected undermore realistic conditions by a factor ≈ De Grandis et al.
ACKNOWLEDGMENTSSimulations were run at CloudVeneto, a HPC facilityjointly owned by the University of Padova and INFN,and at UCL Grace HPC facility (Grace@UCL). Theauthors gratefully acknowledge the use of both facili-ties and the associated support services. RT and RTacknowledge financial support from the Italian MIURthrough grant PRIN 2017LJ39LM.REFERENCES
Abbott, B. P., et al. 2017, ApJL, 848, L13,doi: 10.3847/2041-8213/aa920cAguilera, D. N., Pons, J. A., & Miralles, J. A. 2008, A&A,486, 255, doi: 10.1051/0004-6361:20078786Aubert, J., Aurnou, J., & Wicht, J. 2008, GeophysicalJournal International, 172, 945,doi: 10.1111/j.1365-246X.2007.03693.xBaym, G., Pethick, C., Pines, D., & Ruderman, M. 1969,Nature, 224, 872, doi: 10.1038/224872a0Becker, W. 2009, Astrophysics and Space Science Library,Vol. 357, X-Ray Emission from Pulsars and NeutronStars, ed. W. Becker (Springer), 91,doi: 10.1007/978-3-540-76965-1 6Beloborodov, A. M. 2009, ApJ, 703, 1044,doi: 10.1088/0004-637X/703/1/1044Bilous, A. V., Watts, A. L., Harding, A. K., et al. 2019,ApJL, 887, L23, doi: 10.3847/2041-8213/ab53e7Chaikin, E. A., Kaminker, A. D., & Yakovlev, D. G. 2018,Ap&SS, 363, 209, doi: 10.1007/s10509-018-3393-zChandrasekhar, S., & Kendall, P. C. 1957, ApJ, 126, 457,doi: 10.1086/146413Chugunov, A. I., & Horowitz, C. J. 2010, MNRAS, 407,L54, doi: 10.1111/j.1745-3933.2010.00903.xCiolfi, R., & Rezzolla, L. 2013, MNRAS, 435, L43,doi: 10.1093/mnrasl/slt092Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., &Esposito, P. 2018, MNRAS, 474, 961,doi: 10.1093/mnras/stx2679Cumming, A., Arras, P., & Zweibel, E. 2004, ApJ, 609, 999,doi: 10.1086/421324Dormy, E., Cardin, P., & Jault, D. 1998, Earth andPlanetary Science Letters, 160, 15,doi: 10.1016/S0012-821X(98)00078-8Goldreich, P., & Reisenegger, A. 1992, ApJ, 395, 250,doi: 10.1086/171646Gourgouliatos, K. N., & Cumming, A. 2014a, MNRAS, 438,1618, doi: 10.1093/mnras/stt2300 —. 2014b, PhRvL, 112, 171101,doi: 10.1103/PhysRevLett.112.171101Gourgouliatos, K. N., & Pons, J. A. 2020, arXiv e-prints,arXiv:2001.03335. https://arxiv.org/abs/2001.03335Gourgouliatos, K. N., Wood, T. S., & Hollerbach, R. 2016,Proceedings of the National Academy of Science, 113,3944, doi: 10.1073/pnas.1522363113Greif, S. K., Raaijmakers, G., Hebeler, K., Schwenk, A., &Watts, A. L. 2019, MNRAS, 485, 5363,doi: 10.1093/mnras/stz654Gudmundsson, E. H., Pethick, C. J., & Epstein, R. I. 1983,ApJ, 272, 286, doi: 10.1086/161292Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007,Neutron Stars 1 : Equation of State and Structure, Vol.326Haskell, B. 2008, Classical and Quantum Gravity, 25,114049, doi: 10.1088/0264-9381/25/11/114049Ho, W. C. G., Andersson, N., & Graber, V. 2017, PhysicalReview C, 96, doi: 10.1103/physrevc.96.065801Hollerbach, R., & R¨udiger, G. 2004, MNRAS, 347, 1273,doi: 10.1111/j.1365-2966.2004.07307.xHorowitz, C. J., Berry, D. K., & Brown, E. F. 2007,PhRvE, 75, 066101, doi: 10.1103/PhysRevE.75.066101Kantor, E. M., & Gusakov, M. E. 2007, MNRAS, 381,1702, doi: 10.1111/j.1365-2966.2007.12342.xKarageorgopoulos, V., Gourgouliatos, K. N., &Contopoulos, I. 2019, MNRAS, 487, 3333,doi: 10.1093/mnras/stz1507Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55,261, doi: 10.1146/annurev-astro-081915-023329Kondratyev, I., Moiseenko, S., Bisnovatyi-Kogan, G., &Glushikhina, M. 2019, in High Energy Phenomena inRelativistic Outflows VII, 59.https://arxiv.org/abs/1912.05980Lander, S. K. 2014, MNRAS, 437, 424,doi: 10.1093/mnras/stt1894
D magneto-thermal evolution of neutron stars15