3D code for MAgneto-Thermal evolution in Isolated Neutron Stars, MATINS: thermal evolution and lightcurves

The thermal evolution of isolated neutron stars is a key element in unraveling their internal structure and composition and establishing evolutionary connections among different observational subclasses. Previous studies have predominantly focused on one-dimensional or axisymmetric two-dimensional models. In this study, we present the thermal evolution component of the novel three-dimensional magnetothermal code MATINS (MAgneto-Thermal evolution of Isolated Neutron Star). MATINS employs a finite volume scheme and integrates a realistic background structure, along with state-of-the-art microphysical calculations for the conductivities, neutrino emissivities, heat capacity, and superfluid gap models. This paper outlines the methodology employed to solve the thermal evolution equations in MATINS, along with the microphysical implementation which is essential for the thermal component. We test the accuracy of the code and present simulations with non-evolving magnetic fields of different configurations (all with electrical currents confined to the crust and a magnetic field that does not thread the core), to produce temperature maps of the neutron star surface. Additionally, for a specific magnetic field configuration, we show one fully coupled evolution of magnetic field and temperature. Subsequently, we use a ray-tracing code to link the neutron star surface temperature maps obtained by MATINS with the phase-resolved spectra and pulsed profiles that would be detected by distant observers. This study, together with our previous article focused on the magnetic formalism, presents in detail the most advanced evolutionary code for isolated neutron stars, with the aim of comparison with their timing properties, thermal luminosities and the associated X-ray light curves.


INTRODUCTION
Neutron stars (NSs) are born from the gravitational collapse of a stellar core during a supernova explosion.During the first minute of its life as a proto-NS, the star is hot (∼ 10 11 K), opaque to neutrinos, and undergoes a shrinking of its radius and a deleptonization, while neutrinos diffuse outward (Prakash et al. 2001).The proto-NS cools down until, at a temperature of ∼ 10 10 K, the matter becomes transparent to neutrinos.During the subsequent secular evolution, the NS temperature drops due to different processes, which are dominated initially by the leakage of neutrinos, during the so-called neutrino cooling era, and later by the emission of photons from their hot surface, throughout what is known as photon cooling era (e.g.Page et al. 2004Page et al. , 2006;;Yakovlev & Pethick 2004).
Studying the precise cooling history of NSs can in principle pro-★ E-mail: stefano.ascenzi@gssi.itvide us with a wealth of information.For instance, during the neutrino cooling era, various neutrino emission processes, characterized by a different dependence on the temperature, coexist.The presence of a given mechanism and its relative importance with respect to the others is ultimately determined by the microphysical conditions occurring in the interior of the NS, such as the stellar composition, the onset of a superfluid phase transition, or/and the presence of strong magnetic fields.Consequently, the cooling history of NSs can significantly enhance our understanding of their internal state.
Moreover, there exists a tension in the fact that the estimated galactic rate of Core-Collapse Supernovae (CCSNe) is smaller than the sum of the birth rates of the different NS populations (Keane & Kramer 2008).This poses a problem for the scenario that attributes the formation of all the NSs to CCSNe, as they alone could not account for the entire population of NSs in our Galaxy, if we posit that the presence of phenomenological sub-classes inherently signifies distinct types of NSs.Various solutions have been proposed to this issue, one of which suggests that different NS populations are related through evolutionary paths (see e.g.Viganò et al. 2013).If this holds true, we only need to consider a single birthrate for these populations, which helps resolve (or at least alleviate) the tension with CCSNe rates.Indeed, magnetothermal cooling models are capable of establishing this evolutionary link (see Pons & Viganò 2019, for a review).
However, modeling the thermal evolution of NSs is not an easy task.The process is determined by several ingredients, which are at present not very well understood, such as the model of superfluidity, the composition of the core, the neutrino processes, and the model and composition of the blanketing envelope, which lies above the outer crust.Moreover, the intense magnetic field in the NSs interiors influences the transport properties of the star, making the conductivity anisotropic.Therefore, an accurate evolution of a strongly magnetized NS requires a multi-dimensional approach.Furthermore, the dissipation of the electric currents, which leads to a decay in the magnetic field, acts like a heating source in the stellar crust.The magnetic field decay, in turn, is regulated by magnetic diffusivity (in addition to being driven by the Hall effect in the crust and by ambipolar diffusion and other little-understood mechanisms in the core), which is a parameter that depends on the local temperature.It hence follows that the equations describing the thermal and the magnetic evolution cannot be considered independently, but must be solved together in a coupled magneto-thermal framework (see Pons & Viganò (2019) for a comprehensive review).
Extensive investigations spanning several decades have been devoted to unraveling the magnetothermal evolution of NSs.Initial attempts primarily delved into cooling phenomena within 1D models, with limited consideration for the magnetic field's impact (Tsuruta & Cameron 1965;Yakovlev & Urpin 1981;Page & Baron 1990;Page et al. 2004Page et al. , 2009;;Yakovlev & Pethick 2004;Kaminker et al. 2008;Potekhin & Chabrier 2018;Beznogov & Yakovlev 2015a,b).Advancing beyond this, subsequent studies ventured into axisymmetric, two-dimensional (2D) calculations; however, the first studies did not consider a self-consistent coupling in the magneto-thermal evolution, but rather assumed a prescription for the temperature evolution and solved for the magnetic field one (Pons & Geppert 2007) or vice versa (Aguilera et al. 2008).
A step forward occurred with the works by Pons et al. (2009); Viganò et al. (2012Viganò et al. ( , 2013)), who presented a consistent treatment of the coupled magnetothermal evolution in a 2D setting.In particular, the later works successfully accounted for the Hall term in the induction equation, which plays a fundamental role in transferring magnetic energy to smaller spatial scales where the dissipation is more effective (Pons & Geppert 2007), but that had posed a numerical challenge for previous codes.
While all these works have been focused on the magnetic evolution in the crust alone, progress has also been made in the study of the more complex evolution in the core (Graber et al. 2015;Ofengeim & Gusakov 2018;Gusakov 2019;Castillo et al. 2020;Dommes et al. 2020;Wood & Graber 2022).There, a multifluid, multiscale dynamics needs to be studied, which implies the presence of ambipolar diffusion, but with non-trivial complications due to superconductivity and highly uncertain coupling coefficients between the different fluid components.The core evolution is an open problem from a theoretical point of view: despite the advances in understanding the relevant timescales and in implementing specific ingredients (such as ambipolar diffusion for non-superfluid/superconducting matter), no fully self-consistent induction equation nor simulations with all the ingredients exist yet.Nonetheless, the general consensus is that the crustal timescales are generally much shorter than the ones for the core.Therefore, while the core evolution will determine the field in Gyr-old neutron stars (like millisecond pulsars and low-mass X-ray binaries), the coupled magneto-thermal evolution in the crust, object of this code, has a much more direct effect on the observables in young neutron stars (magnetars in particular).
The aim of this paper is to introduce the thermal part of MATINS (MAgneto-Thermal evolution of Isolated Neutron Star), the magnetic component of which has been already presented in Dehman et al. (2023a), with coupled magnetothermal applications in Dehman et al. (2023c).MATINS differs from PARODY in three main aspects: 1) while PARODY is a pseudo-spectral code, namely it employs a finite grid in the radial direction and a spherical harmonic expansion in the angular directions, MATINS uses a finite volume scheme; 2) while in PARODY the stellar background is accounted for by assuming an analytical radial profile for the electron number density (or the chemical potential, equivalently) and by defining a radius for the star and the crust-core interface, in MATINS the stellar background and the composition is obtained by self-consistently solving the Tolman-Oppenheimer-Volkoff (TOV) equation with the option of choosing between different equation of states (EOSs) among the ones present in the public database CompOSE (CompStar Online Supernovae Equations of State; Typel et al. 2015;Oertel et al. 2017) 1 .Moreover, solving for the structure allows us to compute the spacetime metric within the star, and thus include self-consistently the general relativistic effects in the magnetothermal evolution.These effects in PARODY are neglected; 3) The microphysics in PARODY is described by analytic approximations which describe the relevant quantities (e.g.heat capacity, electrical and thermal conductivities, neutrino emissivity) with simplified dependency on the temperature and the magnetic field.In particular, these analytic prescriptions aim to mimic the contribution of the electrons on the conductivity and heat capacity, which however is not always the dominant contribution.In MATINS instead we include all the known contributions, through a numerical treatment based on the public code developed by A. Potekhin2 .This code allows us to include a detailed treatment of the relevant quantities including the contribution of different species (e.g.electrons, ion lattice, free neutrons) to them.Furthermore, in MATINS a proper treatment of the superfluidity is implemented, while in PARODY this phenomenon is neglected.
The paper is organized as follows: in Sec. 2 we introduce the thermal evolution equation in the MATINS code, further focusing on the grid employed and on the description of the microphysical setup.Sec. 3 is devoted to the description of the tests we performed to assess the reliability of the code.In Sec. 4 we present several runs of MATINS with and without the magnetic field evolution.In this section, we also introduce a ray-tracing code that takes as input the NS's surface temperature map obtained by MATINS to produce phase resolved spectra that a distant observer would detect, and present preliminary applications of this code.Finally, in Sec. 5 we summarize our work.

