{'fraction_non_alphanumeric': 0.03134986083908876, 'fraction_numerical': 0.018864034635604578, 'mean_word_length': 4.6344562218672865, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unifying the observational diversity of isolated neutron stars via magneto-thermal evolution models

Observations of magnetars and some of the high magnetic field pulsars have shown that their thermal luminosity is systematically higher than that of classical radio-pulsars, thus confirming the idea that magnetic fields are involved in their X-ray emission. Here we present the results of 2D simulations of the fully-coupled evolution of temperature and magnetic field in neutron stars, including the state-of-the-art kinetic coefficients and, for the first time, the important effect of the Hall term. After gathering and thoroughly re-analysing in a consistent way all the best available data on isolated, thermally emitting neutron stars, we compare our theoretical models to a data sample of 40 sources. We find that our evolutionary models can explain the phenomenological diversity of magnetars, high-B radio-pulsars, and isolated nearby neutron stars by only varying their initial magnetic field, mass and envelope composition. Nearly all sources appear to follow the expectations of the standard theoretical models. Finally, we discuss the expected outburst rates and the evolutionary links between different classes. Our results constitute a major step towards the grand unification of the isolated neutron star zoo.


INTRODUCTION
The great advances in X-ray observations during the last several decades have improved our understanding of the physics of neutron stars (NSs). Isolated NSs, as seen in X-rays, are phenomenologically quite heterogeneous. The high X-ray luminosities, long periods, and the occurrence of bursts and outbursts characteristic of Anomalous Xray Pulsars (AXPs) and Soft Gamma Repeaters (SGRs) (Mereghetti 2008), are interpreted as different facets of the restless dynamics of a strong magnetic field in these (typically) young sources, in agreement with the magnetar model (Thompson & Duncan 1995). Many radio-loud Rotation-Powered Pulsars (RPPs) are also detected in X-rays (Becker 2009). Their emission is usually powered by the conversion of a fraction of their rotational energy into non-thermal radiation; however, in a few cases, they also show thermal emission. In particular, high-B pulsars  appear to connect the X-ray-quiet standard radio-pulsars with the very active magnetars, showing intermediate luminosities, and sporadic magnetar-like activity. The radioquiet X-ray Isolated NSs (XINSs) are a class of relatively old, nearby cooling NSs, with the cleanest detected thermal emission and a relatively large magnetic field. Last, Central Compact Objects (CCOs, Gotthelf et al. 2013) represent a handful of puzzling, radio-quiet sources which in some cases combine a very weak external magnetic field with a relatively large luminosity and evidence for anisotropic surface temperature distribution.
Although isolated NSs have been divided into all these observational classes for historical reasons, a unifying vision considers them as different manifestations of the same underlying physics (Kaspi 2010). In this context, one of the main theoretical tasks is to explain the varied phenomenology of their X-ray emission. X-ray spectra carry precious information about the surface temperature and the physics driving the cooling of the NSs. The detected X-ray flux, if accompanied by a reliable distance measurement, leads to an estimate of the bolometric luminosity, which can be confronted with cooling models to infer properties of dense matter in the NS interior. In addition, timing properties (period and period derivative) provide us with information about the rotational energy loss, which is believed to be regulated mainly by the dipolar component of the external magnetic field.
The theory of NS cooling has been developed since short after the first hints of X-ray emission from the surface of NSs in the Sixties (e.g., Morton 1964;Chiu & Salpeter 1964;Tsuruta & Cameron 1965). During the following decades, microphysical inputs have been refined to build state-ofthe-art 1D cooling models (see Yakovlev & Urpin 1981;Nomoto & Tsuruta 1986;Page & Baron 1990;van Riper 1991;Page & Applegate 1992;Pethick 1992, and references within for pioneering results), but the effects of magnetic fields were barely taken into account (Miralles et al. 1998;Page et al. 2000). Page et al. (2004Page et al. ( , 2009 proposed the socalled minimal cooling scenario, in which all the necessary known microphysical ingredients are consistently included (specific heat, conductivities, effects of superfluidity, neutrino emission processes), but fast neutrino-cooling processes (direct URCA by nuclear or exotic matter) are excluded. The minimal cooling paradigm has been shown to be consistent with the luminosity of the weakly magnetised NSs ( 10 13 G), but magnetars and some of the high-B pulsars show a relatively large X-ray luminosity, thus requiring the presence of additional heat sources, likely linked to magnetic field evolution.
The evolution of the magnetic field in NSs has been studied extensively by a number of authors (Goldreich & Reisenegger 1992;Geppert & Rheinhardt 2002;Hollerbach & Rüdiger 2002, 2004Cumming et al. 2004;Arras et al. 2004;Gonzalez & Reisenegger 2010;Glampedakis et al. 2011). In the solid crust of NSs, the magnetic field, besides being dissipated by the Joule effect, evolves due to the Lorentz force acting on the electron fluid, the so-called Hall term. When this term dominates, electric currents are squeezed in a smaller volume and move through the crust, creating small scale structures and allowing an interplay between poloidal and toroidal components of the magnetic field. The Hall term plays an important role in the evolution of large magnetic fields, B 10 14 G, and the complexity of Hall dynamics has limited previous studies. In Aguilera et al. (2008), a 2D cooling code followed the evolution of the temperature, for a given analytical prescription of a simplified, homologous decay of the magnetic field. In , a spectral code follows the evolution of the magnetic field, including both the resistive and the Hall terms, but the temperature evolution was described by an isotropic, analytical cooling law. However, to describe the evolution of strongly magnetised NSs, the coupling between temperature and magnetic field must be taken into account. The first 2D simulations of the fully coupled magneto-thermal evolution ) included only Ohmic dissipation and, only recently, we have overcome numerical difficulties in the treatment of the important Hall term for large magnetic fields. In  we presented the first finite-difference, 2D magneto-thermal code able to manage arbitrarily large B intensities while consistently treating the Hall term. This is crucial to test whether dissipation of currents in the NS crust in realistic models can explain the observational data.
In this paper, we present results for the magnetic field and temperature evolution of NSs by means of a state-of-theart code for the computation of the 2D magneto-thermal evolution. We couple the code described in Aguilera et al. (2008) for the temperature evolution with the relativistic version of the recently developed code of magnetic field evolution presented by . This allows us to study in detail the expected properties of objects of any age. We apply the code to a systematic study of the magnetothermal evolution of NSs, revising the results from previous studies (Aguilera et al. 2008;Pons et al. 2009) and comparing them with the observations of different classes of isolated NSs. For a meaningful comparison of our models with the best available data, we have re-analized the most accurate observations in the XMM-Newton and Chandra archives for all isolated neutron stars (among radio pulsars, high-B pulsars, magnetars, XINSs, and CCOs) with a detected (noncontroversial) surface thermal emission.
In §2, we summarize the spectral data of the isolated NSs that we have analysed, presenting an updated table for their timing and spectral properties. In §3, we recall the basic equations and the state-of-the-art microphysics setup. In §4, we present the results of the simulations, focusing on the evolution of magnetic field, surface temperature, and luminosity. In §5, we compare timing and spectral properties of a sample of magnetised sources, in the spirit of providing a unified vision of isolated NSs. Conclusions are presented in §6.

DATA ON COOLING NEUTRON STARS
To confront theoretical cooling models with observational data, we need to know simultaneously the age and some quantity related to the thermal emission from the NS surface (luminosity or temperature). The number of sources for which both measures are available is limited, and in most cases subject to large uncertainties. In this section, we will discuss in detail the sample of selected sources (see Table 3 for the complete list), which includes: • Four CCOs (Gotthelf et al. 2013), including the very young NS in Cassiopeia A, and the only three CCOs with a measured value of P andṖ . The discrepancy between the CCOs characteristic age and the estimated age of their SNRs, together with their very slow spin-down, indicates that magnetic braking is much less effective than in other classes, and their spin period is very close to the natal one. We have ignored the other CCOs candidates since they have spectral information with poor statistics and/or a very uncertain age of the associated SNR.
• Eight rotation powered pulsars, including the Vela pulsar and the so-called three Musketeers (PSR B0656, PSR B1055 and the γ-ray-loud, radio-quiet Geminga; De Luca et al. 2005). We have excluded most of the young pulsars, many of which associated with pulsar wind nebulae (see Becker 2009 for an observational review of X-ray pulsars), since in those cases data are compatible with non-thermal emission powered by the rotational energy loss, which is orders of magnitude larger than their X-ray luminosity (i.e. Crab pulsar and RX J0007.0+7303 in SNR CTA1; Caraveo et al. 2010). We also exclude several old pulsars (Zavlin & Pavlov 2004)  of m 2 ), since the temperature of the small hot spots is probably unrelated to the cooling history of the NS.
• The seven X-ray nearby Isolated NSs (XINSs) known as the magnificent seven (Haberl 2007;Kaplan & van Kerkwijk 2009). All of them have good spectra, and in most cases well determined timing properties and good distance determinations (often with direct parallax measurements).
• Four high-B radio-pulsars (HB, see Ng & Kaspi 2011 for a review), with inferred magnetic fields Bp ∼ 10 13 − 10 14 G and good quality spectra. We have included the only RRAT detected so far in X-ray (PSR J1819, McLaughlin et al. 2007). We have excluded the magnetar-like pulsar PSR J1846-0258 since during quiescence its X-ray emission does not show a significant thermal component (Ng et al. 2008;Livingstone et al. 2011), and it is orders of magnitude smaller than its rotational energy loss. We now discuss separately the timing and spectral properties of our sample. All the data presented in the following subsections, with links to abundant references, can also be found in our website 1 . We plan to periodically update and extend it.

On timing properties and age estimates.
If both the spin period (P ) and the period derivative (Ṗ ) of the source are known, the characteristic age τc = P/2Ṗ can be used as an approximation to the real age, with which it coincides only if the initial period was much shorter than the current value and the magnetic field has been constant during the entire NS life. Unfortunately this is not the most common situation, and usually, for middle-aged and old objects, τc is found to be larger than the real age, when the latter has been obtained by kinematic measurements (hereafter, kinematic age τ k ). When the object is located in a supernova remnant, the kinematic age can be inferred by studying the expansion of the nebula (see Allen & Horvath 2004 for a review with particular attention to the magnetar associations). For a few other nearby sources (e.g., some XINSs and few magnetars; Tetzlaff et al. 2011; Tendulkar et al. 2012), the proper motion, with an association to a birth place, can give the kinematic age. We have collected the most updated and/or reliable available information on timing properties and kinematic age from the literature, the ATNF catalog 2 (Manchester et al. 2005), and the McGill online magnetar catalog 3 . We present in Table 1 all the sources in our sample, with their well established associations (Gaensler et al. 2001), the known timing properties with the characteristic age and the alternative estimate for the age, when available. The inferred value of the surface, dipolar magnetic field at the pole Bp, assuming the standard dipole-braking formula, is shown as well. We note that, for RPPs timing properties are stable over a time span of tens of years, and P andṖ are precisely measured, but for magnetars the timing noise is much larger. In some cases, different values ofṖ have been reported, differing even by one order of magnitude (see for instance Table 2 of the online McGill catalog). Consequently, we should take these values with caution, especially for the most extreme objects (largest values ofṖ , see § 5 for further discussion).

On luminosities and temperatures from
spectral data analysis.
Luminosities and temperatures can be obtained by spectral analysis, but it is usually difficult to determine them accurately. The luminosity is always subject to the uncertainty in the distance measurement, while the inferred effective temperature depends on the choice of the emission model 2 http://www.atnf.csiro.au/people/pulsar/psrcat/ 3 http://www.physics.mcgill.ca/∼pulsar/magnetar/main.html (blackbody vs. atmosphere models, composition, condensed surface, etc.), and it carries large theoretical uncertainties in the case of strong magnetic fields. We often find that more than one model can fit equally well the data, without any clear, physically motivated preference for one of them, with inferred effective temperatures differing up to a factor of two. Photoelectric absorption from interstellar medium further constitutes a source of error in temperature measurements, since the value of the hydrogen column density NH is correlated to the temperature value obtained in spectral fits. Different choices for the absorption model and the metal abundances can also yield different results for the temperature. In addition, in the very common case of the presence of inhomogeneous surface temperature distributions, only an approximation with two or three regions at different temperatures is usually employed. Last, in the case of data with few photons and/or strong absorption features, the temperature is poorly constrained by the fit, adding a large statistical error to the systematic one. For all of these reasons, the temperatures inferred by spectral fits can hardly be directly compared to the physical surface temperatures extracted from cooling codes.
Because of the above considerations, the luminosity constitutes a better choice to compare data and theoretical models. Since it is an integrated quantity, it averages effects of anisotropy and the choice of spectral model. The main uncertainty on the luminosity is often due to the poor knowledge of the source distance. In the worst cases, the distance is known within an error of a few, resulting in up to one order of magnitude of uncertainty in the luminosity. In addition, the interstellar absorption acts predominantly in the energy band in which most of the middle age NSs emit (E 1 keV). For this reason, hottest (magnetars) or • The source has been recently discovered in outburst and it could have not reached the quiescence yet at the time of the latest available observation. ⋆ Absorption line(s) gabs included in the fit. † Hard tail detected in quiescence hard X-ray spectra (E 20 keV). Table 3. Emission properties of the thermally emitting neutron stars. f X is the 1-10 keV band unabsorbed flux according to the indicated best fit model (the bbodyrad model flux is chosen whenever multiple models are equally compatible with the data). kT bb and R bb are the temperature and radius inferred by the bbodyrad model. kTnsa is the temperature inferred by the nsa model for the weakly magnetised cases (B 5 × 10 13 G) with acceptable associated radius Rnsa, also indicated. kTrcs is the temperature inferred by the RCS model for strongly magnetised sources. L is the bolometric luminosity from the thermal component(s) of the fit, and assuming the indicated distance (whose references are listed in the online table). The range of L includes both statistical and distance errors; for strongly absorbed sources (i.e., most of magnetars) a minimum arbitrary factor of 50% uncertainty is assumed to account for systematical model-dependent uncertainty. kT cool is either the lower temperature for models including 2 BB, compatible with emission from the entire surface, or the upper limit for cases showing emission from a small spot R bb ∼ few km. In the latter case, L cool is the associated upper limit to the "hidden" thermal luminosity. See text for details on the spectral models. All radii, temperatures and luminosities are the values as measured by a distant observer. closest (XINSs) sources are easier to detect. Similarly to the case of the temperature, the choice of different models of absorption and chemical abundances can yield additional systematic errors on the luminosity. However, for the worst cases, the relative error is about 30%, making it usually a secondary source of error compared with the distance.

Data reduction
In order to control the systematic errors in the luminosity and spectral parameters, we have homogeneously reanalysed all the data, extracted directly from the best available observations from Chandra or XMM-Newton. In Table 2 we list the log of all the observations we used. These observations have been selected with the following criteria.
(i) For comparable exposure times of a given source, we always preferred the XMM-Newton observation, given the larger collecting area with respect to Chandra, hence resulting in a more detailed spectrum. Instead we used Chandra data for sources which have bright nebulae.
(ii) For variable sources we took the longest available observation during the quiescent state of the NS.
(iii) We excluded those objects which were fitted equally well without the addition of a thermal component, hence the latter not being statistically significant in the data.
We processed all XMM-Newton observations (Jansen et al. 2001) listed in Table 2 using SAS version 11, and employing the most updated calibration files available at the time the reduction was performed (November 2012). Standard data screening criteria are applied in the extraction of scientific products, and we used only the EPIC-pn camera data. We extracted the source photons from a circular region of 30 ′′ radius, and a region of the same size was chosen for the background in the same CCD but as far as possible from the source position. 4 We restricted our analysis to photons having PATTERN 4 and FLAG=0.
The Chandra data we used in this work were all taken with the Advanced CCD Imaging Spectrometer (ACIS-S; Garmire et al. (2003). Data were analyzed using standard cleaning procedures 5 and CIAO version 4.4. Photons were extracted from a circular region with a radius of 3 ′′ around the source position, including more than 90% of the source photons, and background was extracted from a region of the same size, far from the source position.

Data analysis
All spectra were grouped in order to optimize the signal to noise in each spectral bin, but even for the fainter objects we grouped the spectra in order to have at least 25 counts per spectral bin, and reliably use the chi-squared statistic to assess the goodness of the fits. The response matrices were built using ad-hoc bad-pixel files built for each observation.
For the spectral analysis we used the XSPEC package (version 12.4) for all fittings, and the phabs photoelectric absorption model with the Anders & Grevesse (1989) abundances, and the Balucinska-Church & McCammon (1992) photoelectric cross-sections. We usually restricted our spectral modeling to the 0.3-10 keV energy band, unless the source was such that a smaller energy range was needed for the spectral analysis. We always excluded bad channels when needed.
Depending on the source, we either used a blackbody model (bbodyrad) alone, or added a power law (power) and/or a second blackbody component if statistically significant. When physically motivated, we also fitted data with an NS atmosphere model (nsa: Zavlin et al. 1996;Pavlov et al. 1995), always fixing the NS mass to M = 1.4 M⊙, the radius to R = 10 km and B to the closest value to the inferred surface value Bp, see Table 3.
We find that for CCOs the atmosphere models provide slightly better fits and radii closer to the typical NS values. In the case of RPPs, they can all be described by thermal component(s), with an additional power law. In particular, the three Musketeers (Geminga, PSR B0656 and PSR B1055) are well fit by a double blackbody plus a power law. An atmosphere model (plus a power law) fits Vela better than a blackbody plus power law, and it is compatible with emission from the entire surface. However, in several other cases atmosphere models provide unphysical emitting radii. On the other hand, XINSs X-ray spectra are well fit by a single blackbody model, with absorption features in a few cases; in a couple of cases an atmosphere model can fit the data as well.
For magnetars we have also used a resonant Compton scattering model (RCS: Rea et al. 2008;Lyutikov & Gavriil 2006), adding a power law component when needed. The lack of a unifying description of the magnetosphere, consistently dealing with both magnetic field configuration and Building on previous works (Page et al. 2004;Aguilera et al. 2008;Page et al. 2009;Pons et al. 2009) we have updated and extended our magneto-thermal evolution code in two major ways: the proper treatment of the important Hall term in the induction equation describing the magnetic field evolution , and updated microphysics inputs. This allows us to follow the long-term evolution of magnetized NSs, a necessary step to understand timing and spectral properties of isolated NSs at different ages. In this section, we briefly summarize the equations, the method, and the updated ingredients of the simulations.

Basic equations.
For our purposes, the small structural deformations induced by rotation and the magnetic field can be safely neglected. To include general relativistic effects, we consider the standard static metric where e ν (r) is the lapse function that accounts for redshift corrections, e λ(r) = (1 − 2Gm(r)/c 2 r) −1/2 is the space curvature factor, and m(r) is the enclosed gravitational mass within radius r.
The star thermal evolution is described by the following energy balance equation where cv is the volumetric heat capacity,κ is the thermal conductivity tensor, Qν are the energy loss by neutrino emission per unit volume, and Q h is the Joule heating rate per unit volume. The conduction of heat becomes anisotropic under the presence of a strong magnetic field, which also controls the heating rate and, secondarily, affects the rate of a few neutrino processes (Aguilera et al. 2008). This is the reason why, as shown in previous studies, the thermal evolution is coupled to the magnetic field evolution. In the crust, ions form a Coulomb lattice, while electrons are relativistic, degenerate and can almost freely flow, providing the currents that sustain the magnetic field. The evolution of the system is governed by the Hall induction equation which, using the same notation as in Pons et al. (2009), has the form: where the conductivity σ takes into account all the (strongly temperature-dependent) electron processes. The first term on the right hand side accounts for Ohmic dissipation, while the second term is the Hall effect. The magnetization parameter ωBτe ≡ σB cene (where ωB = eB/m * e c is the gyration frequency of electrons, with τe and m * e are the relaxation time and effective mass of electrons), affects both the evolution of magnetic field and temperature. In the regime where ωBτe ≫ 1 (strong magnetic fields, 10 14 G, and temperatures 5 × 10 8 K, see Aguilera et al. 2008;Pons et al. 2009; for more details), the conductivity tensorκ becomes very anisotropic, the Hall term dominates, and the induction equation acquires a hyperbolic character. The Ohmic and Hall timescales vary by orders of magnitude within the crust and during the evolution, depending strongly on density, temperature, and magnetic field intensity and curvature.
To overcome the limitations of previous hybrid methods (spectral in angles, finite-differences in the radial direction) used in  and Pons et al. (2009), we use a new finite-difference code described in . The main effect of the Hall term is to transfer part of the magnetic energy from large to small scales, and between poloidal and toroidal components. In the case of strong toroidal components, it also leads to the formation of discontinuities of the tangential components of the magnetic field, i.e. current sheets, where the dissipation is strongly enhanced . This directly affects the thermal evolution through the term Q h in Eq. (2).
The core of NSs (or at least a fraction of its volume) is thought to be a type II superconductor. The dynamics of the magnetic field in the core is not clearly understood.
Here, Ohmic dissipation acts over long timescales ( Gyr) due to the very high conductivity of nuclear matter, but other mechanisms such as the interplay between flux tubes and vortices, magnetic buoyancy, or ambipolar diffusion may operate to expel magnetic flux from the core on timescales shorter or comparable to the cooling timescale ( 10 5 ) yr. The detailed study of these mechanisms goes beyond the purpose of this paper, in which we mostly consider models with the magnetic field confined into the crust or, for the models with magnetic fields permeating the core, we only include Ohmic dissipation. This implies that the core magnetic field is basically frozen during cooling timescales.

Microphysics.
The microphysical inputs that enter in the magneto-thermal evolution equations, Eqs. (2, 3), are the thermal conductivity, the specific heat, the neutrino emissivities, and the electrical conductivity. For a complete review of the physics of the crust, with the thermal and transport properties, see Chamel & Haensel (2008); Page & Reddy (2012).
Unless it is explicitly mentioned, we use the same inputs as in Aguilera et al. (2008) and Pons et al. (2009), and references therein. In our current version of the code, we have updated some of these ingredients as follows. For the specific heat, we use the latest version of the equation of state for a strongly magnetized, fully ionized electron-ion plasma (Potekhin & Chabrier 2010) 6 . For the thermal and electrical conductivities in the crust we use the latest FORTRAN routines which include quantizing effects and thermal averaging corrections accounting for partial degeneration of electrons, and the effect of the finite size of nuclei in the inner crust 7 . We have adopted the functional form of the superfluid gaps proposed in Table 1 of Ho et al. (2012). In addition, the contribution to the thermal conductivity from superfluid phonons (Aguilera et al. 2009) has been suppressed, according to recent results which indicate that entrainment between super-  Tables 1 and 3. The left panel corresponds to models with iron envelopes and the right panel to models with light-element envelopes. Arrows label sources for which τc 50 kyr, and no kinematic age is available, so that the real age is expected to be shorter. An uncertainty of 50% has been arbitrarily taken for the kinematic age when error estimates have not been found in the literature.
fluid neutrons and nuclei is larger than expected, resulting in an increased effective mass of the nuclei in the inner crust (Chamel 2012;Chamel et al. 2012).
We note that, even for 10 15 G, including quantizing effects in the crust does not substantially affect the long-term evolution of the star. However, it may have a strong effect on the envelope models, as we discuss later. It should also be noted that the electrical conductivity of the NS crust at low temperatures is dominated by the interaction with the atomic nuclei of the crystal lattice. The phonon Umklapp process is modeled by an analytical formula including the exponential reduction of the scattering rate (Gnedin et al. 2001). Chugunov (2012) shows that, at low temperatures (T < 10 7 K), such formula overestimates the electrical conductivity by several orders of magnitude. However, here we are interested in the evolution of a NS during the first ∼ Myr of its life, when its temperature is still relatively high and the used approximation is still valid.
Electron-impurity collisions become the dominant process at low temperatures T 10 8 K (Shternin & Yakovlev 2006 and refs. within). Therefore, the conductivity at late times is almost independent of the temperature and is determined by the impurity content of the crust. This is usually parametrised by the quadratic deviation of the atomic number Qimp = i xi(Z 2 i − Z 2 ), which is highly uncertain. In isolated NS, especially in magnetars which remain warmer for longer times, Qimp in the outer crust is expected to be low, unlike in accreting NSs in binary systems, where the outer crust is being continuously replenished by newly processed nuclei.
For the inner crust, things are different. First, the crust is for the most part an elastic solid, comprising a Coulomb lattice of normal spherical nuclei. However, in the innermost layers near the crust/core boundary, because of the large effect of the Coulomb lattice energy, cylindrical and planar geometries can occur, both as nuclei and as bubbles (Ravenhall et al. 1983). These phases are collectively named nuclear pasta (by analogy to the shape of spaghetti, maccaroni and lasagna). More sophisticated molecular dynamics simulations (Horowitz et al. 2005;Horowitz & Berry 2008) have shown that it may be unrealistic to predict the exact sizes and shapes of the pasta clusters, and that the actual shape of the pasta phase can be very amorphous, with a very irregular distribution of charge. This is expected to have implications on the transport properties, in particular on the electrical resistivity. Other simulations indicate that the inner crust could be amorphous and heterogeneous in nuclear charge Jones (2004). The expected values of Qimp, coming from the crust formation process alone, could be of the order O(10), providing large values of resistivity.
To include the effect of pasta phases or disorder, we will extend the impurity parameter formalism to the inner crust in a simplified way. We set Qimp = 0.1, except in the pasta region (the crust where ρ > 6 × 10 13 g/cm 3 ), where we will let the parameter vary between Qimp = 1 − 100. Recently we have discussed in detail the influence of the pasta phase or disorder on the expected evolution of timing properties ).

The envelope.
The envelope (so-called the ocean) is the liquid layer above the solid crust and below the atmosphere. It comprises approximately the outermost hundred meters of the star, and it is the region with the largest temperature gradient. Due to its relatively low density, its thermal relaxation timescale is much shorter than that of the crust, which makes any attempt to perform cooling simulations in a numerical grid that includes both regions simultaneously computationally too expensive. The widely used approach is to use fits to stationary envelope models to obtain a phenomenological relation between the temperature at the bottom of the envelope, T b , and the effective temperature Ts. This T b − Ts relation depends on the surface gravity gs, the envelope composition, and the magnetic field strength and orientation.
After the pioneer work by Gudmundsson et al. (1983), further 1D, plane-parallel simulations have followed the same approach, including neutrino processes and the presence of magnetic field, which affects the transport of heat, causing the heat flux to become anisotropic Potekhin et al. 2007). Therefore, the geometry of the magnetic field is important, since regions permeated by radial magnetic field lines are thermally connected to the interior, while zones with tangential magnetic field (equator in the dipolar case) are thermally insulated. Pons et al. (2009) revised the envelope model with a 2D code, allowing for meridional transport of heat, finding that the most important effect is the attenuation of the anisotropy. For consistency, in this paper we have built our own envelope models, revising the fits done in section 3 of Pons et al. (2009) with the updated microphysical inputs. We have placed the outer boundary of our crust at ρ = 3 × 10 10 g cm −3 , and solved different stationary envelope models to obtain the T b − Ts relations that are then used as an external boundary condition in our simulations.
The luminosity seen by an observer at infinity is obtained by assuming blackbody emission at the temperature Ts from each patch of the NS surface. It is not the purpose of this paper to study other emission models, such as as magnetic atmospheres or condensed surface models. These do not change appreciably the relation between T b and the emitted flux (Potekhin et al. 2007) but can be important to interpret the temperature inferred by a spectral fit (see § 2).

