# Standard cooling of rapidly rotating isolated neutron stars in 2D

Mikhail Beznogov, Jérôme Novak, Dany Page, Adriana Raduta

► To cite this version:

Mikhail Beznogov, Jérôme Novak, Dany Page, Adriana Raduta. Standard cooling of rapidly rotating isolated neutron stars in 2D. The Astrophysical Journal, 2023, 942 (2), pp.72. 10.3847/1538-4357/ac9eb7. hal-03819896

HAL Id: hal-03819896

https://hal.science/hal-03819896

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Standard Cooling of Rapidly Rotating Isolated Neutron Stars in 2D

Mikhail V. Beznogov<fnref n="1" />,<fnref n="2" /> ID, Jérôme Novak<fnref n="3" /> ID, Dany Page<fnref n="2" /> ID, and Adriana R. Raduta<fnref n="1" /> ID

Received 2022 June 9; revised 2022 October 28; accepted 2022 October 29; published 2023 January 12

## Abstract

We study the long-term thermal evolution of axisymmetric rotating neutron stars in full general relativity. To this aim, we develop NSCool 2D Rot, a major upgrade to the 1D neutron stars thermal evolution code NSCool by Page. As a first application of our new code, we address the standard cooling of isolated neutron stars with rotation frequencies up to the mass shedding limit. We investigate the effects of the equation of state (EOS) by considering different combinations of core and crust EOSs. The results indicate complex time-dependent evolution of temperature distribution throughout the whole volume of the star, and in particular, in the crust. We show that most of that complexity can be attributed to the formation of a heat blob in the crust and the latitude dependence of the heat diffusion timescale through the crust.

*Unified Astronomy Thesaurus concepts: Neutron stars (1108); Neutron star cores (1107); X-ray point sources (1270)*

## 1. Introduction

Neutron stars (NSs) are born in core-collapse supernova events (Baade & Zwicky 1934) or in accretion-induced collapses of white dwarfs (Canal & Schatzman 1976) and start their lives as proto-neutron stars (Burrows & Lattimer 1986). Alternatively, an NS can be reborn in the merger of a binary NSs system provided that it survives the event (i.e., does not collapse into a black hole, see Kluźniak & Ruderman 1998). In a core-collapse supernova, the proto-neutron star phase starts after the core bounce and ends a few tens of seconds later when a proto-neutron star becomes an ordinary NS that is transparent to neutrinos. By this time the deleptonization process is already finished, and the entropy has decreased significantly (see, e.g., Pons et al. 1999; Roberts 2012; Nagakura et al. 2020; Pascal et al. 2022) and the star has (almost) reached its final size and composition. During the next few hours to a day, it will cool down to temperatures well below $10^{10}$ K, shrink a little bit more and form a crust (Beznogov et al. 2020). After the end of this so-called neo-neutron star phase, an NS can be considered as composed of cold-catalyzed matter.

NSs provide a unique opportunity to study the properties of cold superdense matter (e.g., Haensel et al. 2007), its equation of state (EOS) being the most wanted characteristic. The mystery of the EOS involves two facets: pressure (P) versus density ($\rho$) dependence and chemical composition. The former can be inferred by different techniques such as, e.g., independent measurements of NS masses and radii (e.g., Miller et al. 2019, 2021; Riley et al. 2019, 2021) and gravitational wave (GW) observations of binary NS mergers (e.g., Abbott et al. 2017, 2019). The question of composition, however, is more elusive and the study of NS thermal evolution is one of the most promising tools to tackle it.

To investigate thermal evolution one needs to solve two types of (in the general case) coupled equations: mechanical

structure equations and thermal transport and balance equations. However, in the particular case of long-term evolution, the situation becomes significantly simplified. As mentioned earlier, within a few hours to a day after the birth of an NS its temperature drops well below $10^{10}$ K. At such temperatures, thermal effects on NS composition and mechanical structure become negligible everywhere, except for the outermost layers of the star, which constitute the heat-blanketing envelope and have a negligible impact on the mechanical properties of the rest of the star (Baym 1981). Thus, the chemical composition and equilibrium configurations may be assessed by employing a barotropic EOS of cold baryonic matter in $\beta$ equilibrium everywhere, except for the outer envelope. This allows one to decouple mechanical structure and thermal evolution equations. The former then has to be calculated only once while the latter are solved repeatedly as time progresses (e.g., Baym 1981; Haensel et al. 2007; Potekhin et al. 2015). In this situation, the surface temperature, which is the only really measurable quantity, is determined by the heat-blanketing envelope via the so-called “$T_s-T_b$” relation (e.g., Gudmundsson et al. 1983; Beznogov et al. 2021). The validity of such an approach to the long-term cooling of NSs was explicitly demonstrated by Dong et al. (2018) for the spherically symmetric (1D) case by comparing the cooling of NSs with dynamic (i.e., non-fixed) and static (i.e., fixed) mechanical structures. For more subtle details of this widely used method see Appendix A. Throughout this work, only the long-term thermal evolution (or cooling) will be addressed. Though for convenience sometimes we shall simply call it thermal evolution (or cooling).

Under the hypothesis that no source of internal heating is present, the thermal history depends upon neutrino emission mechanisms, transport properties, and superfluidity gaps of various particle degrees of freedom (Yakovlev et al. 2001; Yakovlev & Pethick 2004; Page et al. 2006; Page 2009; Tsuruta 2009), all of which being strongly dependent on the EOS (Raduta et al. 2018; Fortin et al. 2021). Stars accommodating at least one direct Urca (DURCA) process cool much faster than stars whose thermal evolution is driven

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹National Institute for Physics and Nuclear Engineering (IFIN-HH), RO-077125 Bucharest, Romania; mikhail.beznogov@nipne.ro, araduta@nipne.ro</footnote>

<footnote n="2">²Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de México, CDMX 04510, Mexico; page@astro.unam.mx</footnote>

<footnote n="3">³Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris-Cité, F-92190 Meudon, France; jerome.novak@obspm.fr</footnote>
by the series of the so-called slow (modified Urca and bremsstrahlung), and in the case of paired species, intermediate (Cooper pair breaking and formation) cooling processes. Whether DURCA processes occur or not depends upon the density dependence of the nuclear symmetry energy; nucleation of heavy baryons and quarks; condensations of pions and kaons. Whether one or more species are paired depends upon the magnitude and density dependence of the pairing gaps of fermions, both controlled by the baryon interactions and the numerous particle correlations accounted for. Uncertainties in nuclear physics beyond the saturation density of isospin symmetric matter render the appearance of fast and intermediate cooling processes, and consequently, the thermal evolution, highly model-dependent.

Thermal evolution equations suitable for studying axisymmetric NSs with arbitrary rotation frequencies were derived from first principles by Miralles et al. (1993), who also investigated the cooling of rotating isolated NSs. In that work, the mechanical structure of the star was computed within the slow rotation approximation (Hartle 1967; Hartle & Thorne 1968), while the isothermal interior approximation was used for the thermal evolution. The latter assumption means that the redshifted temperature is uniform throughout the volume, which makes it not directly affected by rotation. In such a situation the angular dependence of surface temperature exists only due to the dependence of heat-blanketing envelope solution on surface gravitational acceleration (Gudmundsson et al. 1982). It was also shown, albeit on a schematic model, that the cold wave from the core reaches the equator later than the poles (the core is usually cooling down faster than the crust due to more intensive neutrino emission; hence, the cold wave).

More sophisticated investigations were done by Schaab & Weigel (1998), who improved upon both mechanical structure and thermal evolution calculations. Their results confirmed the cooling wave propagation conclusions of Miralles et al. (1993) and also demonstrated that fast rotation delays the crust relaxation phase by a few hundred years.

Cooling of NSs rotating at a constant frequency was later on considered by Negreiros et al. (2012). When the constant frequency assumption is relaxed, e.g., by taking into account the pulsar spindown, modifications of microscopic properties affect the cooling history (Negreiros et al. 2013). Considering that the DURCA process is only allowed above a critical density, $\rho_{\text{DU}}$, whose value is EOS dependent, a fast-rotating NS whose central density $\rho_c$ was initially below $\rho_{\text{DU}}$ may see, due to the decrease of the centrifugal force during the spindown, $\rho_c$ increase sufficiently to eventually surpass $\rho_{\text{DU}}$ and slowly turn on the DURCA process in a time-growing inner core. Simulations accounting for fully time-dependent microphysics, including superfluid gaps and mechanical structure were performed within the *thermorotational* model by Negreiros et al. (2017).

A common feature of these works is that they focus on the surface temperature distribution, which is strongly affected by the value of local surface gravitational acceleration. The figures in Negreiros et al. (2012, 2017) nevertheless suggest that internal temperature distributions are worth more attention.

Thermal evolution of NSs in 2D was also investigated in relation to magnetic fields and/or extra heat sources in the crust. Magnetothermal evolution works usually assuming a spherically symmetric (1D) mechanical structure and an axisymmetric (2D) thermal and magnetic structure (e.g.,

Vigano et al. 2021). Note that mixed 1D and 2D calculations, but this time for the core and the crust, respectively, have been employed for magnetized NSs (Geppert et al. 2004) as well as for extra heat sources in the crust (Kaminker et al. 2014).

Our project aims to generalize the studies of long-term thermal evolution of NSs beyond spherical symmetry. This will allow us to more realistically investigate rotating isolated NSs; accreting NSs that are heated anisotropically; NSs deformed under strong magnetic fields. In order to account as much as possible for the rich microphysics relevant to these topics and take benefit from the existing expertise in the field we made the choice of upgrading NSCool (Page 2016), a code in 1D, extensively used for thermal evolution simulations of a large variety of NSs among which are long-term cooling of isolated stars (Page et al. 2004, 2009), cooling of neo-neutron stars (Beznogov et al. 2020), and heating and cooling of accreting stars in low mass X-ray binaries (Ootes et al. 2016; Degenaar et al. 2021; Page et al. 2022). In this paper, we shall only consider nucleonic rigidly rotating NSs under the standard long-term cooling paradigm.