THERMAL EVOLUTION IN MAGNETIZED ISOLATED NEUTRON STARS
In this section, we present the setup of the thermal evolution part of the MATINS code.First, we introduce the thermal evolution equation.Then, we briefly describe the cubed sphere grid that we employ in MATINS.Finally, we illustrate two different microphysical setups that we have included in the code.For specific details on how the thermal evolution equation is discretized and solved in the code, we refer the reader to Appendix B and Appendix C. The boundary conditions are discussed in Appendix D.

Background Stellar Model
The stellar model is computed by the equation: which describes hydrostatic equilibrium assuming an internal spherically symmetric and static metric: where   ( ) is the lapse function, describing the gravitational redshift, and Ω is the solid angle.The functional form of  with radius depends on the function () (which represents the mass distribution inside the star in the Newtonian limit), since In the Eq. 1,  and  denote the pressure and the density, respectively.The lapse function, instead, is obtained by solving the  component of the Einstein field equations, which, assuming the metric in Eq. 2, can be written as along with the boundary condition  2 () = 1 − 2 /( 2 ), where  is the star gravitational mass and  is the radius of the star.In all the equations mentioned above, the constants  and  denote the speed of light in vacuum and the gravitational constant, respectively.Equations 1 and 4 need to be supplemented with the equation for the function (): (5) Equations 1, 4 and 5 are known as TOV equations (Tolman 1939;Oppenheimer & Volkoff 1939).
Solving the backgroud structure of the star via TOV equation requires the use of an EOS for cold dense matter, which gives the pressure and composition as a function of the density, along with the central pressure  0 .Currently, MATINS includes different cold EOS from the CompOSE database.The central pressure  0 is an input parameter which regulates the star mass .
Finally, above the crust lies a thin (∼ 10 − 100 m) liquid blanketing envelope layer.The steep radial gradients of density, pressure and temperature in this region imply an outward decrease of the local thermal timescales by orders of magnitude, which make it computationally unfeasible to follow the long-term (up to Myr) evolution.For this reason, as in all cooling codes, MATINS solves the equations for the core and the crust, down to a threshold pressure given in input and here set (as usual in these studies) to a corresponding density of ∼ 10 10 g/cm 3 .For the uppermost layers, MATINS allows the use of different envelope models (see also Appendix D), which relate the temperature at the star surface   to the temperature at the threshold pressure (the base of the envelope)   , and that may depend, eventually, also on the local magnetic field B.
The dependence of the results on the envelope models have been presented for instance in Potekhin et al. (2015) and Dehman et al. (2023b), and we will briefly show some results below.

Heat Diffusion Equation
The heat diffusion equation determines the evolution of the thermal internal structure and surface luminosity.It can be written as follows: where  is the local temperature,  v is the heat capacity per unit volume,  represents the sources or losses of energy per unit volume, and F is the heat flux.From the lapse function and the temperature, we can define the redshifted temperature T as T ≡   .
The source/loss term on the right hand side can be written as the sum of two contributions  =  ℎ −   .While   represent the neutrino emissivity per unit volume, the term  ℎ represent the heating rate per unit volume due to Ohmic dissipation of the electric current (i.e. the Joule heating), which writes as: where  is the electric conductivity, which has in the crust a typical value in the range  ∼ 10 22 − 10 25 s −1 , and J is the current density obtained, as standard in magnetohydrodynamics, by the fourth Maxwell's equation neglecting the displacement current: It is worth noticing that other heating mechanisms can in principle be at play in the stellar crust, e.g.crust-cracking (Baym & Pines 1971;Cheng et al. 1992), non-equilibrium reactions (e.g. Haensel 1992;Reisenegger 1995;Fernández & Reisenegger 2005;Flores-Tulián & Reisenegger 2006;Reisenegger et al. 2007;González-Jiménez et al. 2015), superfluid vortex creep motion (e.g.Alpar et al. 1984;Umeda et al. 1993;Larson & Link 1999;Schaab et al. 1999); however, they are not considered in the present work.
Concerning the loss of energy due to the emission of photons at the surface, we include this contribution by the boundary conditions, as described in Appendix D.
Since the thermal conductivity is generally dominated by the contribution of electrons, under the presence of strong magnetic fields it becomes anisotropic, with its contribution in the direction orthogonal to the field  ⊥ quenched with respect to its contribution along the field lines  ∥ .It is common to treat the problem in the relaxation time approximation, where the ratio between the two quantities is expressed by the following formula (Urpin & Yakovlev 1980): where the parameter    0 is the product between the characteristic electron relaxation time  0 and the electron gyrofrequency   , which is written as: where  and  *  are the electron charge and effective mass, respectively, and B is the local magnetic field.This approximation shows that    0 is the parameter governing the heat conduction: when    0 ≫ 1, either because the magnetic field is strong or because the electron collision rate is low, the electrons are tightly anchored to the field lines, and slip along them, while the transverse motion is strongly inhibited.This results in a reduction of the thermal conductivity in the direction orthogonal to the field lines.Vice versa, when    0 ≪ 1, the magnetic field can only poorly constrain the electron motion, and the conductivity becomes isotropic ( ⊥ →  ∥ ).
According to Eq. 9 and following Pérez-Azorín et al. (2006), we can write the heat flux as: where b ≡ B/||B|| is the unit vector in the direction of the magnetic field.Three terms contribute to the heat flux: a term parallel to the temperature gradient, which is the same appearing in Fourier's law (F ∝ ∇ ∇ ∇, valid in the absence of magnetic field); a term parallel to the magnetic field, which is due to the fact that electrons tend to move along the magnetic field lines; and a last term perpendicular to both the previous ones, known as Hall term, which is due to the drifting motion of electrons through field lines due to the presence of the temperature gradient.The magnitude of    0 regulates which of these terms gives the dominant contribution.In the limit where    0 is negligible, the first term is the dominant, while when    0 ≫ 1, the second term dominates.
From Eq. 11 we can appreciate that the flux shows a linear dependence with respect to the temperature derivatives, such that it is possible to re-write the previous equation in matrix form: where we assumed the Einstein sum convention of repeating indices.
The explicit form of the matrix    and its derivation is reported in Appendix E.

Grid
In our code, we use the cubed sphere coordinates as described by Ronchi et al. (1996), implemented in MATINS by Dehman et al. (2023a), and successfully used in several other contexts like geophysics (Breitkreuz et al. 2018;Ding & Wordsworth 2019;van Driel et al. 2021), general relativity (Lehner et al. 2005;Hébert et al. 2018;Carrasco et al. 2018Carrasco et al. , 2019) ) and magnetohydrodynamics (Koldoba et al. 2002;Fragile et al. 2009;Hossein Nouri et al. 2018;Wang et al. 2019;Yin et al. 2022).The spherical volume of the star is characterized using a radial coordinate  and two angular coordinates (, ).Each spherical surface of constant  is divided into six patches, which can be pictured as the six faces of a cube that have been inflated to adhere to a sphere.The six patches are identical, non-overlapping, and, as in a cube, each of them is bounded by four patches.A representation of the cubed sphere grid is reported in Figure 1.The angular coordinates  and  are defined differently in each patch, such that they are not singular and assume values in the range [−/4, /4].Moreover,  and  are non-orthogonal everywhere, except at the center of each patch, which implies non-diagonal terms in the metric, i.e. additional terms in the mathematical operators used in the discretized equations (gradients, scalar and cross products).This grid presents two main advantages: first, it has a radial coordinate, which is desirable since the background is spherically symmetric (see Sec. 2.1).Moreover, the problem we aim to solve is characterized by strong radial gradients, such that a radial coordinate allows us to refine the resolution more in the radial direction than in tangential ones.Second, the coordinates are regular everywhere, in contrast, for example, with the spherical ones, which are singular along the axis.The metric, the line, surface and volume elements and the differential operators for the metric are reported in Appendix A and in more detail by the accompanying paper Dehman et al. (2023a).
In order to consider the derivatives, we use one layer of ghost cells lying behind each edge in the  and  directions.The values of the quantities at these ghost points are obtained by linear interpolation from the cells of the neighbouring patches.This prescription is explained more extensively in Appendix C.
Hereafter, the resolution of the grid will be denoted by the number of radial cells   and the number of angular cells per-patch   , where   is the same for the  and  directions.In line with its 2D predecessors (Viganò et al. 2012(Viganò et al. , 2021)), in MATINS the thermal grid has half the resolution of the magnetic grid.Specifically, the temperature is evolved only at the center of the thermal grid, while the magnetic field is defined and evolved also at its edges and interface centers.Hence, the magnetic field discretized locations are 8 times more than the temperature ones.This choice is guided by two considerations: firstly, the computational needs, since the cost of the implicit scheme adopted for the thermal evolution (see Appendix B) scales non-linearly with the total number of points, while the magnetic evolution is linear.Secondly, because of the nature of the equations: the non-linear Hall term in the induction equation naturally create small structures and a cascade (Dehman et al. 2023a), something not present, by definition, in the diffusion of the heat which thus only needs a more moderate resolution.