Cooling of weakly magnetized neutron stars.
Before we discuss the effects of strong magnetic fields on the cooling history of NSs, we begin by briefly revisiting the cooling of low magnetic field sources, with our microphysical inputs. Comprehensive reviews can be found in Page et al. (2004) and Yakovlev & Pethick (2004).
In Fig. 1 we present the set of cooling curves for nonmagnetised NSs with masses ranging between 1.10 and 1.76 M⊙ (lines from top to bottom), together with the available observational data presented in §2. We show the bolometric luminosity as a function of time. The left/right panels correspond to models with heavy element envelopes (iron) and light-element (hydrogen) envelopes, respectively. Sources with estimates of the kinematic age are shown with the associated error bar on the age. With arrows pointing towards the left we indicate the sources with τc 50 kyr and no kinematic age. This implies a likely overestimate of the real age for the middle-age and old objects; such an age overstimate could be even orders of magnitude for the oldest objects (like SGR 0418+5729).
In general, after ≈ 100 yr, low mass stars (M 1.4M⊙) are brighter than high mass stars. For the high-mass family, M 1.4M⊙, the proton fraction in the center of the star (for our choice of equation of state) is high enough to activate the direct URCA processes, which results in fast cooling before one hundred years. Within the low-mass family, cooling curves are similar at early ages (< 100 yr). The differences at t ∼ 10 2 −10 3 yr are due to the delayed transition of neutrons in the core to a superfluid state, which activates neutrino emission by means of Cooper pair formation and breaking. This effect also depends on the equation of state employed: a stiffer equation of state results in lower densities and the transition to superfluidity is further delayed, assuming the same gap. After the effect of the transition to a superfluid core is finished, at 10 3 yr, cooling curves for low mass NSs tend to converge again, following the same curve independently of the mass. For models with heavy element envelopes (left panel), the luminosity drops below 10 33 erg/s at most a few thousand years after birth, while light-element envelope models (right panel) predict luminosities around 10 33 − 10 34 erg/s for quite a long period (up to several 10 4 yr).
In young NSs, the evolution is dominated by neutrino cooling (the softest decay of the curves), while at late times, the emission of thermal photons from the star surface is the dominant factor (the steep decline of luminosity in the log-log plot). The turning point between the neutrinodominated cooling era and the photon-dominated cooling era happens much earlier for light-element envelope models (a few times 10 4 yr) than for heavy envelope ones (≈ 10 5 yr). During the photon cooling, the evolution of the surface temperature of the star can be roughly described by a power law Ts ∝ t −1/8α , with α ≪ 1, the exact value depending on the envelope model, which relates Ts ∝ T 0.5+α b (Page et al. 2004). For all weakly magnetised models, after ∼ 1 Myr, the luminosity has dropped below 10 31 erg/s, the surface temperature goes below 20 − 30 eV and the star becomes invisible to X-ray observations. Note also that in the photon cooling era NSs with light-element envelopes are much cooler Figure 3. Snapshots of the evolution of model A14 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, whose thickness has been enlarged by a factor of 4 for visualization purposes. Black lines are the poloidal magnetic field lines, while color scale indicates the toroidal magnetic field intensity (yellow: positive, red: negative). Bottom panels: temperature map inside the star. than those with iron envelopes. Page et al. (2011) suggested that a star could in principle have a light-element envelope at the beginning of its life, residual of accretion of fall-back material after the supernova explosion, that is progressively converted to heavier elements by nuclear reactions. From our quick glance at the classical cooling theory of low magnetic field NSs, the most interesting observation is that the high-B objects are systematically hotter than what theoretical non-magnetised cooling curves predict. This provides strong evidence in favor of the scenario in which magnetic field decay powers their larger luminosity. Even considering the (more than likely) overestimate of the NS age by using its characteristic age, which can reconcile some of the objects with standard cooling curves, it is clear that most magnetars and some high-B pulsars stars cannot be explained. Thus it is necessary to add the effects of the magnetic field in the simulations, extending the four parameter family of classical theoretical cooling models (compositional differences in the core, superfluid properties, composition of the envelope, and mass) with one more, the magnetic field. This is the subject of the next sections.