This paper is organized as follows. In Section 2, we briefly present the spacetime metric appropriate for describing rigidly rotating NSs as well as the equations of their thermal evolution. Models of EOSs and the mechanical structure of NSs rotating up to the Kepler limit are considered in Section 3. In Section 4, we present the results of our cooling simulations. Our conclusions are presented in Section 5. Throughout this paper, by matter density ($\rho$) we mean energy density including rest mass divided by the speed of light squared.

## 2. Equations for Thermal Evolution of Rotating Stars

### 2.1. Spacetime Metric

Calculation of the spacetime metric and structure of axisymmetric rigidly rotating NS in full general relativity (GR) is a complicated issue. Relatively simple procedures exist only in the limit of slow rotation, where rotation-induced modifications can be treated as small perturbations of an (already known) static configuration (Hartle 1967; Hartle & Thorne 1968). For rapidly rotating NSs more sophisticated techniques are required for solving the coupled system of mechanical structure and Einstein field equations. One of the most efficient and accurate methods was proposed by Bonazzola et al. (1993); the numerical integration is based on a spectral method. In this work, we employ the LORENE library<fnref n="4" /> (Gourgoulhon et al. 2016), suitable to deal with rotation frequencies from zero up to the mass shedding limit.

Herein, we briefly review the working hypothesis, choice of coordinates, metric potentials, and NS mechanical structure; for more information, see Bonazzola et al. (1993).

The motion of matter is assumed to be (i) axisymmetric (with respect to the rotational axis), (ii) stationary, (iii) purely azimuthal (i.e., no convection, no meridional currents; we also neglect all mass flows induced by heat transport). And (iv) the EOS of the matter is barotropic. For possible relaxation of (iii) see, e.g., Birkl et al. (2011) and of (iv) see, e.g., Camelio et al. (2019).

Following Bonazzola et al. (1993), we use the *quasi-isotropic* gauge, relying on spherical coordinates with the polar angle $\theta$ measured from the rotation axis. Note that the quasi-isotropic coordinates used in LORENE are different from the

<footnote n="4">⁴ https://lorene.obspm.fr</footnote>
conventional Schwarzschild coordinates used in the spherically symmetric Tolman–Oppenheimer–Volkoff equations.<fnref n="5" />

To write down the components of the metric tensor we need to additionally specify the choice of the basis tetrad (basis vectors). Choosing natural tetrad $e_\alpha = \partial_\alpha = \partial/\partial x^\alpha$ with $x^\alpha = (ct, r, \theta, \phi)^T$, we obtain the following expression in the notation of Gourgoulhon et al. (1999)<fnref n="6" />:

$$ g_{\alpha\beta} = \begin{pmatrix} -N^2 + B^2 \left(\frac{N^\varphi}{c}\right)^2 r^2 \sin^2\theta & 0 & 0 & -B^2 \left(\frac{N^\varphi}{c}\right) r^2 \sin^2\theta \\ 0 & A^2 & 0 & 0 \\ 0 & 0 & A^2 r^2 & 0 \\ -B^2 \left(\frac{N^\varphi}{c}\right) r^2 \sin^2\theta & 0 & 0 & B^2 r^2 \sin^2\theta \end{pmatrix}, \quad (1) $$

where $N, A, B$, and $N^\varphi$ are four metric functions that depend on $r$ and $\theta$ but, because of stationarity and axisymmetry, do not depend on $t$ and $\varphi$; $c$ is the speed of light; $\alpha, \beta = 0, 1, 2$, and 3. Note that $N, A$, and $B$ are dimensionless while $N^\varphi$ has the dimension of frequency. The line element reads as

$$ ds^2 = g_{\alpha\beta} dx^\alpha dx^\beta = -N^2 c^2 dt^2 + A^2 (dr^2 + r^2 d\theta^2) + B^2 r^2 \sin^2\theta (d\varphi - N^\varphi dt)^2. \quad (2) $$

## 2.2. Heat Transport Equations

Energy balance and energy transport equations for rotating axisymmetric systems in GR were derived by Miralles et al. (1993) by means of the relativistic Boltzmann equation. Here, we write them down following Schaab & Weigel (1998), with the notation adjusted to match that in Gourgoulhon et al. (1999) (see the previous subsection; see also Negreiros et al. 2012).

The energy balance (or conservation) equation is<fnref n="7" />:

$$ \partial_r \tilde{h}_{(r)} + \frac{1}{r} \partial_\theta \tilde{h}_{(\theta)} = -r^2 \sin \theta A^2 B \Gamma \times \left[ -\left(\frac{N}{\Gamma}\right)^2 Q + C_V \partial_t \tilde{T} \right], \quad (3) $$

where $C_V$ is the specific heat, $Q$ stands for the heat sources/ sinks, and

$$ U = Br \sin \theta (\Omega' - N^\varphi)/N, \quad (4) $$

$$ \Gamma = \frac{1}{\sqrt{1 - \left(\frac{U}{c}\right)^2}}, \quad (5) $$

$$ \tilde{h}_{(i)} = \frac{ABN^2}{\Gamma} r^2 \sin \theta h^i, \quad i = r, \theta, \quad (6) $$

$$ \tilde{T} = \frac{N}{\Gamma} T, \quad (7) $$

with $T$ being the local temperature and $\Omega'$ being the rotation frequency (in radians per second; $\Omega' = 2\pi\Omega$). Redshifted

quantities are denoted with a tilde. In Equation (6) the brackets in the components of the redshifted flux, $\tilde{h}_{(i)}$, emphasize that these do not form a vector. In contrast, the components of the local flux $h^i$ in the same equation do form a vector and this equation defines contravariant components, as indicated by the superscript. Also note that the local flux is defined in the comoving reference frame.

The energy flux is calculated in the diffusion approximation:

$$ \tilde{h}_{(r)} = -r^2 \sin \theta \text{KNB} \partial_r \tilde{T}, \quad (8) $$

$$ \tilde{h}_{(\theta)} = -r \sin \theta \text{KNB} \partial_{\theta} \tilde{T}, \quad (9) $$

where $K$ is the thermal conductivity; in this work we assumed that $K$ is isotropic (scalar).

Equations (3), (8), and (9) can be combined into a single second-order partial differential equation for the evolution of the redshifted temperature $\tilde{T}$:

$$ \partial_t \tilde{T} = \frac{1}{r^2 \sin \theta A^2 B \Gamma C_V} \left[ \partial_r (r^2 \sin \theta \text{KNB} \partial_r \tilde{T}) + \partial_\theta (\sin \theta \text{KNB} \partial_\theta \tilde{T}) \right] + \left(\frac{N}{\Gamma}\right)^2 \frac{Q}{C_V}. \quad (10) $$

Equation (10) can be further recast as

$$ \partial_t \tilde{T} = G_{(rr)} + G_{(\theta\theta)} + \left(\frac{N}{\Gamma}\right)^2 \frac{Q}{C_V}, \quad (11) $$

$$ G_{(rr)} = \frac{1}{r^2 A^2 B \Gamma C_V} \partial_r (r^2 \text{KNB} \partial_r \tilde{T}), \quad (12) $$

$$ G_{(\theta\theta)} = \frac{1}{r^2 \sin \theta A^2 B \Gamma C_V} \partial_\theta (\sin \theta \text{KNB} \partial_\theta \tilde{T}). \quad (13) $$

Equation (11) is our main equation. It has the form of a classical heat diffusion equation. It is parabolic and allows for instantaneous heat propagation. One may question the applicability of such an equation in the framework of GR. The issue has been previously addressed by Dommes et al. (2020). According to these authors, instantaneous propagation of heat and GR are not incompatible as long as the hydrodynamic description is applicable. The NSs we are interested in, i.e., with thermal evolution timescales much longer than hydrodynamic timescales, fall into this category, which makes the use of Equation (11) safe.

The solution of Equation (11) requires initial and boundary conditions, which we comment on in the following.

## 2.3. Initial and Boundary Conditions

The initial temperature profile at the initial time is chosen, for all NS, as follows. Matter in the center of the star and at the bottom of the envelope is given redshifted temperatures of $1.5 \times 10^{10}$ K and $10^{10}$ K, respectively. In the intermediate shells, $\tilde{T}$ is assumed to scale with the logarithm of density. Such a setup corresponds roughly to constant internal flux and improves numerical convergence of the first time step as it makes the initial and boundary conditions more compatible; starting with constant redshifted temperature would have meant zero initial internal flux, which is incoherent with nonzero initial surface flux.<fnref n="8" /> Note that the information on the initial temperature distribution is rapidly lost, which means that the

<footnote n="5">⁵ In spherical symmetry, the quasi-isotropic gauge turns into the isotropic gauge, for which there exists a coordinate transform to the Schwarzschild one. Note, however, that outside spherical symmetry there is no gauge generalizing the Schwarzschild coordinates.</footnote>

<footnote n="6">⁶ In (Bonazzola et al. 1993) some components of the metric tensor are defined differently, while the notation of Gourgoulhon et al. (1999) corresponds exactly to the one used in LORENE.</footnote>

<footnote n="7">⁷ For neutrino losses $Q_\nu < 0$, the sign before the source term is chosen accordingly.</footnote>

<footnote n="8">⁸ For neo-NS modeling a dedicated matching procedure is required to have fully compatible initial and boundary conditions, see (Beznogov et al. 2020).</footnote>
<table><caption>**Table 1**