Microphysics
The stellar microphysics enters the heat diffusion equation through the heat capacity per unit volume  v , the source term , the thermal conductivity  ⊥ and the    0 parameter.In MATINS these quantities are defined at the center of the cells.To compute the value of  ⊥ and    0 at the cell interfaces, needed to calculate the flux (see Appendix E), we interpolate their values at the cell centers.We tested a linear interpolation on the variables  ⊥ and    0 , as well as on their logarithm.We chose the latter since we consider it as more suitable to describe quantities that can exhibit orders of magnitude variations among adjacent cells.Moreover, our numerical experiments show slightly higher stability of the code, although, in any case, the differences in the results are negligible and decrease with increasing resolution.
The microphysics setup employed by the MATINS code is the same as the one already used in the previous 2D code (Viganò et al. 2013(Viganò et al. , 2021)), where the microphysical quantities are calculated numerically exploiting the public code released by A. Potekhin.A detailed review of the microphysics can be found in Potekhin et al. (2015), while here we provide a summary description of the most important microphysical process contributing to the macroscopic quantities entering in Eq. ( 6).
As already described in the previous Sec.2.2, the thermal conductivity under the presence of a magnetic field is anisotropic and it is described by a tensor, whose components are calculated using Potekhin's public code3 .In the core, the conductivity, dominated by electrons, is very high.In the crust, the conductivity is lower, and it is mostly influenced by the interaction between electrons, impurities in the lattice, lattice phonons, and phonons of the neutron superfluid and non-superfluid neutrons (see Fig. 4.5 of Viganò 2013, for a comparison between the different contributions).In the direction of the magnetic field, the electrons are expected to be the main contributor to the heat transfer, while in the orthogonal direction, the phonons can become dominant if the magnetic field is particularly intense.Finally, in our analysis, we neglect the quantizing effects due to the gradual fillings of Landau levels, which leads the thermal conductivity to oscillate with varying density around the classical value.Although these effects are considered in the Potekhin code, their inclusion is particularly expensive from a computational point of view, and since their contribution becomes important only at densities lower than that of the crust (such as in the external envelope), we can safely neglect them.
The heat capacity per unit mass in the outer crust has contributions from the degenerate electron gas and the ion lattice, while in the inner crust also the neutron gas contributes if the temperature is above the neutron superfluidity critical value.For its calculation we employ the public code by Potekhin4 , properly designed for a strongly magnetized, fully ionized electron-ion plasma.
The neutrino emissivity dominates the stellar cooling in the first 10 4 −10 5 yr (neutrino cooling era).After that the drop in temperature quenches the neutrino emission and the star cools mainly through the emission of thermal photons from the surface (photon cooling era).In this setup the neutrino emissivity is described by the formulae provided in Table 1 of Potekhin et al. (2015).
The microphysical setup of our code implements also superfluidity correction to the previous quantities for neutrons (singlet state) in the crust and for neutrons and protons in the core (triplet and singlet states, respectively).The energy gap and the critical temperature are approximated as a function of Fermi momenta by the parametrization provided in Kaminker et al. (2001), whose fit parameters depend on the chosen superfluid gap model.In MATINS, we provide the possibility to choose between the same gap models present in the 2D code of Viganò et al. (2021).The models and the relative parameter values are listed in Table II of Ho et al. (2015).Throughout this work, we employ the same gap model choice, corresponding to the model indicated by Ho et al. (2015) as SFB (Schwenk et al. 2003) for the neutron in the crust, TToa (Takatsuka & Tamagaki 2004) for the neutrons in the core, and CCDK (Chen et al. 1993;Elgarøy et al. 1996) for the protons in the core.

Timestep and workflow
Since the cooling timescale is highly dependent on the temperature itself (mostly due to the neutrino processes), it increases by orders of magnitude in the long term.Therefore, in MATINS the timestep is set proportional to the time itself.Users can customize this timestep by specifying the constant of proportionality, with default value set at 0.1.Additionally, the timestep has a minimum and a maximum value, that typically we set as Δ min = 10 −2 yr and Δ max = 10 4 yr, respectively.With this choice of parameters, MATINS completes a simulation with 1 Myr as total physical time in ∼ 250 thermal evolution timesteps.
The workflow of MATINS closely resembles the one presented by Viganò et al. (2021) in the predecessor 2D magnetothermal code.For a detailed visual representation, we refer the reader to Fig. 2 in their work.At the beginning of the simulation, we import the input parameters and utilize them to construct the stellar background.Additionally, we initialize the temperature and the magnetic field.Subsequently, we initiate the thermal evolution loop, wherein, at each timestep, the microphysics of the star (thermal and electric conductivities, heat capacity, neutrino emissivity) is updated using the new values of temperature and magnetic field.
The magnetic evolution loop is nested within the thermal evolution loop.During each thermal timestep, the magnetic field undergoes evolution, involving the calculation of currents and the Joule heating term, which enters the source term on the right-hand side of Eq. 6.If the magnetic field evolution is disabled, the Joule heating is switched off.Once the nested magnetic loop is complete (typically, a magnetic timestep of fraction of year, meaning hundreds or thousands per each thermal timestep), the temperature is updated, and a new cycle begins.

TESTS
In this section, we describe two problems employed to test MATINS: a test with a simplified microphysics and a known analytical solution (hereafter benchmark test), and a test with realistic microphysics.

A benchmark test
In this section, we describe a benchmark test to MATINS: a problem with a known analytical solution presented in Pérez-Azorín et al. (2006).It consists of a non-stratified spherical shell with inner radius   and outer radius   , with a uniform, non-evolving magnetic field oriented along a given axis.The model assumes a uniform conductivity and heat capacity.Additionally, the general relativistic corrections are neglected.No source term is included, thus the problem tests the anisotropy of the heat diffusion.This problem admits the following analytical solution: where  0 and  0 are the initial values of temperature and time,  is the polar angle calculated from the axis centered on the star and parallel to the magnetic field direction (see Ronchi et al. (1996); Dehman et al. (2023a) for details about transformation between spherical and cubed-sphere coordinates).The problem is expressed in dimensionless physical units, where  v =  ⊥ =  0 = 1,   = 5,   = 10,  0 = 100,    0 = 10.With this choice of units the diffusion timescales, defined as We choose the simulation timestep Δ = 0.01 ≪  diff .Given the unconditional stability of the implicit method we use to solve the heat diffusion equation, the timestep is not subjected to the Courant condition or any similar limitation.As such, the timestep is the same for all the resolutions employed here.For this test, instead of the usual boundary conditions (described in Appendix D), we impose, at each time , the analytical value of  given by Eq. ( 13) at the innermost and outermost radii.
First, we tested our code on this problem by choosing the orientation on the uniform magnetic field along the -axis, such that the magnetic field is perfectly radial at the centre of the north and south patches of our cubed-sphere grid.We performed three simulations at three different angular resolutions: a low (  = 7), medium (  = 14), and high resolution (  = 28).In all these cases the radial resolution is held fixed, at a value of   = 20.In Figure 2 we show at each time step the average L2 relative error, quantifying the volume-averaged deviations from the analytical solution, defined as where the sum is performed all over the computational domain.Here  num denotes the numerical solution, while  (  ,   , ) denotes the analytical solution calculated at the center of the i-th cell.
The results are reported with circle symbols in Figure 2. The label "centre" indicates that in these simulations the magnetic field is radial at the center of a patch.As one can see, improving the angular resolution by a factor of 2, improves the average L2 by almost an order of magnitude.With respect to the three cases with a radial field at the patch's centre (blue, orange and green curves), we note that, although L2 decreases with increasing angular resolution, the increase of L2 in time is faster at higher resolution.This behavior can be ascribed to the fact that at higher resolution the size of the cell in the transverse direction becomes comparable to its radial size (which is fixed in the three cases), causing the radial resolution to constrain the numerical performance.Additionally, the error accumulates nonlinearly with increasing resolution, so that we do not expect the three curves to exhibit the same slope.Nevertheless, despite their faster growth, the L2 at medium and high resolution also exhibit a negative second derivative, indicating a deceleration in their growth rates over time, which is more pronounced in the high resolution case.
The result of this simulation for the high-resolution case is reported in Fig. 3 for three different times.As expected, we notice that the temperatures diffuse mainly along the direction of the magnetic field, while the diffusion in the orthogonal direction is suppressed.In Figure 4 we report also the equatorial cut at the final timestep of the simulation.From this figure, we can appreciate how the axial symmetry is maintained during the simulation.
Finally, we repeat the same test but orienting the symmetry axis of the problem in different directions, in order to check if we are able to obtain the same evolution of the system.We perform four 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0  additional simulations at medium resolution, in which the field is oriented along the -axis, the -axis, the center of the edge between an equatorial patch and the north pole patch, and the corner between two equatorial patches and the north pole patch.The result in terms of average L2 is shown in Figure 2, along with the result for the cases with the field along the -axis discussed before.Here, the errors for the field orientated along the patch centres (namely along the , , or -axis) perfectly overlap and correspond to the orange curve.In the other two cases where the field is radial at the centre of an edge between two patches (red crosses) or at the corner between three patches (purple diamond), the error is, at the end of the simulation, a factor ∼ 2 and a factor ∼ 8 smaller than in the previous cases with the same resolution, respectively.Additionally, also the curve's slopes appear different, even if in all cases we note a flattening in the final steps.The reason for these differences can be, at least in part, ascribed to the differences in the cells, both in terms of size and shape, of the cubed sphere grid within the patch.The higher error in the "centre" case, for example, can be motivated by the fact that the cells at the centre of the patch are the largest ones.Therefore, the hotcolumn (see Figure 3), in the "centre" case (orange curve in Figure 2) becomes less resolved than the other two cases (red and purple curves), leading to a larger L2 error5 .On the other hand, the "corner" case (purple curve) is associated with a lower error than the "edge" case (red curve), even though the corner cells are slightly larger than the cells at the edge's centre (∼ 77% and ∼ 70% of the size of the cell on the patch centre, respectively).In this case, we ascribe the lower error to the difference in the shape of the cells at the corner compared to those at the edge's center.Indeed, the former are more suitable to approximate the circular section of the hot column, than the latter.We thus conclude that the difference in the error among   the three cases is non-trivial, and is determined by a combination of different factors.