Initial magnetic field.
The initial configuration of the magnetic field in a newly born NS is still unknown. MHD equilibria considering a hot and fluid, newly born NS, point to a dipolar configuration with an equatorial torus where the lines are twisted (Ciolfi et al. 2009;Lander 2013, for instance). However, it is not clear whether the complicated dynamics of the proto-NS actually leads to this configuration, or if it is stable. We denote by "initial" the configuration at the epoch of the freezing of the crust (from hours to weeks after the SN explosion, depending on the density).
Lacking robust arguments in favor of specific initial configurations, we consider three generic families with different geometries, illustrated in Fig. 2 and mathematically described in § 2.1 of Aguilera et al. (2008). In the crustal confined model (type A geometry, left panel), the magnetic field lines do not penetrate into the core (Br = 0 at the crust/core interface), and the currents are entirely contained in the crust. The toroidal magnetic field is extended through the crust, and its form is chosen to be dipolar: where k is the normalization. In the type B configurations (central panel), the magnetic field threads the core, where the bulk of the current circulates with very slow dissipation. Type C geometry (right panel) is a hybrid case where a double system of currents supports a large-scale, core-extended dipole and an additional, stronger crustal magnetic field. Each family is parametrised by the intensity of the dipolar component at the pole, B 0 p , and the maximum intensity of the toroidal magnetic field at any point, B 0 t . In Table 4.1 we list the models studied in this work.
During the whole evolution we match the crustal magnetic field with a potential magnetic field at the surface, allowing for a general, time-dependent combination of multipoles (see  for details about the outer boundary condition). In particular, we follow the evolution of the external dipolar component Bp, which determines the timing properties (spin-down, braking index) of the star.