Nuclear and Astrophysical Properties of Our Three EOS Models</caption><thead><tr><th>Model</th><th><i>E</i><sub>b</sub><br>(MeV)</th><th><i>n</i><sub>sat</sub><br>(fm<sup>-3</sup>)</th><th><i>K</i><sub>sat</sub><br>(MeV)</th><th><i>E</i><sub>sym</sub><br>(MeV)</th><th><i>L</i><sub>sym</sub><br>(MeV)</th><th><i>K</i><sub>sym</sub><br>(MeV)</th><th><i>M</i><sub>G,max</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>R</i><sub>14</sub><br>(km)</th><th>Λ̃</th><th><i>n</i><sub>DU</sub><br>(fm<sup>-3</sup>)</th><th><i>M</i><sub>DU</sub><br>(<i>M</i><sub>☉</sub>)</th></tr></thead><tbody><tr><td>APR</td><td>-16.00</td><td>0.160</td><td>266.0</td><td>32.6</td><td>57.6</td><td>-102.6</td><td>2.17</td><td>11.33</td><td>272-299</td><td>0.78</td><td>2.00</td></tr><tr><td>IUF</td><td>-16.40</td><td>0.155</td><td>231.3</td><td>31.3</td><td>47.2</td><td>28.5</td><td>1.95</td><td>12.64</td><td>499-530</td><td>0.61</td><td>1.77</td></tr><tr><td>NL3-ωρ</td><td>-16.24</td><td>0.148</td><td>271.6</td><td>31.7</td><td>55.5</td><td>-8.1</td><td>2.75</td><td>13.82</td><td>1042-1051</td><td>0.50</td><td>2.55</td></tr></tbody></table>

**Note.** Listed are binding energy per nucleon ($E_b$); compression modulus ($K_{sat}$); symmetry energy ($E_{sym}$), its slope ($L_{sym}$), and curvature ($K_{sym}$) at the saturation point of symmetric matter ($n_{sat}$). The remaining columns correspond to the maximum NS gravitational mass ($M_{G,max}$); radius of $M_G = 1.4 M_☉$ NS in Schwarzschild gauge; range of combined tidal deformability for the GW170817 event; particle number density for the DURCA process threshold ($n_{DU}$) and the corresponding NS gravitational mass ($M_{DU}$).

initial thermal state is irrelevant for long-term evolution studies (Beznogov et al. 2020).

Boundary conditions at the origin and on the symmetry axis have mathematical motivation. They consist in requiring that $\tilde{T}$ is regular and axisymmetric. Calculation of $\tilde{T}$ in the limit $r \to 0$ requires that one integrate Equation (11) over a small spherical volume centered on $r=0$.

At variance with this, surface boundary conditions have physical motivation. Following the standard approach in studies of NS cooling, we have chosen the boundary density $\rho_b$ of the heat-blanketing envelope (i.e., its bottom) to be $10^{10} \text{ g cm}^{-3}$ (Potekhin et al. 2015). This is the outer boundary where we impose the boundary condition (see Appendix A for more details).

The surface emission is determined by projecting the heat flux vector $h^i$, taken at the outer boundary, onto the normal to this boundary.

From Equations (6), (8), and (9) we obtain

$$ \frac{\partial \tilde{T}}{\partial r} = - \frac{A N}{K \Gamma} h^r, \qquad (14) $$

$$ \frac{\partial \tilde{T}}{\partial \theta} = - \frac{A N}{K \Gamma} r h^{\theta}. \qquad (15) $$

Now we need to define the unit normal vector (we are only interested in $r$ and $\theta$ covariant components):

$$ n_r = A \left[ 1 + \frac{1}{R^2} \left( \frac{dR}{d\theta} \right)^2 \right]^{-\frac{1}{2}}, \qquad (16) $$

$$ n_{\theta} = -A \left[ 1 + \frac{1}{R^2} \left( \frac{dR}{d\theta} \right)^2 \right]^{-\frac{1}{2}} \frac{1}{R} \frac{dR}{d\theta}, \qquad (17) $$

where $R=R(\theta)$ defines the radial coordinate of the outer boundary of the star, which by definition coincides with the bottom of the heat-blanketing envelope (see Appendix A for details). Since the heat flux is a contravariant vector, we can directly compute the scalar product:

$$ h_s = \vec{n} \cdot \vec{h} = K \frac{\Gamma}{N} \left[ 1 + \frac{1}{R^2} \left( \frac{dR}{d\theta} \right)^2 \right]^{-\frac{1}{2}} \times \left[ -\left. \frac{\partial \tilde{T}}{\partial r} \right|_{r=R} + \left. \frac{\partial \tilde{T}}{\partial \theta} \right|_{r=R} \frac{1}{R^2} \frac{dR}{d\theta} \right]. \qquad (18) $$

On the other hand, the surface flux $h_s$ can be computed via the heat-blanketing solution. In the plane-parallel approximation at

any angle $\theta^9$,

$$ h_s(\theta) = \sigma_{SB} T_s^4(\theta), \quad T_s(\theta) = T_s(T_b(\theta)), \qquad (19) $$

where $T_s(T_b)$ is the $T_s-T_b$ relation. For the latter, we employ the envelope model of Potekhin et al. (1997) with the parameter $\eta$ set to $10^{-18}$. This is equivalent to an iron (non-accreted) heat-blanketing envelope. Combining Equations (18) and (19) one obtains the surface boundary condition.

We have developed a code based on this formalism and called it NSCool 2D Rot to emphasize inheritance from NSCool, a publicly available NSs thermal evolution code in 1D (Page 2016). For more details, see Appendix B.

### 3. Models of Rotating Stars

#### 3.1. EOS

In order to identify the dependence on EOS, simulations are conducted for NSs built upon different models for the core EOS. For the sake of completeness, results corresponding to different models for the crust EOS will be commented as well.

For the core, we consider the following models: APR (Akmal et al. 1998), IUF (Fatoyev et al. 2010), and NL3-ωρ (Horowitz & Piekarewicz 2001; Pais & Providência 2016). Some of their characteristics, translated into values of symmetric saturated nuclear matter parameters and properties of cold-catalyzed NSs, are presented in Table 1. Present uncertainties in the isoscalar and isovector channels, due to limited information from nuclear physics experiments, are accounted for by the ranges spanned by the parameters $K_{sat}$ and $L_{sym}$ and $K_{sym}$, respectively. The models with the largest values of $K_{sat}$ predict maximum NS masses much exceeding the lower limit on NS maximum masses, $\approx 2 M_☉$, while the model with the lowest value of $K_{sat}$ is marginally consistent with this constraint. Radii of canonical mass NS, largely governed by the symmetry energy over $2n_{sat} \lesssim n_B \lesssim 3n_{sat}$, span a 2.5 km interval. The density dependence of the symmetry energy in the supra-saturation regime explains also the largely different density thresholds for the opening of DURCA processes. We nevertheless note that, due to the limited values of $L_{sym}$, none of these EOSs allows for DURCA processes to operate in NSs with canonical masses (Fortin et al. 2016). The less massive NS that allows for nucleonic DURCA process corresponds to the IUF EOS, i.e., the model with the largest value of $K_{sym}$.

Two EOSs are alternatively used for the crust. The first one corresponds to HZD (Haensel et al. 1989) supplemented by NV (Negele & Vautherin 1973) models, covering the density

<footnote n="9">⁹ One should differentiate between local and redshifted quantities (the latter are denoted by a tilde).</footnote>
<table><caption>**Table 2**
Properties of Static and Rigidly Rotating NSs with $M_B = 1.6 M_\odot$ Built upon Different EOS Models</caption><thead><tr><th>Model</th><th>Ω</th><th>ρ<sub>c</sub></th><th>M<sub>G</sub></th><th>R<sub>p</sub><sup>Co/Ph</sup></th><th>R<sub>e</sub><sup>Co/Ph</sup></th><th>ΔR<sub>p</sub><sup>Co/Ph</sup></th><th>ΔR<sub>e</sub><sup>Co/Ph</sup></th><th>ζ<sup>Co</sup></th><th>ζ<sub>cr</sub><sup>Co</sup></th><th>g<sub>s14,p</sub></th><th>g<sub>s14,e</sub></th></tr></thead><tbody><tr><td>APR</td><td>0</td><td>1.034</td><td>1.437</td><td>8.96/12.62</td><td>8.96/12.62</td><td>0.82/1.04</td><td>0.82/1.04</td><td>1.000</td><td>1.000</td><td>1.88</td><td>1.88</td></tr><tr><td>APR</td><td>716</td><td>0.989</td><td>1.446</td><td>8.46/11.93</td><td>9.57/13.39</td><td>0.80/1.02</td><td>1.02/1.28</td><td>0.884</td><td>0.785</td><td>1.92</td><td>1.44</td></tr><tr><td>APR</td><td>1096</td><td>0.890</td><td>1.465</td><td>7.61/10.69</td><td>11.85/16.11</td><td>0.80/1.03</td><td>2.21/2.69</td><td>0.642</td><td>0.363</td><td>1.92</td><td>0.26</td></tr><tr><td>IUF</td><td>0</td><td>0.817</td><td>1.452</td><td>10.15/13.80</td><td>10.15/13.80</td><td>0.94/1.16</td><td>0.94/1.16</td><td>1.000</td><td>1.000</td><td>1.50</td><td>1.50</td></tr><tr><td>IUF</td><td>716</td><td>0.740</td><td>1.463</td><td>9.37/12.72</td><td>11.27/15.13</td><td>0.92/1.14</td><td>1.31/1.60</td><td>0.831</td><td>0.700</td><td>1.53</td><td>0.98</td></tr><tr><td>IUF</td><td>932</td><td>0.657</td><td>1.478</td><td>8.61/11.62</td><td>13.48/17.69</td><td>0.93/1.17</td><td>2.49/2.97</td><td>0.639</td><td>0.375</td><td>1.51</td><td>0.20</td></tr><tr><td>NL3-ωρ</td><td>0</td><td>0.525</td><td>1.463</td><td>11.40/14.98</td><td>11.40/14.98</td><td>1.39/1.69</td><td>1.39/1.69</td><td>1.000</td><td>1.000</td><td>1.22</td><td>1.22</td></tr><tr><td>NL3-ωρ</td><td>716</td><td>0.489</td><td>1.478</td><td>10.05/13.23</td><td>13.14/17.06</td><td>1.34/1.64</td><td>2.19/2.61</td><td>0.765</td><td>0.612</td><td>1.27</td><td>0.62</td></tr><tr><td>NL3-ωρ</td><td>826</td><td>0.470</td><td>1.486</td><td>9.47/12.44</td><td>14.80/18.99</td><td>1.34/1.65</td><td>3.28/3.87</td><td>0.640</td><td>0.409</td><td>1.27</td><td>0.16</td></tr></tbody></table>