Time [diffusion timescale]
Note that, regardless of the (low) resolutions employed here and the magnetic orientation, the relative uncertainty, estimated by the square root of the L2 error, is in the range 1 − 4% and decreases with resolution, therefore perfectly acceptable.
Finally, with the aim of assessing how the accuracy of our numerical solution varies with respect to the angular position on the cubed sphere grid, we report in Figure 5 a map of the relative uncertainty, calculated at the final timestep, on a spherical surface cut at the (arbitrarily chosen) radius  = 0.25( out +  in ).We define the relative error as [ num  −  (  ,   , )]/max[ (  ,   , )], where the denominator is the maximum temperature on the chosen spherical cut.We use the maximum of the temperature in the definition of the relative error, instead of the value point-by-point, in order to avoid large error values where the temperature is low.We observe that the error is largest at the interfaces between the equatorial and the polar patches, but not at the edges between one equatorial patch and the next one.

Tilting the field in a realistic model
In the previous section we studied the effect of tilting the orientation of the magnetic field in the benchmark test.In this section, we repeat the exercise with a realistic model, which, unlike the benchmark test, includes a physical NS background, neutrino emission processes, a non-uniform magnetic field, general relativistic corrections, and the boundary conditions described in Appendix D. The model employed for this study is Sly4-M1.4-B14-L1 reported in Table 1, characterized by a crust-confined dipolar poloidal magnetic field configuration.This configuration is characterized by the presence of two antipodal magnetic poles and it is symmetric for rotation around the magnetic axis, defined as the axis passing through the centers of the poles and the center of the star.The temperature field is initially uniform, with a value of  = 10 10 K. Further details on realistic configurations are provided in the next section.Here we limit ourselves to analyzing the impact of the different orientations of the magnetic axis on the results.This is shown in Figure 6, where, in the top panel, the bolometric thermal luminosity of the star is shown for three different configurations: one configuration in which the magnetic axis passes through the center of the patches (blue points), one in which it passes through the center of the edge between two patches (orange points), and one in which it passes through the corner between three patches (green points).From this figure, we can appreciate how the three curves appear indistinguishable.In the lower panel, we plot the relative error of the second and third configurations, taking the first one as a reference.We can see that before ∼ 3 × 10 5 yr for the edge configuration and ∼ 6 × 10 5 yr for the corner configuration, the error remains below 1 %, for a total simulation time of 1 × 10 6 yr.Beyond this time, the increase of the relative error, which almost reaches 10 % at the end of the simulation, is likely caused by the rapid decrease in luminosity.Even including the increased error at very late times, the simulations are in good agreement with each other well within the observational uncertainties.This test shows that the results keep their consistency under rotation of the magnetic field not only in the simplified benchmark test, but also for a realistic model.

RESULTS
Let's now turn to the application of MATINS to a realistic NS structure, under different magnetic field configurations.All the initial magnetic field configurations here considered are confined to the crust as discussed more in detail in Dehman et al. (2023a).Unless stated otherwise, we run our simulations for a total time of 1 Myr.As a benchmark for our study, we consider a NS characterized by the unified (i.e.crust + liquid core) EOS SLy4 (Douchin & Haensel 2001), which assumes a minimal composition of neutrons, protons, electrons and muons in the core and is obtained in a tabulated form from the CompOSE database.Finally, in most of our simulations,  (blue), on the edge between two patches (orange) and on the corner between three patches (green).The resolution of all the simulations is   = 25 and   = 30.Bottom: relative error of the "edge configuration" and the "corner configuration" with respect to the "centre configuration".
we use the simple, classical envelope model by Gudmundsson et al. (1983) (see Eq. D5 in Appendix D).This decision is based on its simplicity and its adequacy for the objectives of the current study.
We also discuss one case with a magnetized iron envelope model from Potekhin et al. (2015).
Since here we mostly focus on the thermal evolution, in most simulations we keep the magnetic field fixed.However, we will also present one simulation (up to 10 5 yr due to the substantially higher computational cost) with a fully coupled magnetothermal evolution, where also the Joule heating is accounted for.We refer to Dehman et al. (2023a) for the evolution of the magnetic field with simplified treatment of temperature adopted from Yakovlev et al. (2011) and to Dehman et al. (2023c) for the first fully coupled magnetothermal study with MATINS.We show below the results in terms of thermal evolution (i.e.cooling curve) and the corresponding pulsed profile expected from a given temperature map on the stellar surface.Notice that the total (photon) luminosity of the source at a given time that we report in the cooling curve has been calculated by integrating Eq.D4 over the entire stellar surface.For NS that are not too old, this luminosity is emitted mainly in the X-ray band, and thus this quantity can be compared with observations, if the distance and the interstellar absorption below ∼ 1 keV can be well estimated.Moreover, ideally, one can infer the temperature map on the surface of the star and the emission model (here assumed as a blackbody for simplicity) by a simultaneous fit of light curve in different energy bands.

Thermal evolution
We run different simulations, with and without magnetic field, with different EOSs and different field configurations.Our runs are summarized in Table 1.The first column denotes the name of the simulation, the second the gravitational mass of the star, and the third reports the EOS that we used.The remaining columns specify the configuration of the magnetic field.In MATINS the initial magnetic field configuration is set up by defining the multipolar weights of the expansion of the poloidal and toroidal scalar functions Φ and Ψ, -1,0,1 -2,-1,0,1,2 -7,-2,0,2,5 -6,-5,0,3,10 Table 1.Simulations setup.The first four columns report the name of the simulation, the mass of the star, the EOS, and the magnetic field intensity averaged over the crustal volume.The last four columns describe the initial magnetic field configuration.This is constructed by decomposing the magnetic field into poloidal and toroidal components, which are obtained from a poloidal and a toroidal scalar functions as described in Eq.F1 of the Appendix F. The poloidal and toroidal scalar functions are defined in turn by a spherical harmonic expansion, with the expansion coefficients as input parameters.The columns report the degree  and the order  of the non-vanishing spherical harmonic weights of the initial poloidal and toroidal field expansion.For more details on the initial magnetic field configuration we refer the reader to Appendix F or to Dehman et al. (2023a).The input parameters that define the configurations reported here are summarized in the supplementary Table F1.
respectively, as described in detail in Appendix B of Dehman et al. (2023a) (see their Eq.B3).The fourth column reports the average value of the magnetic field in the crust, the fifth and sixth (seventh and eight) columns report the degree  and the order  of the non-vanishing multipoles in the expansion of the poloidal function Φ (toroidal function Ψ), respectively.As in virtually any existing magneto-thermal simulation in literature, in all these simulations we impose potential boundary condition at the surface.The radial component of the magnetic field at the crust-core interface is kept at zero.We refer to Dehman et al. (2023a) for more details about boundary conditions, and to the Appendix F and Dehman et al. (2023a) for more details on the initial magnetic field configuration.
The first simulation that we show, the Sly4-M1.6-B14-L1,involves a magnetized star.The result of this run is presented in Figure 7.
Here the luminosity of the different neutrino emission processes is represented with colored lines.Continuous lines denote processes occurring in the core of the star and the dashed line denotes processes in the crust (the same process may involve both the crust and the core and in those cases, the same color appears in the plot with both a continuous and a dashed line).The continuous and the dashed black light denote the total neutrino luminosity for processes in the core and in the crust, respectively.Finally, the dotted black line represents the luminosity in photons.From Fig. 7, it is evident that the cooling process is primarily governed by neutrino emission until ∼ 25kyr, beyond which photon emission from the surface becomes the dominant cooling mechanism.
With respect to the neutrino emission, we can appreciate how the core provides the main contribution to the neutrino energy loss being 2-3 orders of magnitude larger than the crustal one, at all times.This reflects that most of the mass (typically ∼ 99%) is in the core, and that it is denser than the crust, although we need to acknowledge that the processes in the two regions are rather different, and as such they cannot be directly compared.In the core, the dominant mechanisms are represented by the Modified URCA processes and the neutron-neutron Bremsstrahlung until ∼ 100 yr, when the core becomes superfluid.After that, neutron Cooper pair formation is the main neutrino cooling mechanism.The onset of the superfluidity in the core is also followed by an increase in the slope in the photon luminosity.The crust, instead, is characterized by a much lower neutrino luminosity with respect to the core.Here the dominant mechanism is the neutrino-synchrotron emission, by virtue of the strong fields present in the crust with this setup.We stress that in the legend we indicate all the neutrino processes included in the simulation, even if some of them are not effectively active in this particular simulation and as such, they do not appear in the figure .For example, the neutrino pair-production is a process that becomes active at low densities and high temperature (it dominates at  ≲ 10 10 g/cm 3 and  > 10 9 K, see Fig. 3 of Potekhin et al. (2015)), but it becomes strongly suppressed, due to the reduced number of positrons, when the temperature drops below the Fermi value.The absence of direct URCA in the core can be instead ascribed to the low abundance of protons in the core that characterizes this setup.In Appendix G we show the activation of this important cooling process for a simulation employing a different EOS, which allows direct URCA above a given threshold mass.
These results do not depend on the 3D nature of the code and are in agreement with our current knowledge about the cooling of NSs.Nevertheless, they show how the interplay among different neutrino emission processes plays an important role in the modeling of the cooling curve of a NS.In this sense, to our knowledge, MATINS is the most advanced 3D code in terms of including all of these essential contributions.