Evolution of a strongly magnetized NS.
We begin the discussion of our results by describing the general evolution of model A14. We run simulations up to several Myr, monitoring the temperature and magnetic field inside the star and at its surface. In the top panels of Fig. 3 we show three snapshots at t = 10 3 , 10 4 , 10 5 yr. For each snapshot, we show in the left hemisphere the surface temperature distribution, and in the right hemisphere the evolution of the crustal magnetic field. We now concentrate on the latter.
The first effect of the Hall term in the induction equation is to couple the poloidal and toroidal magnetic field, so that even if the latter is zero at the beginning, it is soon created. After ∼ 10 3 yr, the poloidal dipolar magnetic field has generated a mainly quadrupolar toroidal magnetic field, with a maximum strength of the same order of the poloidal component, with being B φ negative in the northern hemisphere and positive in the southern hemisphere. An important effect of the Hall term is that currents are gradually drained towards the inner crust, where they can last more. After ∼ 10 5 yr, the toroidal magnetic field and the distorted poloidal lines are concentrated in the inner crust. In the outer crust, the magnetic field approaches a potential configuration that matches the external magnetic field.
In the bottom panels of Fig. 3 we show three snapshots of the crustal temperature evolution. At t = 10 3 yr, the equator is hotter than the poles by a factor of three. As the evolution proceeds and currents are dissipated, the temperature reflects the change of geometry of the poloidal lines (see Fig. 3) and the anisotropy becomes weaker.
The presence of strong tangential components (B θ and B φ ) insulates the surface from the interior. For a dipolar geometry, the magnetic field is nearly radial at the poles, and these are thermally connected with the interior, while the equatorial region is insulated by tangential magnetic field lines. This has a twofold effect: if the core is hotter than the crust, the polar regions will be hotter than the equatorial regions; however, if Ohmic dissipation heats up the equatorial regions, they remain hotter than the poles. In addition, the insulating effect of the envelope must be taken into account. Depending on the local conditions (temperature, magnetic field strength, angle of the magnetic field with the normal to the surface) it may invert the anisotropy: in the surface temperature (see Fig. 3), the equator is cooler than the poles.
In Fig. 4 we show the evolution of the total magnetic energy, compared with the fraction stored in the different components. We also compare our results to the purely resistive case, in which the Hall term has been switched off (and therefore there is no creation of toroidal magnetic field). When the Hall term is included, ∼99% of the initial magnetic energy is dissipated in the first ∼ 10 6 yr, compared to only the ∼ 60% of 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 component. The poloidal magnetic field is supported by toroidal currents concentrated in the equator, which is hotter (see right bottom panel of Fig.  3), while the toroidal field is supported by poloidal currents that circulate in a higher latitude region, where the temperature is lower. Since the resistivity in this regime is strongly dependent on temperature, this causes the more rapid dissipation of the poloidal magnetic field. As a result, at late times most of the magnetic energy is stored in the toroidal component.