**Note.** Different rotation frequencies (Ω) in hertz, specified in the second column, are assumed. The largest value corresponds, in each case, to 99% of the Kepler frequency. Listed are central density (ρ<sub>c</sub>) in units of 10<sup>15</sup> g cm<sup>-3</sup>; gravitational mass (M<sub>G</sub>) in units of (M<sub>c</sub>); radius (R) and crust thicknesses (ΔR) in kilometers; surface gravitational acceleration (g<sub>s14</sub>) in units of (10<sup>14</sup> cm c<sup>-2</sup>). Subscripts “p” and “e” stand for “polar” and “equatorial,” respectively. Superscripts “Co” and “Ph” stand for “coordinate” and “physical” (radii and thicknesses). The polar to equatorial axis ratio is denoted by ζ while ζ<sub>cr</sub> stands for polar to equatorial crust thickness ratio. As a reminder, quasi-isotropic coordinates differ from Schwarzschild ones. See the text for details.

domains below and above the neutron drip density, respectively. The second one corresponds to recent calculations by Mondal et al. (2020) and Viñas et al. (2021); it is based on the finite range Gogny interaction D1M and accounts for shell and pairing corrections. These two EOS models are almost identical in regard to the P(ρ) dependence but provide significantly different chemical compositions in the inner crust. More precisely, nuclei are more neutron-rich in D1M than in HZD +NV, which entails a correspondingly lower density of unbound neutrons in the first model with respect to the second. Differences in the treatment of shell and pairing effects in addition to in-medium modification of cluster energy functionals are reflected also in different density dependencies of the atomic numbers though to a lesser extent. Core and crust EOSs are smoothly matched at n<sub>B</sub> ≈ n<sub>sat</sub>/2. At this point, we note that the calculation of fast-rotating stars by LORENE requires the EOS to be sufficiently smooth.

## 3.2. Mechanical Structure

Before studying thermal evolution let us consider equilibrium configurations of rigidly rotating NSs. We choose to compare NSs with equal values of baryonic mass, a quantity that is conserved if accretion and mass loss are disregarded. Note that other choices, e.g., equal values of gravitational mass or central mass density, are also possible and have been used in the literature.

The value of choice for the baryonic mass in this prospective study is $M_B = 1.6 M_\odot$, which results in gravitational masses close to the canonical value of $1.4 M_\odot$. Two values of the rotation frequency are considered. The first one, 716 Hz, corresponds to the fastest-spinning millisecond pulsar known to date (PSR J1748–2446ad; see Hessels et al. 2006). The second one corresponds to 99% of the Kepler frequency and is EOS dependent.

Table 2 lists some of the properties of NSs. For completeness characteristics of static stars are provided as well. NS radii and crust thicknesses are reported along the equatorial (“e”) and polar (“p”) directions. Two values are provided. The first corresponds to the quasi-isotropic spherical coordinates used by LORENE (superscript “Co”); the second corresponds to what a local observer would measure (superscript “Ph”), and thus, refers to proper (physical) lengths. Note that (i) coordinate lengths measured in quasi-isotropic spherical coordinates

cannot be compared to those derived in Schwarzschild coordinates, and (ii) proper (physical) lengths are the same regardless of the choice of coordinates. Proper radial distances (i.e., with t, θ, ϕ = const) can be computed as