Simulations with different magnetic configurations
Starting with the aforementioned setup for the EOS, the superfluidity, and the envelope model, we study the thermal evolution of a NS characterized by a mass of  = 1.40  ⊙ , varying the configuration of its magnetic field as described in Table 1.
Except for the dipolar cases Sly4-M1.6-B14-L1 and Sly4-M1.4-B14-L1,where we aimed at an axisymmetric poloidal configuration, all other non-axisymmetric configurations have been chosen arbitrarily.For a realistic initial magnetic field configuration inspired by  The same process (marked with a given color) may involve both the core and the crust.Black lines represent the total neutrino luminosity for processes involving the core (continuous line) and the crust (dashed line), respectively.The black-dots report the surface photon luminosity.It is worth noticing that while the legend includes all the possible neutrino processes included in the simulation, some of them are not effectively active in this particular simulation, as such, they do not appear in the figure.Moreover, in case magnetic fields are present in the core, additional processes can become relevant (Kantor & Gusakov 2021).
In Figure 8 we report the cooling curves of the case Sly4-M1.4-B14-L1(dipole), characterized by an axisymmetric poloidal field configuration, and the non-axisymmetric cases Sly4-M1.4-B14-L2(quadrupole), Sly4-M1.4-B14-L5 and Sly4-M1.4-B14-L10.As one can see, all the configurations have a similar cooling curve.This is due to the fact that the cooling processes are the same for all the configurations, while the Joule heating, which is determined by the (time-dependent) magnetic field intensity and geometry, is naturally absent since the field is constant in time.
A 3D rendering of the magnetic field lines and the temperature at the crust-envelope interface   (, ) evaluated after 2726 yr of evolution is reported in Figure 9 for the four different configurations just described 6 .From the figure, it is possible to appreciate the impact of the configuration of the magnetic field on the temperature 6 For the cases Sly4-M1.4-B14-L2(central left), Sly4-M1.4-B14-L5(central right) and Sly4-M1.4-B14-L10(right) we report here a case with enhanced angular resolution   = 41.
distribution.The cold areas on the map correspond to regions with a stronger magnetic field oriented perpendicular to the radial directions.In these areas, the magnetic field acts as an insulating barrier between the outer crust's surface and the hotter interior.As a result, while these regions lose thermal energy through neutrino and photon emission, there is no replenishment by heat diffusion, leading to the formation of cold spots.Vice versa, when the magnetic field is weaker, or it has a significant radial component, the outer regions of the crust are coupled and the heat lost can be replenished by diffusion, thus generating a hot area.Clearly, the complexity of the field reflects the complexity of the temperature map: a magnetic field characterized by small-scale structures leads to smaller hot-spots at the crust-envelope boundary.
In Fig. 8 we have also included (orange dashed line) a variation of the run SLy4-M.1.4-B14-L2,where the envelope described in Appendix B of Potekhin et al. (2015) is used instead of Gudmundsson's model.Similarly to Gudmundsson's model, the envelope in Potekhin et al. (2015) comprises heavy elements such as iron; however, unlike the former, the latter also incorporates the effect of the magnetic field in the modeling.As depicted in the figure, the cooling curves for the two models are comparable, with the magnetized envelope resulting in a slightly lower luminosity.This outcome will be discussed in more detail in Sec.4.3.Finally, we present the run SLy4-M.1.4-B14-L2-alt,which has the same multipole components of the run SLy4-M.1.4-B14-L2,but with a different field configuration, in which the contribution of the toroidal multipoles with  = 2 has been enhanced.This run is presented in two versions: a purely thermal one, where the evolution of the magnetic field is frozen like in the cases previously studied, and a magnetothermal one SLy4-M.1.4-B14-L2-alt-MT,where the magnetic field is evolved and Joule heating is correspondingly included.Since a full magneto-thermal simulation is considerably more expensive from a computational point of view than a simulation with a non evolving magnetic field, we restricted this run to a maximum time of 10 5 yr only.
In the magneto-thermal case, the effect of Joule heating can be appreciated in Figure 8.Here, the brown curves represent the magnetothermal run SLy4-M.1.4-B14-L2-alt-MT,which can be compared with the purple curve representing the same configuration with a non-evolving magnetic field.In this case, the heating produced by the dissipation of the currents has the effect of increasing the luminosity of the star by a factor ∼ 2 − 3 after the first ∼ 100 yr of evolution.
A 3D rendering of this run (with an enhanced angular resolution of   = 41) at time 2726 yr is reported in the bottom row of Fig. 10.The figure represents, from left to right, the temperature   , the magnetic field, and the electric currents.The upper row represents the same run with the non-evolving magnetic field (i.e.SLy4-M.1.4-B14-L2alt)at the same evolutionary time, with the same spacial orientation and with the same color range for the represented quantities, so that the images in the upper and bottom rows are directly comparable.As we can see, the temperature at the base of the envelope, in this case, is higher than in all the previous cases, due to the presence of the heating source.The other point to note is that, apart from an overall scaling, the angular distribution of   is radically different with respect to the equivalent case with the non-evolving magnetic field.It appears in fact that the two distributions are almost inverted, with the hottest regions of one roughly corresponding to the coldest regions of the other.The reason may rely on the fact that the coldest regions in the case with non-evolving magnetic field, which are those with an intense non-radial field that provides an insulating effect, are also the ones that in the magnetothermal case develop stronger currents and consequently more Joule heating.This effect is well known from previous 2D studies like Pons et al. (2009) and Viganò et al. (2013).This deep change in the   angular distribution is particularly important because it is expected to have a major impact on the pulsed profile observed from the source, as we will see in the next section.