Dependence on mass and relevant microphysical parameters.
It is worth to briefly explore the dependence on the NS mass and some of the main microphysical parameters, such as the impurity content of the innermost part of the crust (pasta region). The top panel of Fig. 5 illustrates the sensitivity of the time evolution of Bp to microphysical parameters: we use the same initial magnetic field (model A14), but with different NS masses, and different value of Qimp (impurity parameter in the regions of the crust above ρ > 6 × 10 13 g cm −3 ). Comparing the evolution of the dipolar component of the magnetic field (top panel) for different NS masses and the same Qimp in the pasta region, we do not see significant differences. The small quantitative differences are caused by the different thickness of the crust, which varies with the NS mass (between ∼ 0.5 km for M=1.76 M⊙ and ∼ 1 km for M=1.10 M⊙). In the inner crust, the resistivity is likely dominated by electron-impurity scattering at low temperatures. As a consequence, the parameter Qimp controls the evolution at late time (t 10 5 yr), when the star is cold enough ). If instead of assuming an amor- Figure 6. Same as the upper panels of Fig. 3, but for models: A15 (first row), A14T (second), B14 (third) and C14 (fourth). In models B14 and C14, the magnetic field lines penetrate in the core, but they are not shown here.
phous pasta phase, a low value of Qimp in the pasta region is set, the dissipation at late times is much slower (compare models with Qimp = 0, 25, 100 in the pasta region for the same 1.4 M⊙ NS). In this low temperature, low resistivity regime, the magnetic field evolution becomes dominated by the Hall effect. The most extreme case (Qimp = 0 in the pasta region, cyan line in Fig. 5) shows a sort of stationary oscillatory modes (see also , with small twisted structures propagating meridionally. Another consequence of a low impurity parameter is that, since the magnetic field dissipation timescale is much longer than the star age, a similar asymptotic value of the magnetic field is reached for all NSs born as magnetars ). On the contrary, in models with a high Qimp in the pasta region, the magnetic field dissipation is maintained for millions of years.
Another important effect is that, because of the extra heat provided by the dissipation of magnetic energy, the dependence of the luminosity on other parameters is strongly reduced, compared with non-magnetised NS cooling models. In the bottom panel of Fig. 5, we show that the cooling curves of models with the same initial B 0 p = 10 14 G, but different masses or impurity content, are almost indistinguishable for at least 10 4 yr. Compared with non-magnetised models (Fig. 1), differences are strongly reduced under the presence of Joule heating. We also note that varying the superfluid gaps has a visible effect only when the effect of magnetic field decay is negligible, i.e., for low values of Qimp in the pasta region or weak magnetic fields.
Hereafter we will set the impurity paramter Qimp = 100 in the pasta region (in the regions of the crust where ρ > 6 × 10 13 g cm −3 ), leaving Qimp = 0.1 elsewhere in the crust, and the NS mass to 1.4 M⊙ for the rest of the discussion.

