An Algorithm for the Simulations of the Magnetized Neutron Star Cooling

The model and algorithm for the cooling of the magnetized neutron stars are presented. The cooling evolution described by system of parabolic partial differential equations with non-linear coefficients is solved using Alternating Direction Implicit method. The difference scheme and the preliminary results of simulations are presented.


Introduction
: Observational data of age and temperature of compact stars. Squares -AXPs and SGRs neutron stars with magnetic field from 10 14 G up to 10 15 G. Triangles -radioquiet stars with field around 10 13 G. Circlesradio-pulsars with field ≤ 10 12 G.
Modelling of magnetized neutron stars cooling process is an interesting subject for the investigations of the internal structure of such objects. The main aspect is the investigation of the nuclear matter properties in extremely high densities. The magnetic field of compact object when it is more than 10 14 G has a significant influence on the heat transport inside the star and can give observational effects on the surface temperatures. The model of cooling evolution is given by a system of parabolic partial differential equations with non-linear temperature depending thermal coefficients and cooling (via neutrinos) and heating (due to magnet field decay) sources.
Due to the axial symmetry of the field and spherical symmetry of the matter distribution inside, the problem can be realized with the choice of spherical coordinates in spatial 2D. For the solution we choose the Alternating Direction Implicit (ADI) method [1,2]. In the difference scheme we use non-constant spatial steps and self-correcting a e-mail: hovikgrigorian@gmail.com arXiv:1511.05380v1 [physics.comp-ph] 17 Nov 2015 EPJ Web of Conferences time step. We investigate the special boundary conditions in the center and on the surface of the star configuration. At the current stage of the study we are going to present some preliminary results of the cooling simulations to demonstrate the efficiency of our 2D cooling algorithm.
The nowadays observational data of the magnetized neutron star (Fig. 1, the table of data with the corresponding references see in [3]) could not be explained by 1D cooling simulations neglecting magnetic field inside. Inclusion of the magnetic field is necessary but even not sufficient yet [3,4]. In the current work we focus on the algorithm of 2D simulation of cooling of magnetized neutron star.

Equations for Thermal Evolution of the Magnetized Compact Star
The neutron star is bounded by gravitational field, which in the approximation of slow rotation of the stars can be described by a metric tensor of the space-time manifold [5] ds 2 = −e 2Φ dt 2 + e 2Λ dr 2 + r 2 dΩ 2 . (1) The compact star configuration is constructed with the use of equation of state (EoS) of the stellar matter based on knowledge of nuclear matter [6]. The coefficients of the metric tensor as well as characteristics of the matter distribution are self-consistent solutions of the Einstein equation (specially for this case TOV [5,6]) where the temperature effect on internal structure could be neglected because of high density inside the star.
On the other hand the thermal evolution of the star can be described with the use of energy balance equations, which has parabolic equation form: where c v is the specific heat per unit volume and Q is the energy loss/gain by neutrino emission, Joule heating, accretion heating, etc. The vector F is corresponding to the heat flux. The thermal properties of the stellar matter are investigated in frame of different choices of nuclear matter cooling scenarios [7][8][9][10][11].
The flux is created due to the temperature gradient and tends to equilibrate the temperature Hereκ is the thermal conductivity, which is a tensor in the case of enough strong magnetic field of the star.
Hereafter we will use the notationsT ≡ e Φ T andF ≡ e 2Φ F. The flux components (in spherical coordinates) in the case of dipole magnetic field are: when the φ-component is zero due to the axial symmetry. Influence of the magnetic field, which effects the conductivity, leads to anisotropy along and orthogonal to the field. The relation between conductivity components is defined in terms of the magnetization parameter ω B τ: where τ is the particle relaxation time [12], and ω B is the classical gyrofrequency corresponding to the magnetic field strength B: Mathematical Modeling and Computational Physics 2015 m is the effective mass and e is a charge of a particle involved with magnetic field, here c is the speed of light.
With the use of unit vector b of magnetic dipole the flux can be expressed in the form: When the magnetic field is parallel to the axis z one has b = (cos θ, − sin θ, 0). So for the components of the conductivity we have: For simplicity we will use the following notations: ξ ≡ cos(θ),Q ≡ e 2Φ Q.
The introduction of the new variables ξ and the accumulated mass m = 4π ρr 2 dr (here ρ is the energy density of the matter) instead of θ and r leads the thermal evolution equation (2) to the following form The sought-for function is u = log(T ) and the coefficients are