Light curves
As we saw from the previous results, the magnetic-induced anisotropy in the thermal conductivity generates an inhomogeneous temperature distribution on the surface of the NS, presenting hot spots and colder regions.While the star undergoes rotation, the observer's view of the hot spots on the stellar surface may become periodically obstructed.This results in a periodic modulation of the stellar X-ray flux, which is eventually observable and is indeed observed in several sources (see e.g.Kaspi & Beloborodov 2017;Esposito et al. 2021).
Starting from our simulations, we compute the energy dependent flux (or, equivalently the phase-resolved spectrum) as detected by a distant observer.This is an important step because it allows us to link the results of our simulations with quantities that can be directly observed.
To this aim, we employ the ray-tracing code presented in Perna & Gotthelf (2008) and Viganò et al. (2014) and generalized to work with a non-axisymmetric temperature map.The ray-tracing code takes as input a temperature map from MATINS, and maps it via interpolation onto an internal, equally-spaced spherical coordinate grid.The temperature map is used to calculate the local emission spectrum on the surface of the NS.The phase dependent spectrum is computed (neglecting absorption from the interstellar medium) by integrating the specific intensity allover the stellar surface (Page 1995): Here  and  ∞ denote, respectively, the energy of photons at the star surface and the redshifted energy detected by an observer at infinity. ∞ is the radius of the star at infinity,  is the distance from the source, while the parameter  represents its rotational phase.We assume the local spectrum to be a blackbody   (), and also that the emission at the surface is isotropic.Both assumptions neglect all the possible spectral and anisotropic effects introduced by the passage of radiation trough the atmosphere of the NS, which further complicates the modelling (see e.g.Sec 2.3 of Özel 2013).While we make this assumption for the sake of simplicity in this first work and in order to highlight the dependence of the lightcurve on the temperature map, we note that the spectral and anisotropic distortions introduced by the atmosphere are second order effects.The angles  and  in Eq. ( 15) are the latitude and longitude angles on the star surface measured relatively to the axis aligned with the line-of-sight (LOS), while  ≡ sin , where  is the emission angle of the photon with respect to the normal to the stellar surface.At a given latitude  only the photons emitted with a given angle  will reach the observer.In our ray-tracing code, the relation between  and  is provided by the formula derived by Beloborodov (2002), which is valid for a Schwarzschild spacetime, where   and  are the Schwarzschild and the NS radius, respectively.It is important to acknowledge that the equation presented above serves as an approximation of the geodesic equation, and its accuracy may diminish for significant values of  (refer to the detailed discussion in La Placa et al. 2019).Nevertheless, within the context of our specific application, where the range of  is confined to  ≤ /2, this equation maintains an accuracy level better than 1%, which proves to be satisfactory for our intended purposes, while significantly reducing the computational cost of the code.The flux in Eq. ( 15) depends on the orientation of the temperature map with respect to the LOS, which also depends on the phase .This orientation (hereafter geometry) is defined by two angles: , that is the angle between the LOS and the rotation axis, and , that is the angle between the z-axis in MATINS and the rotation axis.The phase-resolved spectrum can be calculated once  and  are defined.
Here we present the phase-resolved spectrum for three selected cases: Sly4-M1.6-B14-L1,Sly4-M1.4-B14-L2 and Sly4-M1.4-B14-L10.For all the cases we chose the temperature map at  = 2726 yr, a geometry defined by  =  = /2, and a reference distance of  = 4 kpc.In all the ray-tracing simulations, the resolution of the grid is set to 150 points both in  and  directions.The temperature map obtained by MATINS is interpolated on this grid via a barycentric coordinate interpolation.
The first case that we present in Fig. 11 is the one referring to the simulation Sly4-M1.4-B14-L1.Here the top panel displays the phase-resolved spectrum, in an energy-phase plane.The color code represents the specific flux measured in photons/(s cm 2 keV), and the red continuous lines mark iso-flux contours.The bottom panel shows the flux integrated over the energy range 0.1 − 10 keV, measured in erg/(s cm 2 ).Since most of the emission occurs in this energy window, this flux coincides approximately with the bolometric thermal flux.The label here denotes the pulsed fraction (PF), that is the magnitude of the pulsed emission, computed as the difference between the maximum and the minimum of the flux divided by the sum of these.
This run is characterized by a dipolar configuration of a purely poloidal magnetic field.Here the magnetic dipolar moment is coincident with the -axis in MATINS.As such, the angle , in this case, is the inclination angle i.e. the angle between the rotation and the magnetic moment.The temperature map is characterized by a cold equatorial belt and hotter poles, which acts like two broad antipodal hot spots, centered on the magnetic axis.These features, along with the chosen geometry, are reflected on the pulsed profile, which shows two symmetric peaks within a single rotational phase.With the chosen geometry, which is the one that maximizes the PF, the LOS passes through the poles at  = 0, 0.5, where we have the maximum of the emission, while the minima at  = 0.25, 0.75 coincide with the LOS pointing at the equatorial cold belt.Note that in this specific case, due to the peak symmetry given by the magnetic configuration and inclination, the spin period inferred by the X-ray light curves would be half the real value.
In Figure 12, we display the same quantities as in Figure 11, but for the MATINS simulation Sly4-M1.4-B14-L2.This case is characterized by a more complex magnetic field, which results in a temperature map that is no longer axisymmetric like the previous case.As we can see from the figure, the phase-resolved spectrum and the pulsed profile show only a single slightly-asymmetric peak during the rotational phase.The PF is higher than the previous case, reaching a value of ∼ 15%.In the bottom panel of the figure we also show (orange curve) the pulsed profile for the variant of the same run, computed with the magnetized envelope of Potekhin et al. (2015).For this case we note that, although the overall flux is lower, the PF is slightly higher, and a second small peak appears in proximity of the flux minimum.
In order to illustrate the origin of the pulsed profiles, which are less straightforward to interpret than the previous one, we report in Fig. 13 the surface brightness of the star at three selected phases corresponding to the minimum (left), the maximum (center), and an intermediate value of the flux (right).The surface brightness has been computed by integrating the black-body specific intensity into the energy interval 0.1 − 10 keV at each point of the visible surface of the star and then projecting each point into the observer taking into account the propagation of photons in a Schwarzschild spacetime.In the top row we report the result for the standard run, while in the bottom row we report the result for the variant with the magnetized envelope.From the figure, we can see that the temperature map is characterized by several small spots of semi-aperture angle ∼ 10−15 • and a bigger, brighter spot of semi-aperture angle ∼ 35 − 45 • .The latter, due to its enhanced size and brightness, dominates the emission producing the peak at  ∼ 0.5, when it is oriented face-on with respect to the observer.The minimum of the emission in turn occurs at  ∼ 0, when the big spot is completely out of sight.Comparing the two runs characterized by different envelopes we can appreciate how, in the case of the magnetized envelope, the gradients in surface brightness (i.e. in the surface temperature) are larger with respect to the Gudmundsson envelope case; this is in agreement with the slightly higher PF observed in this case.Where the field is tangential to the surface, the magnetized envelope model provides an enhanced screening; viceversa when the field is radial, the surface temperature is higher with respect to the Gudmundsson case, as also reported in Dehman et al. (2023b).This is the case of the two quasi-antipodal hot spots visible in the bottom row of Fig. 13, which are hotter with respect to the equivalent regions in the top row.It is worth noticing that the magnetized envelope enhances the temperature only for these two spots and not for the other existing hot regions.The reason is that these two spots are characterized by a field with a substantial radial component, and as such they are poorly screened by the magnetic envelope.On the other hand, the other hot areas have a high temperature not because of the field orientation but due to its weak local intensity; the field however is tangential to the surface in those locations, such that they are subjected to a higher screening when a magnetized envelope is used.
Note that, unlike the dipolar case, in this case the -axis does not pass through the hottest region.As such it is not guaranteed that this geometry ( =  = 90 • ) is the one that maximizes the PF.Nevertheless, we can see from Fig. 13 that since the widest and hottest region passes close to the center of the figure at the maximum of the pulse profile, while it completely disappears at the its minimum, it is reasonable to expect that this configuration may not be too far from the one with highest PF.
The final case, depicted in Fig. 14, corresponds to the model Sly4-M1.4-B14-L10.Similarly to the quadrupolar case, this particular scenario exhibits one peak that appears within a single rotational phase.However, in this case the shape of the profile is more complex, reflecting the more complex shape of the hot region, which can also be appreciated in the top right panel of Fig. 9. Around  ∼ 0.8 we notice a step-like shape which has a shift in phase of Δ ∼ 0.5, with respect to the peak.This particular profile suggests the presence of two hot regions positioned at almost antipodal positions and characterized by different brightness levels.
Comparing the PFs of the three cases presented here, it is worth mentioning that the Sly4-M1.4-B14-L10configuration has the lowest PF.This is due to the fact that this temperature map is characterized by a high number of hot regions, so that while some of them are occulted by the star rotation some other are revealed, thus reducing the PF.For this reason, even if also in this case, like for the quadrupole one, the choice of  =  = 90 • may not necessarily be the one maximizing the PF, we argue that the maximum PF is probably not too far from the value reported here.
The second case with the lowest PF, comparable to Sly4-M1.4-B14-L10, is the dipole.In this case the low value of the PF can be ascribed to the wideness of the polar hot regions which can never be completely occulted, unlike for the quadrupolar case.
To conclude, we consider the run Sly4-M1.4-B14-L10,with the inclusion of the full magnetothermal evolution.Also for this run, we consider a temperature map at i.e.  = 2726 yr, which gives enough time for the magnetic field to experience some Hall and Ohmic evolution.The result is presented in Fig. 15.The top and the middle panel represent the phase-resolved spectrum expressed in photon counts for the case with non-evolving (top) or evolving (middle) magnetic field.The bottom panel represents the pulsed profile of the two cases.Hereafter we will refer to the first case as thermal case, and to the second case as magnetothermal case.In this plot, we can notice two interesting features: first of all, as expected, the Ohmic dissipation contributes to increasing the surface temperature, and consequently the flux, which increases by almost a factor 5, while the peak of the spectrum moves to higher energies.Secondly, there is a phase-shift of Δ ∼ 0.5 of the profile maximum between the two cases, such that the peak of the pulsed profile in the case with the non-evolving field coincides with the minimum of the magnetothermal profile.This shift denotes a switch in the hot regions with and without the Ohmic dissipation, as already noted in the previous section.Finally, we note that the magnetothermal case is characterized by a PF of ∼ 30 %, which is a value higher than the PFs of all the other cases studied here, where the field is not evolving.