Dependence on initial magnetic field geometry and strength.
We continue the discussion by considering how the initial magnetic field configuration affects the evolution. In the top row of Fig. 6, we present the results for the evolution of model A15. The qualitative behaviour is the same as in model A14, but since the Hall timescale goes as ∝ 1/B, the dynamics is accelerated by a factor of ten. The more relevant role of the Hall term for stronger magnetic field leads to the formation of a discontinuity in B φ . In this particular case the current sheet is located at the equator (see the compressed lines and the "colliding" toroidal magnetic field of opposite sign in the second panel), where the dissipation is strongly enhanced ). Note that the average surface temperature is higher than for model A14 (Fig. 3). The second row of Fig. 6 shows the evolution of model A14T. The magnetic energy initially stored in the toroidal component is 98% of the total energy. Compared with model A14 (Fig. 3), there are significative differences. In particular, the symmetry with respect to the equator is broken because the vertical drift of the toroidal component acts towards the north pole. The initial poloidal magnetic field is distorted and bent at intermediate latitudes, with the formation of a small scale bundle of magnetic field lines, which has been proposed to be necessary for the radio pulsar activity (Geppert et al. 2013). This again locally increases the dissipation rate. At about 10 3 yr, the northern hemisphere is, on average, warmer than the southern hemisphere, and a characteristic hot polar cap is observed. Later, the geometry of the magnetic field becomes more complicated, with the formation of localized bundles of nearly radial magnetic field lines. This results in a temperature pattern with hot and cool rings.
In the third row of Fig. 6, the evolution of the coreextended configuration (model B14) is shown. The magnetic field lines are penetrating inside the core although, for clarity, the figure shows only the configuration in the crust (enlarged for visualization, see Fig. 2 for the real scale). The magnetic field in the core is almost frozen due to the high electrical conductivity of the interior (as discussed above). Some weak Hall activity is developed at the bottom of the crust, but the maximum value of the toroidal field generated is about one order of magnitude weaker than in model A14, and the poloidal lines do not suffer any significant bending. In addition, the reduced heat deposition in the crust results in a much cooler surface compared to all the other models, and similar to the low field cooling models.
Between the two extremes (mostly crustal currents vs. mostly core currents), we also consider the intermediate case of model C14, shown in the bottom row of Fig. 6. The presence of crustal currents activates the Hall dynamics and leads to a similar evolution as for model A14. The only relevant difference is that at very late times there will be a long-lasting magnetic field in the core.
In the top panel of Fig. 7 we compare the evolution of Bp for four models, all with B 0 p = 10 14 G. The difference between crustal-confined (models A14 and A14T) and the core-extended configuration (model B14) is evident. In the latter, the magnetic field is almost constant during 1 Myr. The configuration with an extremely strong initial toroidal component shows a larger luminosity (up to one order of magnitude at some age), due to the larger internal energy reservoir, but the higher temperatures also lead to a much faster decay of Bp than for model A14. Moderate values of toroidal magnetic fields (e.g., with half of the magnetic energy stored in) lead to an intermediate decay, between models A14 and A14T. In the hybrid case, model C14, Bp dissipates the crustal currents before oscillating around the value given by the currents circulating in the core. In the bottom panel of Fig. 7, we compare the luminosity for the same four models. Model B14 is similar to the non-magnetised case (Fig. 1), because most of the current is in the core and does not decay. On the other hand, when currents are totally or partially placed in the crust, extra heat is deposited in the outer layers, resulting in higher luminosities. One important result, not discussed in previous works, is that the Hall term causes the self-regulation of the ratio of the toroidal and poloidal components. The differences due to particular initial ratios are much less important than for the Ohmic case ), except when the initial toroidal component is significantly larger than the poloidal magnetic field. If the initial energy stored in poloidal and toroidal components are similar, the difference between models with and without initial toroidal fields is much smaller.
To conclude this section, we discuss the dependence of the surface temperature on magnetic field strength. In Fig. 8, we compare the surface temperatures at the pole and equator for B 0 p = 10 13 , 10 14 , 10 15 G. At early times, little difference is seen in the evolution of the temperature at the pole, regardless of the magnetic field strength. For the stronger magnetic fields, the polar temperatures can be kept high for longer times. The largest differences, however, are seen in the equatorial temperatures. For weakly magnetised models, the anisotropic blanketing effect of the magnetised envelope leads to much cooler equatorial regions. Conversely, for the strongly magnetised models, the high energy deposition rate by magnetic field decay is able to compensate this effect and the equator is actually warmer than the pole after a few hundred years. A dedicated work will explore the temperature anisotropies and their observational consequences (Perna et al. 2013).

Expected outburst rates.
The evolution of the magnetic field in the NS causes continuous stresses on the crust. Generally, magnetic stresses are balanced by elastic stresses. However, during the evolution, local magnetic stresses can occasionally become too strong to be balanced by the elastic restoring forces of the crust, which hence breaks, and the extra stored magnetic/elastic energy becomes available for powering the observed bursts and flares (Thompson & Duncan 1995. The maximum stress that a NS crust can sustain has been estimated with analytical arguments and through molecular dynamics simulations. In particular, Chugunov & Horowitz (2010) obtained the fit where Γ = Z 2 e 2 /aT is the Coulomb coupling parameter, a = [3/(4πni)] 1/3 is the ion sphere radius, ni the ion number density, Z the charge number, e the electron charge, and T the temperature.  combined results from the magnetothermal evolution of NSs ), with a computation of the magnetic stresses exerted on the NS crust. This allowed them to estimate the frequency and the location (depth and latitude) of starquakes, i.e. NS crust failures. The released energies were also estimated from the stored elastic energy. Their work showed that the classification of objects as AXPs, or SGRs, or high-B NSs, or 'normal' radio pulsars, etc. does not correspond to an underlying intrinsic physical difference: outbursts can occur also in objects outside of the traditional magnetar range, albeit with a lower probability. The follow up study by  further highlighted the importance that the toroidal magnetic field has on the NS observed phenomenology: two objects with similar inferred dipolar B-field (as measured by P andṖ ) can display a very different behaviour depending on the strength of the (unmeasured) internal toroidal component. The stronger the latter, the higher the luminosity of the object, and the more likely starquakes are to occur. In order to follow numerically the long term evolution of NSs, these works used an approximation to treat the Hall term, which was only included when studying the short term evolution of the magnetic field. In the current work, the above approximation has been released, and the effect of the Hall term is fully implemented. Hence, in the following we update their results by using a similar formalism to compute the outburst statistics from the results of our last simulations. In Fig. 9 we show the evolution of the predicted outburst rate for models A15 (red), a model of type A with B 0 p = 3 × 10 14 (black), and the model A14T (blue). While young NSs may undergo energetic outbursts every few years (for model A15) or tens of years (for the model with a weaker magnetic field), middle age sources are expected to show transient activity only every few thousands of years, and with less energy involved. Other models with initial lower magnetic fields, like model A14, do not give an appreciable event rate (less than 10 −3 events per year even for young NSs). However, in presence of strong toroidal magnetic field, model A14T, the outburst rate is strongly enhanced.
Assuming a NS birth rate of 10 −2 per year, there must be ≈ 10 4 NSs in the Galaxy with ages of 1 Myr. If 10% of them are born with B 10 14 G, a naive extrapolation of the estimated event rate at that age ( 10 −4 per year) leads to an outburst rate by old magnetars of one every few years. Therefore, we expect that more and more objects of this class will be discovered in the upcoming years. It is interesting to note that such old magnetars are expected to be detected only via their outburst activation, being too faint in the quiescent state even for very deep X-ray surveys. Since the launch of the Swift satellites in 2004, our capability of detecting and monitoring galactic transients has largely increased, and 5 new magnetars have been discovered in the past years, two of those being relatively old and displaying low inferred magnetic field: SGR 0418  and Swift J1822 Scholz et al. 2012).
On the other hand, the observational detection of an outburst event depends on the quiescent flux. The brightest sources, with L 10 35 erg/s, will experience a barely detectable flux enhancement, while for dimmer magnetars (L ∼ 10 33 − 10 34 erg/s) the flux increase by 2 or 3 orders of magnitude during the outburst and are much easier to detect .