2D Scheme for Mixed Implicit and Explicit Methods
For the numerical calculation we use the most convenient from stability point of view ADI method [1,2] and discretize the direction m taking into account the energy distribution ρ(r). This function is the solution of a single compact star for given central density or given total mass. So, corresponding to each given distance from the center r j we have points of grid m j . The mass for the given j is the same for all angles, therefore ξ could be considered as an independent coordinate. The steps of discretization are m j±1/2 = ± m j±1 − m j and m j = (1/2) ∆m j+1/2 + ∆m j−1/2 = (1/2) m j+1 − m j−1 . Similar notations we use for the angle ξ as well as for any function F(r, θ): So for the first order partial derivatives one has and for the second order onces we introduced R and S operators correspondingly for the same and mixed variables: Finally after discretization of the time with the steps ∆t n the differential equation of the problem can be written as: The introduced parameter σ stands for ADI method: when σ = 1 we apply implicit method in m direction and when σ = 0 in ξ direction.
The tri-diagonal elements and right hand side of the algebraic equations can be written in the following form where we use the following notations: is the main diagonal vector. For solving this tri-diagonal algebraic system the Thomas algorithm is used [13,14]. In the formula (6) the K (1) i, j and K (2) i, j are the centred mean values of To complete the tri-diagonal algebraic system we put the set of boundary conditions to the physical requirements, namely in the center and on the surface of the neutron star [15], as follows. All fluxes of heat from all directions have to vanish at the center point, which leads to X [i][0] = 0. Therefore, the temperature at the center point could be set as mean value of temperatures in very close area of the center.
At the surface we assume that the flux of the energy towards the radial direction has to be set equal to the flux of energy of the photon emission: where J is the maximal value of index j, σ is the Boltzmann constant and T s is the surface temperature, which is connected to the T m -mantle temperature under the core (see Fig. 2): For details one can see [7].
From the symmetry reasons the tangential fluxes on the polar and equatorial planes also are zeros L ξ (i = 0, j) = L ξ (i = I, j) = 0:

Some Features of Simulations and Conclusions
To demonstrate the work of our algorithm we provide simulations with model coefficients, which are taking into the account the main features of physical properties of thermal coefficients. Particularly, we use cooling mechanism of neutrino emission -so called modified URCA process (see for details [17]). We also vary the magnetic field, however we don't include the Joule heating from the decay of magnetic field. As it can be seen from the Fig. 3a due to absence of additional anisotropic heating sources the temperature distribution in all directions becomes isotropic after 50 years of evolution. And we can notice as well that this property is independent from the value of initial anisotropic distribution of the temperature.
One can notice that if temperature in some direction initially less then the temperature at the beginning of the isotropic era it heats up during the evolution (see red lines on Fig. 3a). It means that heat conducting process dominates over neutrino cooling.
The next, we note that this period 50−100 years are short enough to be able to compare this results with observations, because all data we have as it shown in Fig. 1 are for objects older than at least 500 years. On the other hand one can see that the temperature values at the beginning of the isotropic era are around the same what are the observational temperatures of already cold magnetars.
It means that long stay of the heat inside the star for such longer time could not be just the effect of anisotropy due to the magnetization. It is even hard to argue that the observed temperatures of red point in Fig. 1 are estimated only for small area around the poles of the stars.
On the Fig. 3b we demonstrate the automatic choice of the time-step of the numerical algorithm. The top line on this figure is the predicted time-step, the bottom curve is the modified time-step to provide stability of the algorithm. The automatic choice of time-step is needed mainly due to nontrivial mass distribution inside the star.