SUMMARY
The present work, together with Dehman et al. (2023a) introduces MATINS (MAgneto-Thermal evolution of Isolated Neutron Star), a three-dimensional code designed to investigate the secular magnetothermal behavior of a NS, featuring a self-consistent coupling between thermal and magnetic processes, a realistic star structure based on different possible cold dense matter EOS, and detailed local microphysics.While the previous study by Dehman et al. (2023a) described how the code solves the magnetic evolution equation with a particular emphasis on the cubed-sphere grid, our focus here is to illustrate the formalism, testing and some applications of the thermal evolution part.The results from the fully coupled induction and heat diffusion equations have already been presented in Dehman et al. (2023c), which considered a realistic magnetic field configuration derived from a proto-NS dynamo simulation (Reboul-Salze et al. 2021), along with the most updated envelope model (Potekhin et al. 2015).
The MATINS code utilizes a finite volume scheme to solve the thermal evolution (energy conservation) equation in the cubed sphere grid coordinate system (Ronchi et al. 1996).The equation is solved via an implicit scheme, which is optimal to treat the stiffness of the neutrino emissivity, and allows us to use larger time steps, due to the unconditional stability of the method.
We included in our code a numerical implementation of the microphysics that exploits the public code by A. Potekhin7 , which provides a detailed description of microphysical parameters in the outer and inner crust, and in the core of the NS, along with the inclusion of all the relevant neutrino emission processes in the crust and in the core.The structure of the NS is obtained by solving the TOV equation at zero temperature.This requires the choice of an EOS that can be selected from the public online database CompOSE8 .
We applied MATINS to various scenarios, including different EOSs and masses, and also to a case with a fully coupled magnetothermal evolution.We observed the impact of different magnetic field configurations on the distribution of temperature inside the star and in particular on the interface between the envelope and the crust, which is directly related to the temperature map at the surface of the star.We appreciated how the increasing complexity of the field translates into the complexity of the temperature map, where field configurations with small-scale structures result in the formation of smaller hot spots.Moreover, we coupled the output of MATINS to a generalrelativistic ray-tracing code to calculate the phase-dependent spectrum observed by a distant observer from a thermally emitting NS.This code is essential to link the results provided by MATINS in term of surface temperature maps with the observations.In particular, we showed how the simulations discussed in this paper produce different phase-resolved spectra, pulsed profiles, and pulsed fractions for a given orientation of the LOS with respect to the temperature map and star's rotation axis (geometry).The few cases presented in the paper represent just a proof of concept of the synergy between our ray-tracing code and MATINS, and a more systematic study considering a wider sample of simulations and different geometries will be left for future work.
In conclusion, by considering the interplay between thermal and magnetic processes, MATINS provides a more comprehensive and realistic understanding of the complex thermal histories of NS.The versatility of the code in handling various scenarios, together with its successful application in testing and simulations, offers promising potential for advancing our knowledge of NS characteristics and their evolutionary history.In particular, it is an essential tool to simulate and fit simultaneously the timing properties (controlled by the dipolar surface field, Viganò et al. 2013), the thermal luminosity and the Xray light curves.
where Δ =  +1 −   is the discretized timestep leading from step  to  + 1, and the specific heat and neutrino emissivities are evaluated locally.Note that the redshift factors only depend on the radial direction.Here Φ , , , ( 2  F , , , ) represents the net flux across the volume and it is evaluated along the cell surface  , , , , as follows: where the half-integer values of the indices denotes quantities evaluated at the cell interfaces.Once discretized, it is possible to notice that the net flux can be expressed as a linear combination of the temperature within the cell and the temperature in the first neighbouring cells, such that we can express it in matrix form: Eq. (B1) is solved using an implicit backward Euler method, which is particularly suitable to handle the stiff terms appearing in the neutrino cooling and to increase the timestep preserving the stability of the solution.This implies that the flux is written in terms of the updated temperatures, namely those calculated at the timestep  + 1: where the source   contains terms of the discretized equation which depend only on the local old temperatures, and whose explicit form, along with that of the matrix m   , is reported in Appendix E. This corresponds to expressing our differential equation as a system of 6  (  + 2) 2 algebraic equations, whose unknown values are the temperatures in each cell at the timestep  + 1.In our code we solve the system with the aid of the public library LAPACK (Anderson et al. 1999).
The right hand side of Eq. (B5) corresponds to the old temperatures plus the source term (as in the right-hand-side of Eq.B4), which includes the Joule heating and the neutrino emissivity.The latter contribution has a steep dependence with temperature, making the source term a stiff term in our equation.In order to handle it, it is convenient to perform a linearization, writing:

APPENDIX C: GHOST CELLS AND PATCH EDGES
As we mentioned in Sec.2.3 we introduced a layer of ghost cells at each patch interface.The temperature of these extra cells is used to calculate the heat flux through the cells next to the patch edges; consequently they must enter in Eq.B5.While the other temperatures are determined by a physical equation (Eq.6), the temperatures in the ghost cells are determined by a linear interpolation with the temperature of two cells in the neighbouring patch.This results in an equation that is actually a constrain and writes as: Here  +1 ghost,j is the temperature of a ghost cell a given patch,  +1  and  +1  ′ are the temperatures of two non-ghost cell in the neighbouring patch.The index  is the index of the coordinate parallel to the edge, while  ′ =  +1 for  ≤   /2 and  ′ =  −1 for  >   /2.  are the weights defined as: where    and    are the ghost cell coordinates, defined in Sec.2.6 of Dehman et al. (2023a).Eq.C1 basically states that the temperature of the ghost cell of a given patch is the linear combination of the temperature of the cell at the edge of the neighbouring patch with the same parallel (to the edge) coordinate and the temperature of the cell next to it towards the patch center (to have an insight the reader can refer to Fig. 2 of Dehman et al. (2023a)).Finally, the term   represents the weights of the linear interpolation and it belong to the range [0, 1], where   = 0 at the center of the edge (  =   /2,   /2 + 1) and   = 1 at the edge extremes (  = 1,   ).In this way, at the center of the edge  +1 ℎ ,  takes the contribution only from a single cell, while at the extremes it correspond to the numerical average of the two cells.
A separate treatment is necessary for the cells with ,  = 0,   +1, which corresponds to the ghost corner of the patch.These cells are special because they are ghost cells with respect to two different edges, namely, in our approach, their temperature must be coupled to temperature of two different patches.Similarly to Eq. C1 the temperature of the corner ghost cell is written as: where  +1

APPENDIX D: BOUNDARY CONDITIONS
Our setup requires the definition of boundary conditions on the innermost and the outermost radial layers, represented by the interfaces with the core and the envelope, respectively.The core is described by a single radial layer because, due to its high conductivity, it becomes perfectly isothermal (precisely,    ( ) is constant) after at most only a couple of centuries, a timescale that is irrelevant for our purposes and without observational constraints (see e.g.Viganò et al. 2021).For this reason, in our code, we consider a single redshifted temperature for the core, which evolves according to the ratio between the volume-integrated sources and heat capacity, and accounting for the core-crust flux: ) where the neutrino emissivity and the specific heat are integrated over the entire core volume, exploiting the fact that they depend only on the radial coordinate.The flux term Φ  is the net flux at the core-crust interface, written as: which is Eq.B2 integrated along the whole crust-core interface, where we consider only the radial contribution of the flux at this boundary.
Eq. D1 is solved by an implicit backward Euler method.However, it is worth noticing that the flux Φ  is treated explicitly.The reason is that in this term the flux does not depend only on the temperature of the core, but also on the temperature of the cells in the innermost crustal layer ( = 2).We observed empirically by refining the temporal resolution that this approach, alas not fully self-consistent, has a negligible effect on the solution.
From an operational point of view we treat the cells in the innermost layer as ghost cells, in the sense that like the ghost cells at the border of the patches (see Appendix C), their temperature is determined by a constrain, which is written as: whre  +1 core is determined by the solution of Eq.D1.The external boundary condition consists in the addition of a source term, describing the emission from the stellar surface.Such emission is assumed for simplicity to be a blackbody with surface temperature   and is given by the Stefan-Boltzmann law (written in a discretized form) where  is the luminosity emitted by the radiating area  at temperature   , and   is the Stefan-Boltzmann constant.
It is worth noticing that here the surface temperature   is not the temperature in the outermost layer of the computational domain, corresponding to the crust-envelope interface and denoted by   .  is instead the temperature at the top of the envelope.The envelope is not part of our computational domain, since in this region the gradients of temperature and density are so steep and the timescales so short that including this zone makes the thermal evolution computationally prohibitive to tackle.Instead, in MATINS we utilize envelope models, which are analytic effective   =   (  , B) relations, obtained by fitting stationary solutions of heat diffusion in the envelope models computed with different   and B. Most of the simulations presented in this paper utilize the envelope model provided by Gudmundsson et al. (1983): where  14 is the surface gravity in units of 10 14 g cm/s 2 .We chose this model because it is particularly simple, not inclusive of the dependence on the magnetic field, which makes it straightforward to interpret the cooling curve in terms of the   profile.Nevertheless, MATINS accounts also for several other models, including magnetized envelopes with heavy element or light element composition.In addition to Gudmundsson model, we provide also a run characterized by the magnetized iron envelope presented in Potekhin et al. (2015).
To see the impact of the envelope model on the cooling curve we refer to Dehman et al. (2023b).
The term in Eq.D4 is treated implicitly similarly to the neutrino cooling: the source term, evaluated at  +1 , is linearized and the terms dependent on   are included in the source vector   , while those dependent on  +1 are included in the diagonal elements of the matrix m It is worth noticing that while such implicit treatment is essential for the neutrino emissivity, due to the strong stiffness of this term, for the photon luminosity it just represent a second order correction that gives little difference compared to an explicit treatment.