THE UNIFICATION OF THE NS ZOO.
After having thoroughly reviewed the available observations in section §2, and having presented in §4 magneto-thermal evolution models that include the effects of magnetic field in the cooling history of NSs, we now proceed to compare the data with theoretical models.
The spin-down behavior of a rotating NS is governed by the energy balance equation relating the loss of rotational energy to the energy loss due to emission of electromagnetic radiation, gravitational radiation, or winds. In the following we will assume that the dominant contribution to the torque is electromagnetic radiation. The spin-down formula for force-free magnetospheres, and valid for oblique rotators (Spitkovsky 2006), gives the following relation where R is the NS radius, χ is the angle between the rotational and the magnetic axis, c is the speed of light, Ω = 2π/P is the angular velocity, and I is the moment of inertia of the star. In our model with M = 1.4M⊙ and R = 11.6 km, the moment of inertia is I = 1.33×10 45 g cm 2 . If a star with a different mass is considered, the value of I can change by ∼ 30% at most. Note that, compared with the standard formula for vacuum, IΩΩ = − B 2 p R 6 Ω 4 6c 3 sin 2 χ, the torque in Eq. (6) is larger by a factor of 1.5−3 compared with the orthogonal vacuum rotator and it is non-zero even for an aligned rotator (χ = 0).
We have numerically integrated Eq. (6) with the time dependent Bp(t) obtained from our simulations. We assumed an aligned rotator (constant χ = 0), and an initial period P0 = 0.01 s, low enough for the results not to be sensitive to the particular value of the initial period. The integrated period at a given age can change by a factor ∼ 2 if the uncertainties in the inclination angle are considered in the spin-down formula. Here we do not consider two possible effects that can contribute to the evolution of timing properties: the variation with time of the angle, χ = χ(t), and the growth of the superfluid region in the core during the early stages, resulting in a change of the effective moment of inertia, I = I(t).
In Fig. 10 we show the evolutionary tracks in the P −Ṗ diagram for a 1.4M⊙ NS with type A geometry and different initial values of the magnetic field strength, compared to the measured timing properties of X-ray pulsars. Asterisks mark different ages (t = 10 3 , 10 4 , 10 5 , 5 × 10 5 yr), while dashed lines show evolutionary tracks that the star would follow if there were no magnetic field decay. Bp is almost constant during an initial epoch, t 10 3 − 10 5 yr, which depends on the initial B 0 p : stronger initial magnetic fields decay before weaker ones. After that initial phase, tracks bend downwards due to the dissipation of the magnetic field under the combined action of the Hall effect with the large  Figure 10. Evolutionary tracks in the P −Ṗ diagram with mass and radius of our model, with B 0 p = 3 × 10 12 , 10 13 , 3 × 10 13 , 10 14 , 3 × 10 14 , 10 15 G. Asterisks mark the real ages t = 10 3 , 10 4 , 10 5 , 5 × 10 5 yr, while dashed lines show the tracks followed in absence of magnetic field decay. resistivity in the innermost part of the crust. Therefore, an asymptotic value of P is reached. The age at which this vertical drop in the track begins depends on the initial magnetic field, being roughly ∼ 10 5 yr for strong magnetic fields (B 0 p 10 14 G), and up to several Myr for weaker magnetic fields. The specific magnitude of the asymptotic value of P depends on the initial magnetic field, the resistivity in the innermost part of the crust, the mass of the star and the angle χ (see Eq. 6).
An inspection of the data reveals that there may be common evolutionary links among groups of isolated NSs. One of such links roughly includes nine magnetars (1E 1547, SGR 1627, SGR 0501, XTE J1810, CXO J1647, 1E 2259, 4U 0142 and the low B magnetars Swift 1822, SGR 0418) and five of the XINS (excluding only the faint source RX J0420 and RX J1605, that has no measure of P anḋ P ). For our particular model, this track corresponds to a NS born with B 0 p = 3 × 10 14 G. In the upper right corner, we identify a second group of eight extreme magnetars, characterized by a largerṖ (CXOU J1714, SGR 1900, 1E 1048, SGR 0526, CXOU J0100, 1RXS J1708, 1E 1841, SGR 1806. These would be consistent with NSs born with a higher initial magnetic field B 0 p ∼ 10 15 G. In light of these results, we can revisit the classical comparison between cooling curves and data reviewed in §4.1 with the predictions of the full magneto-thermal evolution, for a range of magnetic field strengths. In Fig. 11 we show the luminosity as a function of time for different values of the initial magnetic field up to B 0 p = 3×10 15 G, and we compare these cooling curves with the observational data. Compared with the non-magnetised cooling curves (Fig. 1), the most relevant difference is that the inclusion of the magnetic field allows to generally explain objects with high luminosities. Several authors already pointed out the observed correlation between inferred magnetic field and luminosity or temperature in magnetars (e.g., An et al. 2012). We are now able to confirm and quantify the trend. Magnetic fields above B 10 14 G are strong enough to noticeably heat up the crust and power the observed X-ray radiation. The cooling timescale for strongly magnetised objects is about one order of magnitude larger than for the weakly magnetised NSs.
For clarity of discussion, we now comment separately on the sources classified in three groups: • NSs with initial magnetic field B 0 p 10 14 G. This group includes all RPPs, two high-B PSRs (PSR J0726 and PSR J1119), and two of the XINS (RX J0420 and RX J1605). For these relatively low magnetic fields, the luminosity is not expected to be significantly affected by the presence of the magnetic field.
• NSs with initial magnetic field B 0 p 5 × 10 14 G. The 8 magnetars in the upper-right corner of Fig. 10.
We emphasize that this classification does not reflect intrinsic differences between groups; it is simply an arbitrary grouping that helps to highlight the evolutionary paths as a function of the initial magnetic field strength.
5.1 NSs with initial magnetic field 10 14 G.
As already discussed in Sec. 3, the standard cooling curves can account for all the observed sources by simply varying the star mass. We note that there is actually more freedom in the theoretical models than in the particular NS model shown in Fig. 11, in which the microphysics input (e.g. gaps) and the NS mass has been fixed. Considering the uncertainties in the inferred ages and luminosities, the cooling curves for weakly magnetised NSs are consistent with all the observational data, with a few particular cases which are worth discussing.
First, CCOs have very lowṖ (below the shown range), implying a weak external magnetic field. Due to the small spin-down rate, their periods have not changed appreciably since birth, and no information about their real ages can be inferred from timing properties. The weak inferred magnetic field apparently contrasts with the observed surface anisotropies and high luminosities (significantly higher than for standard RPPs) of the sources in the SNRs Kes 79 and Puppis A. Light-element envelope models (see Fig.  1) can reconcile their age and large luminosity, but not the large anisotropy. The latter can be explained by is the hidden magnetic field scenario (Young & Chanmugam 1995;Muslimov & Page 1995;Geppert et al. 1999;Ho 2011;Shabaltas & Lai 2012;, in which the object has actually a strong subsurface magnetic field, which has been screened by the initial fallback of the supernova debris onto the NS. The slow reemergence of the surface magnetic field happens on timescales of 10 3 − 10 5 yr. This scenario physically justifies the models with light-element envelope. On the other hand, a few sources (Vela pulsar, PSR B2334 and PSR J1740) show evidence for enhanced cooling by direct URCA processes (or alternative exotic fast cooling neutrino processes), as already discussed by many authors. With our equation of state and microphysics setup, this can happen only for high masses, M 1.5M⊙.
Comparing the cooling curves to the group of NSs with initial magnetic field B 0 p ∼ 1−5×10 14 G, we see that, although a few objects could still be compatible with non-magnetised curves, many others are clearly too luminous and need an additional heating source (see Fig. 11). Our results show that most of them can be reconciled with the theoretical models for a narrow range of initial magnetic fields, between B 0 p ∼ 3 − 5 × 10 14 G. Note also that, for the few cases where we have both kinematic and characteristic ages (e.g., some XINSs), the latter overestimates the real age. For this reason, it is likely that the two XINSs and SGR 0418, placed on the right hand side of the cooling curves, have in fact real ages ≈0.5-1 Myr.
For the most luminous object that would belong to this group, 4U 0142, we do not have any alternative estimate of the age. However note that it is quite similar in both timing properties and luminosity to 1E 2259, whose kinematic age inferred for the SNR CTB109 associated with 1E 2259 is about 10 4 yr (Castro et al. 2012), more than one order of magnitude smaller than the characteristic age. We note that for these two objects it is difficult to reconcile the observed timing properties and luminosity, even with more extreme models (very strong toroidal magnetic field). While their luminosity is compatible with very high magnetic fields, the timing properties are more consistent with initial B 0 p ∼ 3 × 10 14 G.
Interestingly, these are the two cases in which some weak evidence for the presence of fallback disks has been reported (see Wang et al. 2006for 4U 0142, and Kaplan et al. 2009for 1E 2259. In both cases, the IR measurements are consistent with passive disks, i.e. disks in which the viscosity has been substantially reduced when they become neutral, as expected after ∼ 10 3 − 10 4 yr (Menou et al. 2001). These disks are no longer interacting with the pulsar magnetosphere, and hence the measured characteristic age would again overestimate the actual NS age, since in the past there was a higher torque than there is today.
For these sources, the cooling curves with iron envelopes for B 0 p = 10 15 G are barely reaching their large luminosities, which may be an indication that these, still young, sources possess light-element envelopes, or perhaps these objects are simply born with even higher magnetic fields of a few times 10 15 G.
We need to make some more considerations about these extreme magnetars. First, we note from the evolutionary tracks in the P −Ṗ diagram that this group does not appear to have older descendants, in contrast to the first group of magnetars, which evolves towards the XINSs and SGR 0418. The expected descendants of the extreme magnetars should have periods of few tens of seconds and would be bright enough to be seen. No selection effect in X-ray observations would prevent to observe such sources.
Second, four of them (SGR 1900, SGR 1806, 1E 1841 and 1RXS J1708) show a strong, non-thermal component in the hard X-ray band (20 − 100 keV), whose contribution to the luminosity is up to 10 36 erg/s (Kuiper et al. 2004(Kuiper et al. , 2006Götz et al. 2006). Together with 4U 0142 and 1E 2259, they are the only magnetars showing this hard X-ray emission persistently (i.e., not connected to outburst activity as in SGR 0501 and 1E 1547; Rea et al. 2009;Enoto et al. 2010). We stress again that the tails seen in soft and hard X-ray are ultimately due to the magnetospheric plasma, which can provide a significant amount of energy to the seed thermal spectrum emerging from the surface.
Third, in most of the objects for which the kinematic age is available, a reverse, unusual mismatch with characteristic age, τ k > τc, is observed (see Table 1). This implies that the current value ofṖ is larger than in the past. Note also that SGR 1900 and SGR 1806 are the magnetars with the most variable timing properties (see Table 2 of McGill catalog).
This facts, i.e. the possible overestimate of both magnetic torque and thermal luminosity, and the absence of descendants, could point to some additional torque, apart from the dipole braking, acting temporarily during the early stages. One possibility discussed in the literature is wind braking (Tong et al. 2013), which can be effective during some epochs of the NS life, explaining the anomalous high values ofṖ of some magnetars. Another compelling possibility, supported by the timing noise, is the contribution of Figure 11. Comparison between observational data and theoretical cooling curves. Models A with B = 0, 3 × 10 14 , 3 × 10 15 G are shown for Fe envelopes (solid) and light-element envelopes (dashed).
magnetosphere. In such scenarios, this group of most extreme magnetars would be born with B 0 p 10 15 G, but they would experience extra torque and luminosity due to temporary effects.

SUMMARY
We have presented a comprehensive study of the magnetothermal evolution of isolated NSs, exploring the influence of their initial magnetic field strength and geometry, their mass, envelope composition, and relevant microphysical parameters such as the impurity content of the innermost part of the crust (the pasta region). Our state-of-the-art, 2D magneto-thermal code, is the first able to work with large magnetic field strengths, while consistently including the Hall term throughout the full evolution. The Hall term plays a very important role in the overall magnetic field evolution, strongly enhancing the dissipation of energy over the first ∼ 10 6 yr of NS life, with respect to the purely resistive case. This is due to two main effects: the generation of smaller structures and currents sheets, and the gradual compression of currents and toroidal magnetic field towards the crust/core interface. Hence, the rate of magnetic field dissipation strongly depends on the resistivity given by the amount of impurities in the innermost region of the crust. A highly impure or amorphous inner crust produces a significant increase in the magnetic field dissipation on timescales 10 5 yr. We also found that, while for weakly magnetized objects low mass stars (M 1.4 M⊙) are systematically brighter than high mass stars, this separation is smeared out for highly magnetized stars. We confirm that light-element envelopes are able to maintain a higher luminosity (up to an order of magnitude) than iron envelopes for a long period of time, 10 4 yr, regardless of the magnetic field strength.
The initial magnetic field configuration plays an important role in the observational properties of the NS. If the currents sustaining the magnetic field flow in the core, their dissipation is negliglible, comparable with models in which (most of) currents flow in the crust. In particular, the presence of an initial strong dipolar, toroidal component in the crust breaks the symmetry with respect to the equator, resulting in a warmer hemisphere, and a characteristic hot ring replaces the traditional hot polar cap. If the magnetic field is allowed to penetrate into the core (rather than remaining confined within the crust), then the reduced heat deposition in the crust results in a much cooler surface compared to the case in which the magnetic field lives in the crust only.
The estimated outburst rate, resulting from breaking of the crust by the strong magnetic stresses, is found to be an increasing function of the initial magnetic field strength and a decreasing function of age (in agreement with the results of . While this qualitative trend is consistent with the data (younger and more strongly magnetized objects tend to have more frequent outbursts), a more quantitative comparison between the simulations and the observations is not possible at this stage, due to the lack of sufficient statistics in the data.
For a given choice of the model parameters summarized above, our simulations allow us to predict, for any age of the NS, the distribution of surface temperature, luminosity, or timing parameters. In order to compare theoretical models with observations, we have studied all currently known isolated NSs with a detected surface thermal emission. These objects are representative of different observational classes, including the classical magnetars (AXPs and SGRs), high-B pulsars, standard radio pulsars, XINSs, and CCOs. The sample includes only the objects with clearly seen thermal emission. We have performed for all of them an homogeneous, systematic analysis of the X-ray spectra, inferring the thermal luminosity, and trying to reduce at the minimum the systematic effects due to different data analysis and modelling. We published on-line the results of our analysis, including detailed references for every source 1 .
A comparison between a range of theoretical models and the observations (both timing and luminosities), has shown that, for the weakly magnetised NSs (Bp 10 14 G in the P −Ṗ diagram), the magnetic field has little effect on the luminosity. These objects, of which the RPPs are the most notable representatives, have luminosities which are compatible with the ones predicted by standard cooling models, with the dimmest ones in the sample (Vela pulsar, PSR B2334 and PSR J1740) requiring an iron envelope and an NS mass 1.4M⊙.
The bulk of the magnetars, with inferred magnetic fields identified in the range Bp ∼ a few ×10 14 G in the P −Ṗ diagram, displays luminosities generally too high to be compatible with standard cooling alone. The magneto-thermal evolutionary models with B 0 p ∼ 3 − 5 × 10 14 G can account for their range of luminosities at the corresponding inferred ages. As these objects evolve and their magnetic fields dissipate, their observational properties (both timing and luminosities) appear compatible with those of the XINSs.
Finally, the most extreme magnetars, endowed with magnetic fields of strength 10 15 G as inferred from their timing parameters, are characterized by the highest luminosities, barely compatible with the B 0 p = 10 15 G cooling curve with an iron envelope. However, for the same initial magnetic field, a light-element envelope is able to account for the luminosity of even the brightest objects. Note that most of these extreme objects show hard X-ray, non-thermal emission and unreliable estimate of the dipolar magnetic fields from timing properties, due to the large timing noise. We suggest that, at least for some of them, the magnetosphere plays an important role in influencing timing properties and luminosity of these objects, whose magnetic field may be not so larger than the bulk of magnetars.
In summary, our magneto-thermal simulations have begun to paint a unified picture of the variety of the observational properties of isolated NSs, and their evolutionary paths. Our results can account for the overall phenomenology, whereas in-depth testing of each specific magnetic field configuration will require detailed spectral and timing modeling of the thermal component of each source. This is reserved to future work.
Further, deeper observations and discoveries of magnetars will allow in the next years to refine the comparison with data and constrain the current uncertainties of the theoretical models, such as the properties of the inner crust and the initial magnetic field geometry.