$$l_r = \int_{r_{\min}}^{r_{\max}} A(r', \theta) dr'. \qquad (20)$$

A word of caution is needed with regard to the radii and crust thicknesses mentioned in Table 2. Both are measured up until the bottom of the envelope, i.e., up until ρ<sub>b</sub> = 10<sup>10</sup> g cm<sup>-3</sup>. See Appendix A for more details.

Let us first focus on static (i.e., nonrotating spherically symmetric) configurations. The following remarks are in order: (i) the highest (lowest) central density value is obtained for the most (least) compact configuration, which in turn corresponds to the softest (stiffest) EOS model; note that characterization of the EOS as soft or stiff in relation to $M_B = 1.6 M_\odot$ NS refers to the slope of the P(ρ) curve over the density domain span in these stars, and in general, differs from what is obtained when much wider domains are considered, as it is the case of the maximum mass configuration (see Table 1); (ii) as already discussed in relation to Table 1 uncertainties in the density dependence of the symmetry energy in the supra-saturation domain lead to substantial variations in the ranges of predicted radii; (iii) the scatter between the values of gravitational mass is very limited, less than 4%; (iv) even if the same EOS model is used for the crust, its thickness changes from one model to another; (v) the EOS model that provides the largest (smallest) value for the radius produces the thickest (thinnest) crust.

Let us now consider the case of the mass shedding limit. The largest (smallest) frequency (Ω<sub>K</sub>) corresponds to the APR (NL3-ωρ) EOS. These are the EOSs with the most (least) compact static configuration. The dispersion of Ω<sub>K</sub> is significant while that of the polar to equatorial axis ratio (ζ) is insignificant. Indeed, the Kepler frequency of the APR EOS model exceeds that of the NL3-ωρ EOS model by 25%, while the difference between the axis ratios is just 0.3%. One notes that the APR and NL3-ωρ EOSs provide for both static and maximally rotating stars the smallest and largest values of radii, respectively. Central densities of maximally rotating configurations are lower by 10% (NL3-ωρ EOS) to 20% (IUF EOS) than those of the static configurations and no clear correlation between Ω<sub>K</sub> and ρ<sub>c</sub> can be established. Being dilute the crust is
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Mechanical structure of rotating NSs built upon the APR EOS. Left (right) panel: $\Omega = 716$ Hz (1096 Hz) model. The contour lines show the density distribution as a function of coordinate radius $r$ (in kilometers) and polar angle $\theta$ (in radians). Values on the contours indicate the logarithm of density. The light yellow-shaded area indicates the crust. The crust-core boundary is depicted by a magenta line "14.23" and blue contour "10" depicts the bottom of the heat-blanketing envelope, which is the outer boundary of the star.</figcaption></figure>

much more affected by rotation than the core: at the equator, it increases by 130% (NL3-$\omega\rho$ EOS) to 160% (APR EOS); the crust at the pole is much thinner than the crust at the equator.

NSs rotating with 716 Hz feature properties intermediate between those of static and maximally rotating stars. As in nonrotating and maximally rotating NSs, the largest (smallest) value of central density corresponds to the APR (NL3-$\omega\rho$) EOS and the most (least) deformed configuration is obtained for the NL3-$\omega\rho$ (APR) EOS. The dispersion among the values of $\zeta$ in NSs built upon different EOSs is rather important at 716 Hz, while practically no dispersion is obtained when approaching Kepler frequency. Quite interestingly, the ratios between crust thickness at the pole and at the equator, $\zeta_{cr}$, produced by various EOSs are in a different relation than their counterparts at the Kepler limit.

Further insight into the geometrical properties of rotating NSs and matter distribution is shown in Figures 1 and 2. The considered EOS model is APR; the rotation frequencies are mentioned in the key legend and figure caption. Contour lines in Figure 1 show the spatial distribution of matter density; values on the plot correspond to the logarithm of density. The yellow-shaded area indicates the crust; it extends from the crust-core boundary (magenta line marked "14.23") to the bottom of the heat-blanketing envelope, which coincides with the outer boundary of the star (blue contour "10"). The star rotating at 716 Hz (1096 Hz) deviates little (significantly) from spherical symmetry. Consequently, crust thickness is marginally (sizably) larger at the equator than at the pole. All these features are in agreement with the data in Table 2. We anticipate that these differences in the crust thickness will translate into different heat diffusion timescales that will modify the cooling pattern in the equatorial direction compared to the polar direction.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Radial profiles of density along the polar (solid lines) and equatorial (dashed lines) directions for a static NS (blue line) as well as NSs rotating at 716 Hz (red lines) and 1096 Hz (green lines). The considered EOS model is APR. The horizontal dotted line marks the crust-core transition density.</figcaption></figure>

Figure 2 first confirms what we have already seen in Table 2. Namely, that rapid rotation makes NSs less compact. It also shows that density gradients along the equatorial direction are smaller than those along the polar direction.

Before concluding this section some comments are in order. First, we assume that $\Omega = \text{const}$ throughout the cooling
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Redshifted temperature at the bottom of the envelope (left panel) and on the surface (right panel) as a function of time (cooling curves). Results correspond to static NSs with $M_B = 1.6 M_\odot$ built upon the EOS models indicated in the legend.</figcaption></figure>

process, which is equivalent to discarding any spindown, e.g., due to magnetic braking (Negreiros et al. 2017). It is a reasonable simplification as magnetic fields in isolated NSs are not necessarily strong. The second comment deals with the highest rotation frequency we have considered.

For an NS rotating at near-Kepler frequency any small deviation from axisymmetric configuration immediately leads to rapid loss of angular momentum due to GW emission. These deviations can be caused by various reasons: r-mode instability, f-mode-instability, thermal, or magnetic mountains, etc. Both theoretical computations and the fact that there are currently no known pulsars rotating faster than 716 Hz, suggest that the spindown due to GW emission is extremely efficient and that NSs rotating at the near-Kepler limit cannot maintain their rotation frequency long term. Thus, studying long-term cooling of such stars is academic. We do not expect to find any. Still, we included them in our analysis as a limiting case opposite to the limiting case of nonrotating spherically symmetric stars. We want to check to the maximum possible extent the effects caused by rotation.

Mechanical properties of NSs built upon the D1M EOS of the crust are identical to those obtained when the HDZ+NV EOS is used. The reason is that the $P(\rho)$ dependence of both EOSs is virtually the same. Note nevertheless that differences in chemical composition lead to slight differences in the thermal evolution, see Section 4.

## 4. Cooling of Rotating Isolated NSs

We shall consider here the standard long-term cooling of purely nucleonic isolated NSs with $M_B = 1.6 M_\odot$ built upon the EOS models described in Section 3.1. The results presented in this section are obtained using the HZD+NV EOS for the crust. We have also performed cooling simulations with the D1M crust EOS and we shall briefly comment on them at the end of this section. As mentioned in Section 1, we solve mechanical structure equations only once for each model (using LORENE, see Section 2.1) and then solve only thermal evolution equations at each time step (using NSCool 2D Rot, see Appendix B) (keeping the structure fixed throughout the evolution).

Moderate values of $L_{\text{sym}}$ and NS mass ensure that none of the considered NS models allows for the DURCA process. Nucleonic pairing in the crust and core is disregarded for simplicity. This means that neutrino cooling by the formation and breaking of Cooper pairs does not occur either. So, the main driving force behind the thermal evolution is the modified Urca process, and at later stages, surface photon emission. The neutrino emission rates, thermal conductivities, and specific heat capacity are standard and described by, e.g., Yakovlev & Urpin 1980; Yakovlev et al. 2001; Haensel et al. 2007; Shternin & Yakovlev 2007; for implementation see (Page et al. 2004, 2011). Magnetic fields and heat sources are disregarded.

To start with let us consider the limiting case of static NSs. Redshifted temperatures at the bottom of the envelope, $\tilde{T}_b$, and on the surface, $\tilde{T}_s$, are plotted as a function of time in Figure 3, left and right panels, for the different EOS models. The left panel shows that during the phases of relaxation from initial conditions,$^{10}$ (stage I with $\log t[\text{yr}] \lesssim -0.5$), and early plateau (stage II with $-0.5 \lesssim \log t[\text{yr}] \lesssim 1.5$), the various stars have almost identical values of the bottom temperature. The crust relaxation phase lasts for the longest (shortest) period for the NS built upon the NL3-$\omega\rho$ (APR) EOS. The similarity between $\tilde{T}_b$ can be explained by the identical composition of the crusts, while the duration of the crust relaxation phase reflects the crust thickness. Indeed, as shown in Table 2, static stars built upon the NL3-$\omega\rho$ and APR EOSs have the thickest and thinnest crusts, respectively, which means that the cold wave from the core needs less time to pass through the crust in the case of the NS with the APR EOS than in the one of the NS with the NL3-$\omega\rho$ EOS. The model that is the first (last) one to leave the early plateau phase is characterized by the lowest (highest) $\tilde{T}_b$ during the subsequent crust-core thermalization (stage III with $1.5 \lesssim \log t[\text{yr}] \lesssim 2.5$) and neutrino cooling (stage IV for $t[\text{yr}] \gtrsim 2.5$) stages.

Comparing the left and right panels, one can judge the impact of the heat-blanketing envelope solution that is used as

<footnote>$^{10}$ We will refer to everything happening before the early plateau phase as relaxation from initial conditions. Note, however, that the actual relaxation happens much faster and is followed by some other (sub-)phases before the early plateau phase begins. See (Beznogov et al. 2020) for details.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Redshifted temperature at the bottom of the envelope (left panel) and on the surface (right panel) as a function of time (cooling curves). Results correspond to the APR EOS and different rotation frequencies specified in the legend. Considered are values at the pole ($\theta = 0$) and at the equator ($\theta = \pi/2$).</figcaption></figure>

an outer boundary condition in the cooling calculations. $\tilde{T}_s(t)$ produced by our NS models are identical in stages I and IV; in stage II the highest (lowest) values of $\tilde{T}_s$ correspond to the APR (NL3-$\omega\rho$) EOS; in stage III the ranking is inverted. The rearrangement of cooling curves is due to the dependence of the envelope solution on surface gravitational acceleration, see Table 2.

The modification of $\tilde{T}_{b/s}(t)$ under rotation is investigated in Figure 4 for the particular case of the APR EOS. Redshifted values of the temperature at the bottom of the envelope (left panel) and on the surface (right panel), corresponding to $\Omega = 716$ and $1096$ Hz, are confronted at $\theta = 0$ (pole) and $\pi/2$ (equator). Also considered is the case of a static star. From the left panel, it turns out that $\tilde{T}_{b,p}(t)$ is affected by rotation only in the photon cooling stage (V) and for extreme rotation frequencies. This is not very surprising as the polar regions are less affected by rotation and heat transfer there is mostly radial (unlike equatorial regions, see below). At variance with this, a series of modifications affect $\tilde{T}_{b,e}(t)$: the redshifted value of the temperature at the bottom of the envelope is slightly enhanced in stage I; stage II is delayed and lasts longer; $\tilde{T}_{b,e}^{(II)} < \tilde{T}_{b,p}^{(II)}, \tilde{T}_{b,e}^{(III)} > \tilde{T}_{b,p}^{(III)}$, stage IV is delayed. The increase in the duration of stage II is related to the increased thickness of the crust in equatorial regions: a thicker crust can store more energy and also has a larger heat diffusion time (which goes as the square of the length scale). As a result, the cold wave from the core needs more time to punch through the crust. This explanation was previously put forward by Miralles et al. (1993) and Schaab & Weigel (1998).

The right panel of Figure 4 demonstrates that gravitational acceleration on the surface is accountable for much stronger modifications of $\tilde{T}_{s,e}$ under rotation. At variance with this, no rotation-induced effect is visible on $\tilde{T}_{s,p}$. Overall and over the whole time span redshifted temperature at the equator is lower than the one at the pole and the difference between the cooling curves increases with $\Omega$.

A word of caution is in order here. For nonrotating NSs, the redshifted surface temperature is the one that—with the exception of effects related to atmosphere structure, interstellar absorption, etc.—a distant observer would see. Some authors

assume that this is also the case for rotating NSs. This is however not obvious. Indeed, Vincent et al. (2018) demonstrated that the observed spectra depend on the line of sight and that the spectral distortion, i.e., deviation from color-corrected blackbody spectrum, is small for uniform surface temperature distribution. This obviously does not apply for a fast-rotating NS. A rigorous approach to study this question would be to perform general-relativistic ray tracing, similar to that performed in Vincent et al. (2018), but this is beyond the scope of the current study.

A detailed look at the internal redshifted temperature distribution is offered in Figure 5. The considered stars are built upon the APR EOS; they rotate at $\Omega = 716$ Hz (left column) and $\Omega = 1096$ Hz (right column). Notice that the radial coordinate is the logarithm of density, which hinders information on star deformation. Different instants in the cooling evolution are considered: $t = 10$ yr (top), $t = 150$ yr (middle), and $t = 300$ yr (bottom panels).

One can see that at the age of 10 yr (stage II) the temperature gradient is mostly radial and only slightly affected by the rotation. At the age of 150 yr (stage III) the situation is very different. In the core and innermost crust (densities around $10^{14}$ g cm$^{-3}$) the temperature gradient is still mostly radial while at lower densities ($10^{10} \lesssim \rho \lesssim 10^{12}$ g cm$^{-3}$) in the equatorial zone the temperature gradient is mostly angular. A *heat blob* appears at the equator. It acts almost as a *point source* of heat that flows in all directions from it. For the faster rotating star the effect is more pronounced and the temperature distribution more complex. At the age of 300 yr (stage III) the situation is qualitatively similar to the one at 150 yr. However, for the slow-rotating star, the temperature distribution is already almost uniform, as one would expect from a star that enters the neutrino cooling phase; for the fast-rotating star, the effects of the heat blob are still present in the equatorial regions producing both radial and angular temperature gradients.

Despite complicated temperature distributions, the main idea behind the heat blob is rather simple. During phases I–III thermal evolution of the crust is governed both by its heat conduction and neutrino emission. Rotation causes the crust
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Redshifted temperature distributions in rotating NSs built upon the APR EOS. Left (right) panels correspond to Ω = 716 Hz (1096 Hz). The radial coordinate is translated into the logarithm of density. Top, middle, and bottom panels correspond to the ages of 10, 150, and 300 yr, respectively.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Profiles of redshifted temperature in an $M_B = 1.6 M_⊙$ NS built upon the APR EOS. Results corresponding to different rotation frequencies, indicated in the legend, are confronted at the pole and at the equator at $t = 150$ yr. Representations as a function of coordinate radius (logarithm of density) are provided in the left (right) panel. The two vertical dotted lines in the right panel indicate the neutron drip density and the crust-core transition density, respectively.</figcaption></figure>

to be thicker in the equatorial regions. Thus, the total heat capacity of the crust in this region increases due to the increase in its volume. Total neutrino losses increase for the same reason. Yet the heat exchange can only increase proportionally to the increase in the area of interfaces between the equatorial region and the rest of the star. Obviously, the increase in volume dominates over that of the area. Thus, the equatorial region becomes a reservoir of heat simply because it can store more thermal energy and cannot transfer it to the rest of the star fast enough (heat conduction is limited by the area).

Figure 6 further investigates the distribution of redshifted temperature in rotating stars built upon the APR EOS at $t = 150$ yr, previously considered in the middle panels of Figure 5. Results along the polar and equatorial directions for $\Omega = 716$ and 1096 Hz are confronted with those of static stars. Representations as functions of radial coordinate (density) are shown in the left (right) panel. The core is already isothermal: there is no radial temperature gradient and the values at $\theta = 0$ and $\theta = \pi/2$ coincide. The situation at the crust is the opposite: a non-monotonic radial profile is obtained and $\tilde{T}_p < \tilde{T}_e$. The temperature along the polar direction depends only loosely on the rotation frequency while a strong rotation frequency dependence manifests along the equatorial direction. In the left panel, the displacement on the X-axis between curves corresponding to different directions and rotation frequencies is evocative of NS deformation. The fact that the temperature along the equator is higher than the one at the pole is obviously due to the heat blob. But, as we shall see later, even higher redshifted temperatures are obtained at certain instants over $0 < \theta < \pi/2$. Their values as well as those of $\theta$ obviously depend on the EOS, rotation frequency, and time. The general feature is that the dispersion in $\tilde{T}(r, \theta, t)$ increases with the rotation frequency.

The angular distribution of the redshifted temperature at the bottom of the envelope and different instants is depicted in Figure 7; the same rotating stars as in Figure 6 are considered. The instants for which the curves are plotted are different in the

left and right panels; they correspond either to the early plateau (stage II) or crust-core thermalization (stage III) stages and show the same patterns. For early and mid-stage II, the temperature at the equator is lower than at the pole. As soon as the stars enter stage III the situation becomes inverted. At intermediate instants, $\tilde{T}_b(\theta)$ shows a complex structure and the highest temperature is obtained for $0 < \theta < \pi/2$. The excursion of $\tilde{T}_b(\theta)$ augments with the rotation frequency; the time the NS needs to enter the neutrino cooling phase (stage IV) does the same. All these features agree with what we have shown in Figure 4.

The ratio between redshifted temperatures<fnref n="11" /> at the pole and at the equator is represented as a function of time in Figure 8 (see also Figure 2 of Schaab & Weigel 1998). A wide time span is considered. Left (right) panels depict values at the bottom of the envelope (on the surface). In addition to the NS built upon the APR EOS, we have also plotted the curves for other EOSs here. Similarly to the previous figures results at 716 Hz and 99% of the Kepler frequency are depicted. Both panels confirm the complex behavior of redshifted temperature as a function of time, rotation frequency, and polar angle as previously shown in Figure 7, which has considered a much narrower time span. They also show a significant EOS dependence. Again, in agreement with Figure 7 the amplitude of the variation of $\tilde{T}_{b/s,p}/\tilde{T}_{b/s,e}$ increases with the rotation frequency. All curves manifest two minima and one maximum; the instants when these are obtained are the same for the quantities represented in the left and right panels. The explanation relies on the fact that the $T_s-T_b$ relation is instantaneous. The two minima are obtained for $\log t[\text{yr}] \approx -1.5$ and $2 \lesssim \log t[\text{yr}] \lesssim 3$ corresponding to phase I and phase III, respectively. They are characterized by $\tilde{T}_{b,p}/\tilde{T}_{b,e} < 1$. The maximum occurs over $0 \lesssim \log t[\text{yr}] \lesssim 1$ and is characterized by $\tilde{T}_{b,p}/\tilde{T}_{b,e} > 1$. Early ($\log t[\text{yr}] \lesssim -3$) and late ($\log t[\text{yr}] \gtrsim 4$) in the thermal

<footnote n="11">¹¹ The ratio of local temperatures is exactly the same, as the redshift factor $N/\Gamma$ (see Equation (7)) does not depend on $\theta$ along the isocontours of density (for barotropic EOSs this corresponds to the first integral of mechanical structure, derived from the conservation of energy-momentum).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Redshifted temperature at the bottom of the envelope as a function of polar angle $\theta$ at different instants, as shown in the legends. Results corresponding to the $M_B = 1.6 M_\odot$ NS built upon the APR EOS rotating at 716 Hz (1096 Hz) are illustrated in the left (right) panel. Note that the values on the Y-axis are normalized to the corresponding values at the pole ("p" subscript).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Ratio of polar to equatorial redshifted temperature as functions of time for NS with $M_B = 1.6 M_\odot$ rotating with different frequencies that are indicated in the legend. Left (right) panels correspond to the bottom (top) of the envelope (surface). Various EOS models, introduced in Table 1, are employed. The horizontal black-dotted line in the left panel marks the unity. The remaining horizontal dotted lines in the right panel correspond to the ratio $(g_{s,p}/g_{s,e})^{1/4}$; their color legend is the same as for the EOS models.</figcaption></figure>

evolution $\tilde{T}_b$ is the same at the pole and at the equator. At early moments, this reflects that the initial conditions are the same, see Section 2.3. The situation at late moments is due to the fact that the star has become isothermal.

The fact that, in the crust-core thermalization phase, $\tilde{T}_{b,p} < \tilde{T}_{b,e}$ has been previously noted by Miralles et al. (1993) and Schaab & Weigel (1998), who explained it in terms of the time needed for the cold wave traveling from the core to reach the crust and then punch through it. We remind that the core cools faster than the crust because it emits more neutrinos; in rotating stars the crust is thicker at the equator than anywhere else, see Table 2 and Figures 1 and 2.

Concerning the EOS dependence we note that, around the extrema, the ordering (in amplitude) of $\tilde{T}_{b,p}/\tilde{T}_{b,e}$ at 716 Hz is opposite to the one at 99% of Kepler frequency. Analysis of data in Table 2 reveals that the values of this ratio are correlated

with $\zeta_{cr}$. This finding highlights once more the role of the crust in the thermal evolution. We also note that, at 716 Hz and the near-Kepler limit, the maximum at $0 \lesssim \log t$ [yr] $\lesssim 1$ and the minimum at $2 \lesssim \log t$ [yr] $\lesssim 3$ are reached at the earliest (latest) times in NS built upon the APR (NL3-$\omega\rho$) EOS, for which the crust thickness is the smallest (largest), Table 2.

A comparison of the left and right panels allows one to see the influence of the heat-blanketing envelope and its dependence on the surface gravitational acceleration (Beznogov et al. 2021). For better visibility the ratio between surface gravitational acceleration at the pole and at the equator, to the power 1/4, is represented with dotted horizontal lines in the right panel; the same color convention as for $\tilde{T}_{s,p}/\tilde{T}_{s,e}$ is used. $(g_{s,p}/g_{s,e})^{1/4}$ depends on the oblateness of the star, and thus, on EOS stiffness; for correlation among them, see Table 2 and the
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The ratio of polar to equatorial redshifted temperature as a function of time at different densities inside the NS. The considered EOS model is APR. The star rotates at 1096 Hz. The considered values of density are indicated in the legend (in log). The curve marked "14.5" corresponds to the core while all the others correspond to the crust.</figcaption></figure>

comments in Section 3.2. The right panel in Figure 8 shows that, in most of the considered circumstances (i.e., EOS models and rotation frequencies) the redshifted surface temperature at the pole is higher than its counterpart at the equator; though, as shown by the APR and IUF curves at 716 Hz over $2 \lesssim \log t$ [yr] $\lesssim 2.5$ exceptions from this rule can occur.

Information complementary to that shown in Figure 8 (left panel) is provided in Figure 9, where the ratio between redshifted temperatures at the pole and at the equator is reported, as a function of time, for different densities inside the crust as well as one density in the core (labeled "14.5"). Only the case of the NS built upon APR and rotating at the (near) Kepler limit is considered. We note that the core is almost isothermal; the largest deviation from constant redshifted temperature is obtained at the bottom of the envelope; the lower the density, the more complex the temperature profile; high-density shells reach thermal equilibrium earlier than low-density shells. The minimum at $\log t$ [yr] $\approx -1.5$ is obtained only for $\rho = 10^{10} \text{ g cm}^{-3}$. Thus, it is likely some sort of reaction to the surface boundary condition. Taking into account that the minimum appears at the age of $t \ll 1$ yr, we speculate that it might not be physical. Indeed, the heat diffusion timescale through a heat-blanketing envelope with a bottom density of $\rho_b = 10^{10} \text{ g cm}^{-3}$ is a few months to a year (see, e.g., Beznogov et al. 2021) and variations of temperature happening on a shorter timescale are not reliable in such model.

Before concluding this section, let us comment on the sensitivity to the composition of the inner crust. A comparison of the cooling curves obtained using the D1M and HZD+NV EOSs for the crust indicates that only the crust-core thermalization phase (phase III) is affected. More precisely, our results show that when D1M is used, phase III starts by $\approx 20$ yr earlier.

## 5. Conclusions

We have investigated the standard cooling of isolated nucleonic NSs uniformly rotating with frequencies up to the mass shedding limit. The EOS dependence of thermal evolution was studied by considering different EOS models.

Our results indicate that due to rotation-induced deformation a heat blob is formed in the crust and persists around the equatorial plane over $1 \lesssim \log t$ [yr] $\lesssim 3$, which roughly corresponds to phases II and III. Its exact location and extension depend on the NS deformation, composition, and temperature profiles which, in turn, are determined by the core and crust EOSs, rotation frequency, and time in thermal evolution. As a consequence, the redshifted temperature manifests a complex distribution throughout the interior as well as at the bottom of the envelope. In particular, our results indicate that—depending on time—the maximum of $\tilde{T}_b$ can be obtained at the pole, at the equator or for intermediate values of polar angle. Large-scale structures in the temperature distribution can be explained in terms of the abovementioned heat blob and delays in cooling phases in equatorial regions compared to polar regions and static stars can be attributed to heat diffusion timescales. Both the heat blob and heat diffusion timescales depend on the specific heat, thermal conductivity, and geometry of the crust.

In contrast with the crust, the core is little affected by rotation. Indeed, even for near-Kepler models, the polar to the equatorial relative difference in temperature does not exceed a few percent. This situation can be explained by a combination of different factors: (i) being more compact, the core is less deformed; (ii) the core's heat conductivity is much higher than the crust's one, which means that the core requires less time in order to become isothermal; (iii) the temperature evolution in the core is driven primarily by neutrino emission from its whole volume, which is not affected by rotation.

The redshifted surface temperature $\tilde{T}_s(\theta)$, despite being determined by the redshifted bottom temperature $\tilde{T}_b$, might still show a different angular pattern. The explanation is that the surface gravity acceleration, which depends on the polar angle, enters the $T_s-T_b$ relation.

Our simple explanation in terms of the heat blob is unable to account for some finer details of the temperature distribution and evolution for near-Kepler models. Yet this is not a critical issue as studying long-term cooling of near-Kepler models without taking the spindown into account is a rather academic exercise. It is still useful as it allowed us to set limiting cases of the possible effects of rotation on long-term cooling.

Microphysics and observational issues not addressed in the present work, such as nucleonic pairing in the core and in the crust and the relation between observed and redshifted surface temperatures, will be considered elsewhere.

The authors thank D. Yakovlev for useful discussions. M.B. and A.R. acknowledge support from a grant from the Ministry of Research, Innovation and Digitization, CNCS/CCCDI—UEFISCDI, Project No. PN-III-P4-ID-PCE-2020-0293, within PNCDI III. D.P. and M.B. acknowledge financial support by the Mexican Consejo Nacional de Ciencia y Tecnología with a CB-2014-1 grant No. 240512 and the Universidad Nacional Autónoma de México through a UNAM-PAPIIT grant No. IN109520. M.B. also acknowledges support from a postdoctoral fellowship from UNAM-DGAPA. J.N. gratefully acknowledges the Italian Istituto Nazionale di Fisica Nucleare
(INFN), the French Center National de la Recherche Scientifique (CNRS), and the Netherlands Organization for Scientific Research for the construction and operation of the Virgo detector and the creation and support of the EGO Consortium.
*Software*:

1. NSCool 2D Rot (not publicly available)

2. LORENE (Gourgoulhon et al. 2016, https://lorene.obspm.fr)

3. Python (https://www.python.org) with packages:
   (a) ` matplotlib` (Hunter 2007, https:// matplotlib.org)
   (b) `f90wrap` (Kermode 2020, https://github.com/jameskermode/f90wrap)
   (c) ` numpy` (Harris et al. 2020, https:// numpy.org)
   (d) `scipy` (Virtanen et al. 2020, https://scipy.org)

# Appendix A
## Heat-blanketing Envelope and Long-term Thermal Evolution of NSs

In this paper, we have employed a widely accepted approach to handle long-term thermal evolution of NSs. The main idea was described in Section 1 and here we want to briefly address some finer points regarding the heat-blanketing envelope and surface boundary condition.

As mentioned in Section 1, the very outer layers of an NS are not degenerate enough and their mechanical structure depends noticeably on temperature. A conventional solution to circumvent this issue is to split the star into two domains: the star itself and a heat-blanketing envelope. This split significantly simplifies all calculations. This idea dates back to approximately the 1980s to pioneering works by Gudmundsson et al. (1982, 1983). The envelope goes from the actual surface of the star until its bottom at some density $\rho_b$, which is also called boundary density. The particular choice of $\rho_b$ depends on the specifics of the task at hand, like the required heat diffusion timescale through the envelope (see below), the presence of magnetic fields, etc. We adopted the *standard value* of $\rho_b = 10^{10} \text{ g cm}^{-3}$. Other choices are possible, but beyond the scope of the discussion here.

The envelope is thin, contains $\sim 10^{-7} M_\odot$ of matter, and has a negligible impact on the mechanical structure of the star. However, unlike the rest of the star, the thickness of the envelope depends strongly on the temperature and varies from $\sim 10$ m in cold stars to a few hundred meters in very hot stars with near-Eddington luminosity (Beznogov et al. 2020, 2021). Thus, to have a temperature-independent notion of the star's radius it is rather common to measure it from the center until $\rho_b$. The same reasoning applies to crust thickness. Since the envelope is thin, the change in radius is marginal; crust thickness is somewhat more affected. The values of radii and crust thicknesses presented in Table 2 are calculated this way.

While the envelope has a negligible impact on the mechanical structure of the star, it plays an important role in thermal evolution. Thus, the question of heat transport through the envelope should be addressed. Under the conditions discussed in, e.g., Gudmundsson et al. (1983), which are usually satisfied for long-term cooling, heat transport through the envelope can be computed *separately* in quasistationary plane-parallel approximation employing a temperature-dependent EOS of Coulomb plasma (Potekhin & Chabrier 2010). The result of such *separate* calculation is the relation between the temperature at the top of the envelope and the temperature at its bottom ($T_b$). Now, by

definition, the temperature at the top is the surface temperature ($T_s$). Thus, we have obtained the $T_s-T_b$ relation. Two comments are necessary here as they were mentioned in Section 4: (i) the functional dependence $T_s(T_b)$ (or its inverse $T_b(T_s)$) is, of course, instantaneous. One should keep in mind, however, that in reality heat needs some time to diffuse through the envelope. For $\rho_b = 10^{10} \text{ g cm}^{-3}$, the heat diffusion timescale is $\sim 1$ yr. Thus, one cannot study temperature variations occurring on a timescale faster than this; (ii) as shown in Gudmundsson et al. (1982, 1983) to a very high accuracy $T_b$ is not actually a function of $T_s$, but rather of $T_s/g_s^{1/4}$ ($g_s$ being the surface gravity acceleration). This is very convenient as it allows calculating the $T_s-T_b$ relation for some value of $g_s$ and then rescaling it for any other value. We have used this feature extensively, as for rotating star $g_s$ varies with polar angle $\theta$.

Now, from the point of view of the cooling computation, the outer boundary of the star is the bottom of the envelope, i.e., the outer boundary density is $\rho_b = 10^{10} \text{ g cm}^{-3}$. This is where we define the surface boundary condition of the main cooling calculation. And this is another reason to measure radii and crust thicknesses up until the outer boundary of the main calculation. To relate this outer boundary to the actual surface one employs the $T_s-T_b$ relation (see surface boundary condition equations in Section 2.3).

The last question to consider here is whether $\rho_b$ is high enough that thermal effects on the mechanical structure can be ignored at densities $\rho \ge \rho_b$. The answer is straightforward: Fermi temperature of electron gas at $\rho = \rho_b = 10^{10} \text{ g cm}^{-3}$ is $T_F \approx 10^{11} \text{ K}$ (for ions with mass-to-charge numbers ratio equal to 2). The contribution of thermal effects is of the order $(\pi T/T_F)^2$ (Baym 1981). Thus, even at the very beginning of the cooling calculation with the initial temperature at the bottom of the envelope set to $10^{10} \text{ K}$, this contribution is less than 10%. After a few more days the temperature drops more than twice and thermal contribution becomes totally negligible. Obviously, a few days are completely irrelevant for cooling on the timescales of years or more.

For more details on the envelope and how the $T_s-T_b$ relations are calculated, see the pioneering works by Gudmundsson et al. (1982, 1983) and/or the recent review by Beznogov et al. (2021) and references therein.

# Appendix B
NSCool 2D Rot

NSCool 2D Rot is an NS thermal evolution code developed by two of us (M.B. and D.P.). It is a major upgrade to the one-dimensional thermal evolution code NSCool by Page (2016). Our code is written fully in FORTRAN 90. All functions and subroutines *inherited* from NSCool were rewritten in FORTRAN 90. Our code makes extensive use of FORTRAN 90 derived types and employs parallelization via OpenMP in the main solver.

The numerical scheme is based on a finite difference approach. The spatial discretization uses a second-order central differences scheme with physical quantities (e.g., temperature, density, etc.) defined in the grid points, first spatial derivatives defined in the middle between the grid points and second spatial derivatives<fnref n="12" /> defined once again in the grid points; this

<footnote n="12">¹² Since thermal conductivity is assumed to be scalar in the present work, there are no mixed second derivatives as can be explicitly seen from Equations (11)–(13).</footnote>
approach is similar to the staggered grids used in computational fluid dynamics. Unlike other existing 2D NS thermal evaluation codes, the grid in NSCool 2D Rot explicitly includes points on the symmetry (i.e., rotation) axis as well as the grid point in the star center.

The temporal discretization is done with a fully implicit backward Euler scheme. Despite being only first order (in time), this method is L-stable (see, e.g., Hairer & Wanner 1996), which is crucial for solving the very stiff and highly nonlinear Equation (11). We have tested trapezoidal and midpoint (semi-)implicit A-stable second-order schemes and found that they produce small oscillations in the solution during photon cooling phase (which was not totally unexpected, see, e.g., Bui 1979). While some second or higher-order L-stable methods do exist (see, e.g., Bui 1979; Bulatov et al. 2011), their implementation is significantly more involved than for the backward Euler scheme. Taking into account that NSCool also relies on backward Euler temporal discretization and has proven to have no issues with accuracy, we decided that this method is sufficient for our purposes. Nevertheless, we have performed accuracy tests by comparing calculations done with different time steps.

Nonlinear algebraic equations obtained after spatial and temporal discretization are handled by our custom-written solver. It employs the Newton-Raphson method with a full Jacobian, which is computed by numerical differentiation via a second-order central difference scheme and taking into account its structure (i.e., only nonzero elements are computed and stored). The linearized equations at each Newton iteration are handled by our custom-written modification of generalized minimal residual solver (GMRES; Saad & Schultz 1986) with the tridiagonal part of the Jacobian as the preconditioner. This rather complicated setup with nested iterations (outer Newton iterations and inner GMRES iterations) was chosen for two main reasons: (i) the size of the Jacobian makes direct linear solvers too slow and its complicated structure prevents the employment of specialized banded or sparse matrix solvers (due to integration needed for the boundary condition at the origin, nine-diagonal structure of the Jacobian is supplemented by a column and a row) (ii) other solution strategies like simple fixed-point iterations (see, e.g., Hoffman & Frankel 2001) or the alternating-directions implicit method (see, e.g., Peaceman & Rachford 1955) were either not converging or converging very slowly during our initial tests.

Before using NSCool 2D Rot for practical applications, we performed extensive tests of the solver itself, checked the overall robustness of the spatial and temporal discretization, and compared the results with NSCool for spherically symmetric configurations.

To check the solver we exported the function that generates a set of discretized nonlinear algebraic equations at each time step to Python using `f90wrap` tool. Then we solved exported equations employing `scipy_optimize_root` at exactly the same time steps as in NSCool 2D Rot and compared the resulting 2D temperature distributions. We set the tolerances (maximum $\delta\bar{T}/\bar{T}$, where $\delta\bar{T}$ is a Newton-Raphson correction to $\bar{T}$) identically in both solvers (i.e., in Python and NSCool 2D Rot) and equal to $10^{-10}$. The maximum relative difference in $\bar{T}$ between two solvers solving exactly the same equations did not exceed $3 \times 10^{-10}$ for all time steps before the age of ~50 yr (further in time Python solver was not converging).

Then, we compared cooling curves computed by NSCool 2D Rot with different spatial grids and time steps. Doubling the number of grid points in both the $r$ and $\theta$ directions (i.e., 4 times more grid points in total) resulted in a maximum relative change in $\bar{T}_{s,e}$ of ~0.8%. This is on a par with the behavior of NSCool where doubling the number of grid points leads to an ~0.7% maximum relative change in surface temperature. Doubling the number of time steps resulted in an ~1.6% maximum relative change in $\bar{T}_{s,e}$. Again, this is on a par with NSCool and was expected as the requested temperature change between consecutive time steps<fnref n="13" /> was 5% and 2.5%. So it is not surprising that the two calculations differ from each other by a comparable amount.

Finally, we compared cooling curves produced by NSCool 2D Rot and NSCool for identical spherically symmetric systems. The maximum relative difference in $\bar{T}_s$ was ~6.4%. On the scale of, e.g., Figure 3 such difference would not be seen (the curves would be on top of each other). For any practical applications, this difference is irrelevant. Overall, we consider that the performance of the solver is excellent and employed spatial and temporal discretization scheme is on a par with NSCool.

The current implementation of NSCool 2D Rot was designed to work together with LORENE, but this is not mandatory as the codes do not communicate with each other directly. The results of LORENE (metric functions, matter distribution, etc.) are simply stored in files that are read by NSCool 2D Rot.

## ORCID iDs

Mikhail V. Beznogov https://orcid.org/0000-0002-7326-7270

Jérôme Novak https://orcid.org/0000-0002-6029-4712

Dany Page https://orcid.org/0000-0003-2498-4326

Adriana R. Raduta https://orcid.org/0000-0001-8421-2040

## References

* Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, PhRvL, 119, 161101

* Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, PhRvX, 9, 011001

* Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, PhRvC, 58, 1804

* Baade, W., & Zwicky, F. 1934, PNAS, 20, 254

* Baym, G. 1981, ApJ, 248, 767

* Beznogov, M. V., Page, D., & Ramirez-Ruiz, E. 2020, ApJ, 888, 97

* Beznogov, M. V., Potekhin, A. Y., & Yakovlev, D. G. 2021, PhR, 919, 1

* Birkl, R., Stergioulas, N., & Müller, E. 2011, PhRvD, 84, 023003

* Bonazzola, S., Gourgoulhon, E., Salgado, M., & Marck, J. A. 1993, A&A, 278, 421

* Bui, T. D. 1979, JACM, 26, 483

* Bulatov, M. V., Tygliyan, A. V., & Filippov, S. S. 2011, CMMPh, 51, 1167

* Burrows, A., & Lattimer, J. M. 1986, ApJ, 307, 178

* Camelio, G., Dietrich, T., Marques, M., & Rosswog, S. 2019, PhRvD, 100, 123001

* Canal, R., & Schatzman, E. 1976, A&A, 46, 229

* Degenaar, N., Page, D., van den Eijnden, J., et al. 2021, MNRAS, 508, 882

* Dommes, V. A., Gusakov, M. E., & Shternin, P. S. 2020, PhRvD, 101, 103020

* Dong, J. M., Wang, L. J., & Zuo, W. 2018, ApJ, 862, 67

* Fattoyev, F. J., Horowitz, C. J., Piekarewicz, J., & Shen, G. 2010, PhRvC, 82, 055803

* Fortin, M., Providência, C., Raduta, A. R., et al. 2016, PhRvC, 94, 035804

* Fortin, M., Raduta, A. R., Avancini, S., & Providência, C. 2021, PhRvD, 103, 083004

* Geppert, U., Küker, M., & Page, D. 2004, A&A, 426, 267

<footnote n="13">¹³ This user-defined parameter is a part of an adaptive time step selection algorithm.</footnote>
Gourgoulhon, E., Grandclément, P., Marck, J.-A., Novak, J., & Taniguchi, K. 2016, LORENE: Spectral methods differential equations solver, Astrophysics Source Code Library, ascl:1608.018

Gourgoulhon, E., Haensel, P., Livine, R., et al. 1999, A&A, 349, 851

Gudmundsson, E. H., Pethick, C. J., & Epstein, R. I. 1982, ApJL, 259, L19

Gudmundsson, E. H., Pethick, C. J., & Epstein, R. I. 1983, ApJ, 272, 286

Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Neutron Stars. 1, Equation of State and Structure (Astrophysics and Space Science Library), Vol. 326 (New York: Springer)

Haensel, P., Zdunik, J. L., & Dobaczewski, J. 1989, A&A, 222, 353

Hairer, E., & Wanner, G. 1996, Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems (2nd edn.; Berlin: Springer)

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hartle, J. B. 1967, ApJ, 150, 1005

Hartle, J. B., & Thorne, K. S. 1968, ApJ, 153, 807

Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Sci, 311, 1901

Hoffman, J., & Frankel, S. 2001, Numerical Methods for Engineers and Scientists (2nd edn.; Boca Raton, FL: Taylor & Francis)

Horowitz, C. J., & Piekarewicz, J. 2001, PhRvL, 86, 5647

Hunter, J. D. 2007, CSE, 9, 90

Kaminker, A. D., Kaurov, A. A., Potekhin, A. Y., & Yakovlev, D. G. 2014, MNRAS, 442, 3484

Kermode, J. R. 2020, JPCM, 32, 305901

Kluźniak, W., & Ruderman, M. 1998, ApJL, 505, L113

Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJL, 887, L24

Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJL, 918, L28

Miralles, J. A., van Riper, K. A., & Lattimer, J. M. 1993, ApJ, 407, 687

Mondal, C., Viñas, X., Centelles, M., & De, J. N. 2020, PhRvC, 102, 015802

Nagakura, H., Burrows, A., Radice, D., & Vartanyan, D. 2020, MNRAS, 492, 5764

Negele, J. W., & Vautherin, D. 1973, NuPhA, 207, 298

Negreiros, R., Schramm, S., & Weber, F. 2012, PhRvD, 85, 104019

Negreiros, R., Schramm, S., & Weber, F. 2013, PhLB, 718, 1176

Negreiros, R., Schramm, S., & Weber, F. 2017, A&A, 603, A44

Ootes, L. S., Page, D., Wijnands, R., & Degenaar, N. 2016, MNRAS, 461, 4400

Page, D. 2009, in Neutron Stars and Pulsars (Astrophysics and Space Science Library) ed. W. Becker, Vol. 357 (Berlin: Springer), 247

Page, D. 2016, NSCool: Neutron star cooling code, Astrophysics Source Code Library, ascl:1609.009

Page, D., Geppert, U., & Weber, F. 2006, NuPhA, 777, 497

Page, D., Homan, J., Nava-Callejas, M., et al. 2022, ApJ, 933, 216

Page, D., Lattimer, J. M., Prakash, M., & Steiner, A. W. 2004, ApJS, 155, 623

Page, D., Lattimer, J. M., Prakash, M., & Steiner, A. W. 2009, ApJ, 707, 1131

Page, D., Prakash, M., Lattimer, J. M., & Steiner, A. W. 2011, PhRvL, 106, 081101

Pais, H., & Providência, C. 2016, PhRvC, 94, 015808

Pascal, A., Novak, J., & Oertel, M. 2022, MNRAS, 511, 356

Peaceman, D. W., & Rachford, H. H., Jr. 1955, J. Soc. Indust. Appl. Math., 3, 28

Pons, J. A., Reddy, S., Prakash, M., Lattimer, J. M., & Miralles, J. A. 1999, ApJ, 513, 780

Potekhin, A. Y., & Chabrier, G. 2010, CoPP, 50, 82

Potekhin, A. Y., Chabrier, G., & Yakovlev, D. G. 1997, A&A, 323, 415

Potekhin, A. Y., Pons, J. A., & Page, D. 2015, SSRv, 191, 239

Raduta, A. R., Sedrakian, A., & Weber, F. 2018, MNRAS, 475, 4347

Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJL, 887, L21

Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJL, 918, L27

Roberts, L. F. 2012, ApJ, 755, 126

Saad, Y., & Schultz, M. H. 1986, SJSC, 7, 856

Schaab, C., & Weigel, M. K. 1998, A&A, 336, L13

Shternin, P. S., & Yakovlev, D. G. 2007, PhRvD, 75, 103004

Tsuruta, S. 2009, in Neutron Stars and Pulsars (Astrophysics and Space Science Library) ed. W. Becker, Vol. 357 (Berlin: Springer), 289

Vigano, D., Garcia-Garcia, A., Pons, J. A., Dehman, C., & Graber, V. 2021, CoPhC, 265, 108001

Viñas, X., Gonzalez-Boquera, C., Centelles, M., Mondal, C., & Robledo, L. M. 2021, Symm, 13, 1613

Vincent, F. H., Bejger, M., Rózańska, A., et al. 2018, ApJ, 855, 116

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Yakovlev, D. G., Kaminker, A. D., Gnedin, O. Y., & Haensel, P. 2001, PhR, 354, 1

Yakovlev, D. G., & Pethick, C. J. 2004, ARA&A, 42, 169

Yakovlev, D. G., & Urpin, V. A. 1980, SvA, 24, 303