APPENDIX E: MATRIX ELEMENTS
In this Appendix we write extensively the matrix elements    and    and we show how to derive them.First of all we start from Eq. 11 and we explicit all the three term on the right hand side of the equation.We have first the temperature gradient, which according to Eq. (A11) writes in cubed sphere grid as: The second term of interest is the scalar product b • ∇ ∇ ∇ T, which writes as: Finally, we report the three component of the vector product b×∇ ∇ ∇ T in the last term: It is useful to write the three terms in eq.11 as follows: where the terms , , , , , , , , , , , , ,  are the following functions: It is straightforward to show that the matrix elements    in Eq. 12 are: Concerning the matrix elements    in Eq. (B5), the first step is to consider the net flux in Eq. (B2) and make explicit its dependence on the temperature.In its discretized form this equation writes as: where the six terms on the right hand side are the contributions from the six faces of the cell.We can now write explicitly the discretized fluxes in the previous equation exploiting Eq. ( 12).We evaluate the temperature derivatives as cell centered difference between the values at the first neighbours.They are located either at the centers (for the derivatives associated to the direction normal to the surface, namely the diagonal terms of the conductivity) or at the middle of the edges of the 3D cells (associated to the transverse conductivity, namely the off-diagonal terms).The latter temperatures are obtained as the average between the values at the centers of the four closest cells.

APPENDIX F: INITIAL MAGNETIC FIELD AND INPUT PARAMETERS OF OUR SIMULATIONS.
In this section, we detail the formalism used to define the magnetic field initial configuration and, for the sake of reproducibility, the input parameters that define the configurations presented in Table 1 in Section 4. Generally, a magnetic field B can be decomposed into a poloidal component B pol and a toroidal component B tor , which can be expressed via two scalar functions Φ(, x) and Ψ(, x) as follows: where k is an arbitrary vector, which in our coordinate system is conveniently identified with the radial versor k = r.To define the initial magnetic field configuration it is convenient, following Geppert & Wiebicke (1991) where  = 1, .. max is the degree and  = −, ...,  is the order of the multipole.
The initial configuration of the magnetic field is fixed by choosing a set of spherical harmonics.In this way, for example, we can construct a general dipolar field by choosing Φ  , Ψ  ≠ 0 for  = 1 and Φ  = Ψ  = 0 for  ≠ 1.
Regarding the radial dependence of the coefficients Φ  and Ψ  , we impose the radial profile of the dipolar poiloidal scalar function Φ =1, () as in eq. ( 8) of Aguilera et al. (2008) (eq.B9 of (Dehman et al. 2023a)): and  is a parameter related to the curvature of the field, calculated for a given stellar radius .This choice of radial dependence allows for a smooth match with the external potential boundary condition.
On the other hand, the toroidal field and the  > 1 poloidal field multiples are confined inside the crust of the star, with the following radial dependence: where   and   are normalization input parameters that are set at the beginning of the simulation.
Finally, the poloidal component of the field is normalized to the maximum (absolute) value of the radial field at the surface and rescaled by the input parameter  pol,init , while the toroidal component is normalized to the average root mean square within the volume and rescaled by the input parameter  tor,init .In this way the input parameters  pol,init ,  tor,init ,   and   determine the configuration and the magnitude of the magnetic field in MATINS.In Table F1 we report the input parameters defining the configuration of our simulations.The names of the simulations are the same as the ones reported in Table 1 .

APPENDIX G: SIMULATION WITH A DIFFERENT EQUATION OF STATE
In this Appendix we present two simulations of a non-magnetized star using the BSk24 EOS (Goriely et al. 2013): BSk24-M1.4-B0and BSk24-M1.8-B0,characterized by  = 1.47  ⊙ and  = 1.87  ⊙ , respectively.Since these simulations do not present any magnetic field, we adopted for these runs a reduced angular resolution   = 7, while   = 30.The aim of this run is to show how with an appropriate EOS increasing the mass above a threshold value allows for the activation of Direct URCA processes, which leads to more efficient cooling.This result is shown in Figure G1.In the top panel, the secular evolution of photon luminosity is reported.The blue line represents the case BSk24-M1.4-B0,while the orange line represents the case BSk24-M1.8-B0.The bottom left panel represents the neutrino luminosity for the case BSk24-M1.4-B0,while the bottom right panel is for the case BSk24-M1.8-B0.Solid curves represent processes located in the core, while dashed lines represent those ones occurring in the crust of the star (note that the same process can occur both in the core and in the crust, such as n-n Bremsstrahlung).From Fig. G1 we can appreciate how the direct URCA processes are the main contributor to the neutrino luminosity in the BSk24-M1.8-B0run, while it is completely absent both in the BSk24-M1.4-B0and in the Sly4-M1.6-B14-L1case in Fig. 7 that will be discussed next.It is worth noting that in this case most of the energy is radiated by the process within the first years of the life of the NS.In contrast, in BSk24-M1.4-B0,the dominant mechanism is the modified URCA, which radiates a much smaller amount of energy.The enhanced efficiency of the direct URCA with respect to the modified URCA reflects in the photon luminosity, which drops by more than 2 orders of magnitude with respect to the lower mass case before the first 100 yr of evolution, in agreement with our current knowledge about NS cooling (Page & Applegate 1992).

Figure 1 .
Figure 1.Visualization of the cubed sphere grid at a given radius.The edges between the patches have been highlighted with a bold line.Within each patch, thin lines mark the  and  directions.

Figure 2 .
Figure2.Relative L2 deviation between the analytic and numerical solution as a function of time (measured in diffusion timescale units).The circles represent simulations with the magnetic field oriented along the z-axis with low (  = 7,   = 20; blue), medium (  = 14,   = 20; orange), and high (  = 28,   = 20; green) resolution.Red crosses and purple diamonds are simulations at medium resolution with different magnetic field orientations.

Figure 3 .
Figure 3. Meridional cut of the internal temperature profile (in dimensionless units) of the benchmark test.In this test, the microphysics is uniform, with dimensionless value of  v =  ⊥ = 1, the magnetic field is uniform and oriented along the -axis, and with    0 = 10.The initial temperature profile has a cylindrical symmetry around the magnetic field axis and evolves in time according to Eq. 13.The figure shows the run with high resolution (  = 20,   = 28) at three different times ( = 0 diff ,  = 2 diff and  = 4 diff ).

Figure 4 .
Figure 4. Equatorial cut of the internal temperature profile (in dimensionless units) of the benchmark test with the field oriented along the z-axis and high resolution (  = 20,   = 28) at time  = 4 diff .

Figure 5 .
Figure 5. Map of the relative error in the benchmark test on a spherical shell at  = 0.25( out +  in ).The patches' edges are underlined with a red curve.

Figure 6 .
Figure6.Top: Cooling curve for three different configurations of the Sly4-M1.4-B14-L1model, where the magnetic axis falls on the center of the patch (blue), on the edge between two patches (orange) and on the corner between three patches (green).The resolution of all the simulations is   = 25 and   = 30.Bottom: relative error of the "edge configuration" and the "corner configuration" with respect to the "centre configuration".

Figure 7 .
Figure7.Evolution of neutrino and photon luminosities of all the different emission processes over the NS history for the simulation Sly4-M1.6-B14-L1.This run, which uses the SLy4 EOS, is characterized by a mass of  = 1.6  ⊙ and a dipolar poloidal magnetic field with a crust-averaged value of  avg = 6 × 10 14 G.The chosen resolution is   = 25,   = 30.Solid-colored lines represent different neutrino emission processes occurring in the core.Dashed lines represent different neutrino emission processes occurring in the crust.The same process (marked with a given color) may involve both the core and the crust.Black lines represent the total neutrino luminosity for processes involving the core (continuous line) and the crust (dashed line), respectively.The black-dots report the surface photon luminosity.It is worth noticing that while the legend includes all the possible neutrino processes included in the simulation, some of them are not effectively active in this particular simulation, as such, they do not appear in the figure.Moreover, in case magnetic fields are present in the core, additional processes can become relevant(Kantor & Gusakov 2021).

Figure 10 .
Figure 10.3D rendering at 2726 yr of the model SLy4-M.1.4-B14-L2-alt.Top: Simulation with non-evolving magnetic field.The two panels show, from left to right, the temperature and the magnetic field lines.Bottom: full magnetothermal simulation, which includes also the evolution of the magnetic field and Joule heating.The three panels show, from left to right, the temperature, the magnetic field lines, and the electric currents lines (related to the Ohmic heating).All the simulations adopt a mass of  = 1.4  ⊙ , the EOS SLy4, and a resolution of   = 30 and   = 31.

Figure 13 .
Figure 13.Surface brightness of the star projected on the observer plane at three different rotational phases.Top: simulation with Gudmundsson envelope.Bottom: simulation with the iron magnetized envelope.The rotational axis is horizontal in each image.The images refer to the three (per simulation) rotational phases highlighted in the right panel of Fig. 12.The color of the frame around each image matches the color of the dot in the previous figure to which the image refers.Note that the phases in the top row and in the bottom row are close but not exactly equal (refer to Fig. 12).
the non-ghost cells at the corner of the two patches adjacent to the original patch at the given corner.For example, if we are considering the ghost corner ,  =   + 1 of patch ,  +1 1 and  +1 2 are the temperatures in the cells  =   ,  = 1 of patch  and  = 1,  =   of patch  , respectively.