© The Authors 2026

1. Introduction

From both the theoretical and experimental point of view, the properties of matter at very high densities, above the nuclear saturation density ρ₀ = 2.8 × 10¹⁴ g cm⁻³, are still unknown and remain as an open problem for fundamental physics. However, they can be studied through the thermal evolution of neutron stars (NSs) since the matter in their interior can reach several times ρ₀. By contrasting thermal evolution models with temperature measurements of NSs of different ages, we can learn, for example, about the composition and potential superfluid and superconducting states of very dense matter.

According to standard cooling models, NSs are born with a very high interior temperature, ∼10¹¹ K, and rapidly cool by emitting thermal neutrinos produced in their core. Neutrinos act as the main cooling agent and control the internal temperature for the first ∼10⁵ yr. In fact, several X-ray observations of young NSs have been made to contrast their temperatures with cooling curves and study the neutrino emission process, as well as the role of the composition, mass, and the state of matter of the NS (for a review see, e.g., Yakovlev & Pethick 2004).

Beyond the neutrino cooling age, photon emission from the NS surface becomes the main cooling mechanism. By ∼10⁷ yr, the temperature of passively cooling NSs is expected to drop below ∼10⁴ K, becoming undetectable by current telescopes. However, observations with the Hubble Space Telescope (HST) have revealed plausibly thermal ultraviolet (UV) radiation from the following objects:

• Millisecond pulsars (MSPs) PSR J0437−4715 (hereafter J0437; Kargaltsev et al. 2004; Durant et al. 2012) and PSR J2124−3358 (hereafter J2124; Rangelov et al. 2017), with characteristic ages τ ≡ P/(2Ṗ) > 10⁹ yr, where P is the rotation period and dots denote time derivatives, and

• Classical pulsars (CPs) PSR B0950+08 (hereafter B0950; Pavlov et al. 2017; Abramkin et al. 2022) and PSR J0108−1431 (hereafter J0108; Abramkin et al. 2021), with τ > 10⁷ yr.

Thus, unless the true ages of these pulsars are much less than their characteristic ages (which is very unlikely), this suggests the presence of internal heating mechanisms operating in very old NSs. Additionally, the nondetection of the CP PSR J2144−3933 (hereafter J2144; Guillot et al. 2019), which makes it the coolest NS known so far, also provides a constraint on potential heating mechanisms.

Gonzalez & Reisenegger (2010) analyzed several proposed NS heating mechanisms, namely rotochemical heating (Reisenegger 1995), vortex creep (Alpar et al. 1984), dark matter accretion (Kouvaris & Tinyakov 2010; de Lavallaz & Fairbairn 2010), crust cracking (Baym & Pines 1971), and magnetic field decay (Goldreich & Reisenegger 1992; Pons et al. 2025). They found that only rotochemical heating and vortex creep can plausibly explain the inferred surface temperature of J0437, the only pulsar mentioned above that had been detected by HST at that time (see also Fujiwara et al. 2024b).

For the younger, more strongly magnetized CPs, magnetic field dissipation might also give an important contribution. However, models for this process are not yet fully developed (e.g., Pons et al. 2025, and references therein), as the magnetic dissipation depends on many uncertain assumptions, such as the initial structure of the magnetic field (e.g., Becerra et al. 2022; Dehman et al. 2023), the impurity content of the crust (Cumming et al. 2004), and superfluid and superconducting properties of the core (Gusakov et al. 2020). Therefore, we do not include it in the present study.

Another potential heat source we do not consider here are relativistic particles from the magnetosphere hitting the NS surface (e.g., Harding & Muslimov 2002). As discussed by Kargaltsev et al. (2004), it is unlikely that a substantial fraction of this heat gets distributed over the whole NS surface in order to contribute significantly to the UV emission.

One of the processes we do consider is vortex creep (Alpar et al. 1984). It is based on the idea that neutrons can become superfluid in the deep layers of the NS crust, forming quantized vortices that dissipate heat by friction with the crustal nuclear lattice as they move outward as a result of the NS spin-down (for recent work on this process, see Fujiwara et al. 2024a; Nam & Sekizawa 2025).

Rotochemical heating is produced by non-equilibrium beta reactions in the NS core due to the departure from beta equilibrium produced by the continuous spin-down of the NS. Gonzalez & Reisenegger (2010) considered the case of rotochemical heating in NS cores with normal (non-Cooper-paired) matter, as in the study of Fernández & Reisenegger (2005). However, important theoretical progress has been made since then. Petrovich & Reisenegger (2010) modeled rotochemical heating in the presence of Cooper pairing, that is, with superfluid neutrons and/or superconducting protons, considering a uniform and isotropic Cooper pairing gap. The latter suppresses reactions, allowing chemical imbalances to grow more than in the case of normal matter, which eventually leads to a higher heat input for old fast rotating NSs (Reisenegger 1997). González-Jiménez et al. (2015) extended this study to the case of a nonuniform and anisotropic energy gap, finding that the evolution is similar to the case with a uniform gap equal to its minimum value. Yanagi et al. (2020) made it even more general by considering nonuniform gaps for both neutrons and protons. More recently, Kantor & Gusakov (2021) took into account that the proton superconductor in the NS core can be threaded by a magnetic field, either through quantized flux tubes with normal cores (type-II superconductor) or through normal domains of uncertain geometry (type-I superconductor). This allows for a higher rate of Urca reactions within these regions, as previously suggested by Schaab et al. (1998).

Furthermore, Gusakov et al. (2015) proposed a variant of rotochemical heating in which pulsars that have previously accreted a substantial fraction of their crust can be heated by pycnonuclear reactions in the deep layers of the crust. The continuous NS spin-down increases the departures from nuclear equilibrium, triggering nuclear reactions as the density increases. This mechanism only operates in MSPs and, as presented in the original paper, would produce luminosities comparable to those observed in J0437 and J2124. Recent insights into the equation of state of accreted crusts (Gusakov & Chugunov 2020) are likely to strongly reduce this effect. Since this reduction has not yet been calculated, we use the original prediction as an upper limit.

The goal of this work is to revisit the problem of heating mechanisms in old NSs. We calculate their thermal evolution curves including vortex creep, crustal heating, and rotochemical heating with and without Cooper pairing (assuming isotropic gaps and neglecting the effect of magnetic flux penetration in the proton superconductor) and compare them to the HST observations mentioned above.

Section 2 describes the observational data available on the five old pulsars used in this study. Sect. 3 presents a general discussion of the thermal evolution of NSs, whereas Sect. 4 we discuss the three internal heating mechanisms considered in the present work. Sect. 5 we discuss the numerical simulations and their results and compares them to the observational data, and Sect. 6 we present our conclusions. In Appendix A, we analyze thermal oscillations that occur when the energy gap is large and uniform and their suppression in the case of reduced reaction rates.

2. Observational data

We considered surface temperature measurements derived from HST observations of old CPs and MSPs, collected from previous studies, and corrected as described below. The measurements are summarized in Table 1. Our sample includes all pulsars with characteristic age of τ > 10⁷ yr that have available HST observations.

In general, the thermal emission from NSs is expected to be reprocessed by an atmosphere. However, for the relatively strong magnetic fields (B ∼ 10¹² G) and cool surface temperatures (T_s ∼ 10⁵ K) of old CPs, atmosphere models are rather uncertain (a condensed surface might even be present; see, e.g., Potekhin 2014). Nevertheless, state-of-the-art hydrogen atmosphere models generally follow a trend in which for increasingly strong magnetic fields, the UV tail of the spectrum at a given effective temperature approaches the Rayleigh–Jeans tail of a blackbody spectrum of the same temperature (e.g., Lloyd 2003)¹. In order to simplify the analysis of HST observations, we consistently considered surface temperatures of CPs derived from blackbody spectral fits.

We reanalyzed and updated the upper limit on the surface temperature of J2144, accounting for the latest recalibrations of HST’s Solar Blind Channel (SBC) instrument (Avila et al. 2019) in imaging mode as explained on the Space Telescope Science Institute webpage². Following Guillot et al. (2019), we performed blackbody fits to the photometric UV data, accounting for the uncertainties in the distance, dust extinction, and NS radius, and we considered the resulting temperature as an upper limit, since we did not include the potential contribution of a nonthermal component from the pulsar magnetosphere. The temperatures for B0950 and J0108 are used as reported by Abramkin et al. (2022) and Abramkin et al. (2021), respectively, who also fitted blackbody models and included the ACS/SBC recalibrations.

In the case of MSPs, their relatively low magnetic field, B ∼ 10⁸ G, greatly simplifies the modeling of the atmospheric thermal emission (see, e.g., Zavlin et al. 1996). For J0437, we consider spectral fits of non-magnetic hydrogen atmosphere models to the cold thermal component as in the recent analysis of González-Caniulef et al. (2019), but using UV spectroscopic data taken with the Space Telescope Imaging Spectrograph (STIS; used by Kargaltsev et al. 2004) instead of the SBC instrument (used by Durant et al. 2012), as the latter requires recalibration that is not yet available for the spectroscopic mode³. We also include the soft X-ray emission of J0437 as observed by ROSAT (see also Guillot et al. 2016). Including a Gaussian prior on the extinction and accounting for the effect of the hot polar caps, as described in González-Caniulef et al. (2019), we obtain the best fitting parameter values $T_s^{\infty }=2.{36}_{-0.14}^{+0.15}\times {10}^5\mathrm{K}$, $R=13.{00}_{-0.96}^{+1.16}\mathrm{km}$, N_H = (1.4 ± 0.3)×10²⁰ cm⁻², and E(B − V) = 0.01 ± 0.01, very similar to those reported in that paper. The resulting surface temperature is more tightly constrained than those obtained by Kargaltsev et al. (2004), Durant et al. (2012), and Bogdanov et al. (2019), which were used in previous comparisons with theoretical models (such as Fernández & Reisenegger 2005; Gonzalez & Reisenegger 2010; Petrovich & Reisenegger 2010; González-Jiménez et al. 2015; Kantor & Gusakov 2021). We note that Kargaltsev et al. (2004) and Durant et al. (2012) assumed blackbody models to infer the temperature, which do not appear to be consistent with realistic extinction values and are known to yield lower temperatures than hydrogen (or helium) atmosphere models (González-Caniulef et al. 2019). In addition, the NICER data analysis of Bogdanov et al. (2019) used a very crude background model, which could explain the difference in temperature deduced from the spectral analysis. Therefore, our new result might be more trustworthy than previous ones.

We also note that J0437 is the only pulsar in our sample that is in a binary system, which has allowed for a precise measurement of its mass, M = 1.418 ± 0.044 M_⊙ (Reardon et al. 2024), through radio timing, and its coordinate radius has been measured from its NICER X-ray lightcurve and spectrum to be $R=11.{36}_{-0.63}^{+0.95}\mathrm{km}$ (Choudhury et al. 2024), consistent with the value given above (see Hoogkamer et al. 2025 for a different analysis). Furthermore, the temperature of the white-dwarf companion implies an age 6.0 ± 0.5 Gyr (Durant et al. 2012), close to the pulsar’s characteristic age, τ = 6.6 Gyr, which means that, at the end of the mass transfer stage, the pulsar must have been spinning much faster than at present. Assuming dipole spin-down with a constant magnetic dipole moment (braking index n = 3), the initial period must have been in the range 0.7 < P₀/ms < 2.4 (at 68% confidence), compared to the current period P = 5.75 ms.

For J2124, the available data do not yet allow us to distinguish what fraction of the UV emission is thermal versus nonthermal power law, which only makes it possible to establish an upper limit on the surface temperature. Following Rangelov et al. (2017), we perform a blackbody fit to the recalibrated UV photometric data, accounting for the uncertainties associated with the distance and extinction and yielding a model-dependent upper limit T_s, BB^∞ < 1.8 × 10⁵ K. To this result, we applied the conversion formula to a hydrogen atmosphere provided by González-Caniulef et al. (2019) at the end of their Sect. 3, which yields an upward correction by a factor 2.7, and thus an absolute upper limit of T_s^∞ < 4.9 × 10⁵ K.

3. Thermal evolution models

The two classes of NSs studied in this paper (CPs and MSPs) have substantially different properties and evolutionary histories. CPs are born in the core collapse of a massive star, which leads to a supernova explosion and a remaining proto-NS with a very high initial core temperature T_c ∼ 10¹¹ K, a strong external magnetic dipole field B ∼ 10^{11 − 13} G (in addition to the likely presence of higher multipoles), and a moderate initial rotation period, P₀ ∼ 10 ms − 0.5 s (e.g., Faucher-Giguère & Kaspi 2006; Igoshev et al. 2022; Du et al. 2024). The magnetic dipole torque causes the star to progressively spin down, while the cooling and heating processes studied in the present paper cause the temperature to evolve over time. We assumed that the magnetic field is constant in time, ignoring the possibility of a substantial evolution during this phase (e.g., Goldreich & Reisenegger 1992).

If a NS is in a binary system, it might accrete ∼0.1 M_⊙ of matter (Özel et al. 2012; Antoniadis et al. 2016), heating it up to a core temperature T_c ∼ 10^8 − 9 K (Nava-Callejas et al. 2025) and reducing its rotation period to ∼1 − 10 ms (Sautron et al. 2026), thus turning it into a MSP. Either the accretion process itself (Bisnovatyi-Kogan & Komberg 1974) or the long intervening time (Cruces et al. 2019) also leads to a reduction of the magnetic dipole field to B ∼ 10^8 − 9 G, making the spin-down of MSPs much slower than that of CPs. In what follows, we consider the “birth” or initial state of a MSP to be the time when accretion stops and the NS evolution is no longer affected by its companion (which might by then have disappeared).

The thermal evolution of a NS depends on the reactions and the state of matter in the core, whose outer part is composed of degenerate neutrons (n), protons (p), electrons (e), and muons (μ). The composition of the inner core is uncertain, so we assumed the whole core to has this composition. The redshifted core temperature, T_c^∞, can be taken to be uniform because the typical thermal evolution time-scales are much longer than the diffusion time (Reisenegger 1995). Its evolution can be computed from (Thorne 1977):

$$\begin{array}{c}\hfill {\dot{T}}_c^{\infty }=\frac{1}{C}(L_H^{\infty }-L_{\nu }^{\infty }-L_{\gamma }^{\infty }),\end{array}$$(1)

where C is the total heat capacity of the star, $L_H^{\infty }$ is the total power injected by heating mechanisms, $L_{\nu }^{\infty }$ is the power emitted as thermal neutrinos and antineutrinos from the NS interior, and

$$\begin{array}{c}\hfill L_{\gamma }^{\infty }=4\pi R_{\infty }^2\sigma {(T_s^{\infty })}^4\end{array}$$(2)

is the power released as thermal photons from the surface. In Eq. (2), σ is the Stefan-Boltzmann constant, R_∞ is the NS radius, and T_s^∞ is its (effective) surface temperature. Throughout this paper, dots denote time-derivatives, and the superscript or subscript ∞ labels quantities redshifted to infinity, i.e., as measured far away from the NS. The temperatures of the core, T_c^∞, and the surface, T_s^∞, have a one-to-one relation to each other, which depends on the envelope composition (Gudmundsson et al. 1983; Potekhin et al. 1997). However, in the late-time evolution studied in this work, both T_c^∞ and T_s^∞ are low, which implies that the cooling in chemical equilibrium is strongly dominated by photon emission from the surface, which depends only on T_s^∞. Also, at late times, the heat released in the NS core can be considered to be emitted instantaneously⁴,

$$\begin{array}{c}\hfill L_{\gamma }^{\infty }\approx L_H^{\infty }-L_{\nu }^{\infty },\end{array}$$(3)

so the thermal evolution becomes nearly independent of the NS’s heat capacity and core temperature. Thus, the envelope model, although required as an input for the simulations, does not affect the late-time results for T_s^∞.

Two main kinds of weak interaction processes can contribute to the neutrino luminosity, L_ν^∞, and internal heating, L_H^∞, namely direct Urca (Durca) reactions,

$$\begin{array}{c}\hfill n\to p+\ell +{\overline{\nu }}_{\ell },p+\ell \to n+{\nu }_{\ell },\end{array}$$(4)

and modified Urca (Murca) reactions,

$$\begin{array}{c}\hfill n+N\to p+\ell +{\overline{\nu }}_{\ell }+N,p+\ell +N\to n+{\nu }_{\ell }+N,\end{array}$$(5)

where ℓ stands for a negatively charged lepton (electron or muon), and ν_ℓ and ${\overline{\nu }}_{\ell }$ for the corresponding neutrino or antineutrino. The latter reactions differ from the former by the presence of a spectator nucleon N (a neutron or proton), which does not change flavor, but is needed to satisfy the conservation of energy and momentum in the case of low abundance (and therefore small Fermi momentum) of protons and leptons. This allows the Murca reactions to be present in any NS core, whereas the much faster Durca reactions are only allowed above a certain threshold density. Thus, the latter are likely to be present only in the higher-density regions of relatively massive NSs. The luminosities due to both Durca and Murca reactions for normal matter in chemical equilibrium can be written as

$$\begin{array}{c}\hfill L_{\nu ,\mathrm{eq}}^{\infty }=({\stackrel{\sim }{L}}_e+{\stackrel{\sim }{L}}_{\mu }){(T_c^{\infty })}^q,\end{array}$$(6)

where q = 6 for Durca and q = 8 for Murca reactions (Yakovlev et al. 2001). The constants ${\stackrel{\sim }{L}}_{\ell }$ correspond to reactions involving electrons (ℓ = e) and muons (ℓ = μ) and are given in Table 2 for the particular models used in the present paper. For matter out of chemical equilibrium (discussed in Sect. 4.1) and/or with Cooper pairing, these luminosities must be multiplied by reduction factors that depend on the chemical imbalances and/or pairing gaps, respectively (e.g., Petrovich & Reisenegger 2010, and references therein).

Both CPs and MSPs start their evolution at a high temperature, so their cooling is initially dominated by neutrino emission (L_ν^∞), with photon emission (L_γ^∞) taking over ∼10⁵ yr later, at lower temperatures. The heating mechanisms contributing to $L_H^{\infty }$, which become relevant at relatively low temperatures, are explained in Sect. 4.

4. Heating mechanisms

4.1. Rotochemical heating

Since the birth of a NS, its continuous spin-down reduces the centrifugal force, increasing the local pressure and density everywhere inside the star. This changes the Fermi energies of the different particle species in the core, producing chemical imbalances,

$$\begin{array}{c}\hfill {\eta }_{\ell }\equiv {\mu }_n-{\mu }_p-{\mu }_{\ell },\end{array}$$(7)

where μ_i are the chemical potentials of each particle species (labeled as in Eqs. (4) and (5)), and inducing non-equilibrium Urca reactions that act in the sense of restoring the chemical equilibrium. As for the temperature, the fact that the diffusion time is much shorter than the evolution times of interest allows us to consider the redshifted chemical imbalances, η_ℓ^∞, to be uniform throughout the NS core (Reisenegger 1995). These variables evolve as

$$\begin{array}{cc}\hfill {\dot{\eta }}_e^{\infty }& =2W_{\mathit{npe}}\mathrm{\Omega }\dot{\mathrm{\Omega }}-Z_{\mathit{npe}}\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{\mathit{npe}}^{\infty }-Z_{\mathit{np}}\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{np\mu }^{\infty },\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill {\dot{\eta }}_{\mu }^{\infty }& =2W_{np\mu }\mathrm{\Omega }\dot{\mathrm{\Omega }}-Z_{\mathit{np}}\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{\mathit{npe}}^{\infty }-Z_{np\mu }\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{np\mu }^{\infty },\hfill \end{array}$$(9)

where Ω = 2π/P is the angular velocity of rotation of the NS, $\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{\mathit{npe}}^{\infty }$ and $\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{np\mu }^{\infty }$ are the total net rates at which neutrons are converted into protons and electrons and into protons and muons, respectively, and W_npe, W_npμ, Z_npe, Z_npμ, and Z_np are constants that depend on the structure of the NS (Fernández & Reisenegger 2005; Reisenegger et al. 2006; Petrovich & Reisenegger 2010) and whose values for the models studied in the present paper are given in Table 2⁵. The net effect of Urca reactions on the internal energy budget can be written as

$$\begin{array}{c}\hfill L_{\mathit{RH}}^{\infty }-L_{\nu }^{\infty }={\displaystyle \sum _{\ell =e,\mu }}{h}_{\ell }{\eta }_{\ell }^{\infty }\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{np\ell }^{\infty },\end{array}$$(10)

where h_ℓ are dimensionless functions of the temperature, the chemical imbalances, and the NS structure and composition. For normal matter in the limit |η_ℓ|/k_BT_c → ∞, where k_B is the Boltzmann constant, they satisfy h_ℓ → 1/2 for Durca reactions and h_ℓ → 5/8 for Murca reactions, representing the fraction of the chemical energy release that gets deposited as heat in the NS core, while the rest is emitted as neutrinos and antineutrinos, which escape. For |η_ℓ|/k_BT_c ≲ 5, one obtains h_ℓ < 0, corresponding to net cooling (Fernández & Reisenegger 2005).

Since the spin-down depends on the magnetic field, CPs increase their chemical imbalances much faster than MSPs (Gonzalez & Reisenegger 2010). For the strong magnetic field of CPs, the effect of rotochemical heating on temperature appears early (at t ≲ 10⁵ yr), while for the weaker magnetic field of MSPs it appears later (t ≳ 10^6 − 7 yr).

At later stages, the impact of compression due to spin-down on CPs becomes negligible, and their chemical imbalances only decrease slowly due to weak interactions, injecting heat and thus slowing down the cooling. MSPs, on the other hand, might reach a quasi-steady state where reactions counteract the rise in chemical imbalances caused by spin-down:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{np\ell ,\mathrm{qs}}^{\infty }\approx -K_{\ell }\mathrm{\Omega }\dot{\mathrm{\Omega }},\end{array}$$(11)

where

$$\begin{array}{cc}\hfill K_e& \equiv -\frac{2(Z_{np\mu }{W}_{\mathit{npe}}-Z_{\mathit{np}}{W}_{np\mu })}{Z_{\mathit{npe}}{Z}_{np\mu }-Z_{\mathit{np}}^2},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill K_{\mu }& \equiv -\frac{2(Z_{\mathit{npe}}{W}_{np\mu }-Z_{\mathit{np}}{W}_{\mathit{npe}})}{Z_{\mathit{npe}}{Z}_{np\mu }-Z_{\mathit{np}}^2},\hfill \end{array}$$(13)

while heating and cooling effects also balance each other, yielding

$$\begin{array}{c}\hfill L_{\gamma ,\mathrm{qs}}^{\infty }\approx -K\mathrm{\Omega }\dot{\mathrm{\Omega }}{\overline{\eta }}_{\mathrm{qs}}^{\infty },\end{array}$$(14)

where K ≡ K_e + K_μ and

$$\begin{array}{c}\hfill {\overline{\eta }}_{\mathrm{qs}}^{\infty }\equiv \frac{1}{K}{\displaystyle \sum _{\ell =e,\mu }}{K}_{\ell }{h}_{\ell }{\eta }_{\ell ,\mathrm{qs}}^{\infty }\end{array}$$(15)

(Fernández & Reisenegger 2005; Petrovich & Reisenegger 2010). Here and below, the subscript “qs” labels variables evaluated in the quasi-steady state. Thus, the fraction of the spin-down power emitted as thermal photons is

$$\begin{array}{c}\hfill \frac{L_{\gamma ,\mathrm{qs}}^{\infty }}{|{\dot{E}}_{\mathrm{rot}}|}=\frac{K{\overline{\eta }}_{\mathrm{qs}}^{\infty }}{I}\sim 3\times {10}^{-4}\frac{{\overline{\eta }}_{\mathrm{qs}}^{\infty }}{\mathrm{MeV}},\end{array}$$(16)

where E_rot = IΩ²/2 is the kinetic energy of rotation and I is the moment of inertia of the star.

For a given NS structure with given weak interaction processes (Murca, Durca, or others, with or without Cooper pairing) and given NS spin-down parameters at a certain time, Eqs. (11) and (14) determine the three quasi-steady-state parameters η_e, qs^∞, η_μ, qs^∞ and T_s, qs^∞. In the relevant limit k_BT_c^∞ ≪ |η_ℓ^∞|, the reaction rates become independent of T_c^∞, so Eq. (11) determines η_e, qs^∞ and η_μ, qs^∞. Replacing these in Eq. (14) yields the surface temperature T_s, qs^∞, independent of the NS envelope composition.

Conversely, if a NS is in this quasi-steady state, the measured values of its spin-down parameters Ω and $\dot{\mathrm{\Omega }}$ and its photon luminosity can be used in Eq. (14) to infer ${\overline{\eta }}_{\mathrm{qs}}^{\infty }$, without having to assume specific weak interaction processes. For the measured spin parameters and surface temperature of J0437 and the parameter K of our Model M, one obtains ${\overline{\eta }}_{\mathrm{qs}}^{\infty }\approx 1.5\mathrm{MeV}$, but we note that K varies by a factor ∼3 between different equations of state, and the dependence on (T_s^∞)⁴ introduces an additional large uncertainty.

In order to reach the quasi-steady state from an initial chemical equilibrium state through the spin-down of the star (represented by the first term on the right-hand side of Eqs. (8) and (9)), the initial rotation rate Ω₀ must satisfy

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_0^2-{\mathrm{\Omega }}^2\gtrsim \frac{K{\overline{\eta }}_{\mathrm{qs}}^{\infty }}{{\displaystyle {\sum }_{\ell }}(-W_{np\ell })K_{\ell }{h}_{\ell }}\approx \frac{L_{\gamma ,\mathrm{qs}}^{\infty }}{\mathrm{\Omega }\dot{\mathrm{\Omega }}{\displaystyle {\sum }_{\ell }}{W}_{np\ell }{K}_{\ell }{h}_{\ell }}·\end{array}$$(17)

Taking again the observed parameters of J0437 and the theoretical parameters from our Model M, and assuming h_e ≈ h_μ ≈ 1, we obtain that the initial rotation period of this pulsar must satisfy P₀ ≲ 1.6 ms, substantially shorter than its present period P = 5.75 ms, but consistent with the range inferred from the age of its white dwarf companion (see Sect. 2). In any case, we reiterate that there are theoretical and observational uncertainties affecting these estimates.

When |η_ℓ^∞|≫k_BT_c^∞ and |η_ℓ^∞|≫Δ_j^∞ for both ℓ = e, μ and j = n, p (which includes the case of normal matter out of chemical equilibrium, with Δ_j^∞ = 0), the net reaction rates become (Reisenegger 1995; Petrovich & Reisenegger 2010)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{D,np\ell }^{\infty }& \approx \frac{42}{457}\frac{{\stackrel{\sim }{L}}_{D,\ell }}{{(\pi k_B)}^6}{({\eta }_{\ell }^{\infty })}^5\mathrm{for}\mathrm{Durca}\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Gamma }}}_{M,np\ell }^{\infty }& \approx \frac{24}{11513}\frac{{\stackrel{\sim }{L}}_{M,\ell }}{{(\pi k_B)}^8}{({\eta }_{\ell }^{\infty })}^7\mathrm{for}\mathrm{Murca},\hfill \end{array}$$(19)

where ${\stackrel{\sim }{L}}_{D,\ell }$, ${\stackrel{\sim }{L}}_{M,\ell }$, are the constants appearing in the neutrino luminosity for each process in chemical equilibrium, as given in Eq. (6). Setting the net rates from Eqs. (18) and (19) equal to the quasi-steady-state rate for J0437 from Eq. (11) yields chemical imbalances substantially smaller than those required to account for the thermal photon luminosity of this pulsar.

Thus, to account for the observed luminosity of J0437 by rotochemical heating, at least one of the nucleon species in the NS core must be Cooper-paired with a relatively large energy gap. Petrovich & Reisenegger (2010) included the effects of Cooper pairing with uniform energy gaps in the thermal and chemical evolution equations of Fernández & Reisenegger (2005) (slightly corrected by Reisenegger et al. 2006) by introducing reduction factors, which are functions of the chemical imbalances and the temperature of the core. This work confirms the arguments of Reisenegger (1997) in the sense that in the Cooper-paired state, reactions are strongly suppressed as long as the chemical imbalances remain below the threshold energy

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{thr}}=\{\begin{array}{cc}{\mathrm{\Delta }}_n+{\mathrm{\Delta }}_p\hfill & \mathrm{for}\mathrm{Durca}\mathrm{reactions},\hfill \\ min\{{\mathrm{\Delta }}_n+3{\mathrm{\Delta }}_p,3{\mathrm{\Delta }}_n+{\mathrm{\Delta }}_p\}\hfill & \mathrm{for}\mathrm{Murca}\mathrm{reactions},\hfill \end{array}\end{array}$$(20)

where Δ_n and Δ_p are the energy gaps of the neutrons and protons. As a result, the chemical imbalances can grow freely up to this threshold value. Once it is reached, Urca reactions become activated. These reactions stop further growth, rapidly heat the star, and keep it in a quasi-steady state, as in the case of normal matter. Therefore, the chemical imbalances remain near Δ_thr, i.e., at a value that can be much higher than in the normal case, releasing more energy in each reaction and thus leading to a higher quasi-steady temperature, given by Eqs. (14) and (15) with h_ℓ ≈ 1 and ${\eta }_{\ell ,\mathrm{qs}}^{\infty }\approx {\mathrm{\Delta }}_{\mathrm{thr}}^{\infty }\approx {\overline{\eta }}_{\mathrm{qs}}^{\infty }$. This hypothesis can match the parameters of J0437 if Δ_thr^∞ ≈ 1.5 MeV. (The value obtained by Petrovich & Reisenegger 2010 is lower because they used a lower surface temperature for this pulsar.)

González-Jiménez et al. (2015) extended the previous studies by considering the effect of density dependence and anisotropy of the energy gaps. They found that, for the case of the Murca process, the star reaches the quasi-steady state when the chemical imbalances exceed the smallest value of Δ_thr within the core, making the behavior quite similar to the case with uniform gaps.

On the other hand, Petrovich & Reisenegger (2011) found that for Durca reactions with a large uniform and isotropic energy gap, MSPs do not settle to a quasi-steady state, and instead their chemical imbalances and temperature oscillate strongly (on cooling time scales). Thus, in this case, it is not possible to make a reliable prediction about the temperature of these stars or to infer conclusions about the heating mechanisms from the observed temperature. We note, however, that the gaps are not expected to be uniform, but to vary with density. Therefore, the chemical imbalances will grow to the point where they reach approximately the minimum value of Δ_thr, activating reactions in a thin shell where the gap takes this value, thus impeding further growth (González-Jiménez et al. 2015). This reduces the effective volume of the core where reactions are active. We model this effect in Appendix A, where we show that oscillations also happen for Murca reactions with a large enough, uniform energy gap. However, for both Durca and Murca, reducing the reaction rates by a small enough factor f, the oscillations are eliminated. In the simulations reported in the rest of this paper, we introduce such a factor to mimic the effect of a nonuniform gap and eliminate the oscillations.

We note that none of this work considers lepton decay (conversion of electrons into muons, and vice-versa; Alford & Good 2010), a process that is intrinsically much slower than the Urca processes, but that is not suppressed by Cooper pairing, and which tends to reduce the difference between the electron and muon chemical potentials. Because Urca processes make these two chemical potentials evolve in a very similar way, their differences generally remain small anyway. Thus, this process is not expected to substantially change the results (Kantor & Gusakov 2021) and is also ignored in the present work.

Another effect not considered in the work described above, but initially pointed out by Schaab et al. (1998) and recently discussed and evaluated by Kantor & Gusakov (2021), is that, when the protons undergo the transition to the superconducting state, the magnetic flux gets concentrated into small regions in which the protons continue behaving as normal (nonsuperconducting) particles. Since this introduces an additional parameter, namely the fraction of the total volume in these regions, we do not consider it here (effectively assuming perfect flux expulsion from the superconductor), but we expect to address it in a forthcoming paper.

4.2. Vortex Creep

In the inner crust, the superfluid neutrons form an array of discrete quantum vortices, which are pinned to the nuclei of the solid nuclear lattice, generally keeping the superfluid rotating faster than the rest of the star. As the star spins down, the relative velocity between the superfluid and the vortex array grows, inducing a Magnus force. When a critical velocity difference is reached between the superfluid and the vortex array, the Magnus force overcomes the pinning force, and the pinning and unpinning of many vortices releases energy that heats the star, with a total time-averaged power (Alpar et al. 1984)

$$\begin{array}{c}\hfill L_{\mathrm{vc}}=J|\dot{\mathrm{\Omega }}|,\end{array}$$(21)

where J is the excess angular momentum, given by $J\approx I_p{\stackrel{\sim }{\omega }}_{\mathit{cr}}$, I_p is the moment of inertia of the superfluid in the pinning layer, and ${\stackrel{\sim }{\omega }}_{\mathit{cr}}={({\mathrm{\Omega }}_s-{\mathrm{\Omega }}_c)}_{\mathit{cr}}$ is the average over this layer of the critical lag between the angular velocity of the crust, Ω_c, and the macroscopic angular velocity of the superfluid, Ω_s, and $\dot{\mathrm{\Omega }}$ is the time-derivative of the star’s angular velocity (almost the same for the crust and the superfluid). The value of J is expected to be in the range ∼10^{43 − 45} erg s, so

$$\begin{array}{c}\hfill L_{\mathrm{vc}}\sim ({10}^{29}-{10}^{31})|{\dot{\mathrm{\Omega }}}_{-14}|{\mathrm{erg}\mathrm{s}}^{-1},\end{array}$$(22)

with ${\dot{\mathrm{\Omega }}}_{-14}=\dot{\mathrm{\Omega }}/({10}^{-14}{\mathrm{s}}^{-2})$ (Gonzalez & Reisenegger 2010). If, at a late time t ≳ 10⁶ yr in the evolution of a NS, vortex creep is the dominant heating mechanism, its photon luminosity will be given by Eqs. (21) and (22). We note that this was used by Fujiwara et al. (2024a) to obtain values of J for many pulsars. However, some of those values differ by orders of magnitude, which we do not expect from NS models, making it likely that different heating mechanisms are at work.

4.3. Crustal heating

When NSs accrete matter from a binary companion, their crust gets compressed, inducing some nuclear transformations and leaving some layers of matter close to the threshold for such transformations to occur. After accretion has stopped, the decrease of the centrifugal force due to the continuous spin-down of the star causes an additional small compression that can make the matter cross this threshold, inducing further exothermic reactions (electron captures, neutron emissions, or pycnonuclear reactions) at the interfaces between the deep layers of the crust with different dominant nuclei, releasing heat at a rate of Gusakov et al. (2015)

$$\begin{array}{c}\hfill L_{\mathrm{CH}}\approx 7\times {10}^{28}a\frac{R_6^7}{M_{1.4}^2}{\mathrm{\Omega }}_3|{\dot{\mathrm{\Omega }}}_{-14}|{\displaystyle \sum _i}\frac{p_i}{{10}^{31}{\mathrm{erg}\mathrm{cm}}^{-3}}\frac{q_i}{1\mathrm{MeV}}{\mathrm{erg}\mathrm{s}}^{-1}.\end{array}$$(23)

Here R₆ = R/(10⁶ cm), R is the stellar radius, M_1.4 = M/(1.4 M_⊙), M is the stellar mass, Ω₃ = Ω/(10³ s⁻¹), Ω is the angular velocity, and a is a positive, dimensionless parameter of order unity. The various transitions between nuclear species in the crust of the NS, labeled by index i, occur at critical pressures p_i and can release an energy q_i per nucleon, whose respective values are listed in Haensel & Zdunik (1990). We considered all reactions to determine the upper limit of the temperatures that this mechanism can predict, resulting in the following expression:

$$\begin{array}{c}\hfill L_{\mathrm{CH}}\approx 1.32\times {10}^{30}{\mathrm{\Omega }}_3|{\dot{\mathrm{\Omega }}}_{-14}|{\mathrm{erg}\mathrm{s}}^{-1}.\end{array}$$(24)

As mentioned in Sect. 1, Gusakov & Chugunov (2020) suggested that this effect might be much weaker. Therefore, we regard Eq. (23) as an upper limit. We consider this mechanism only for MSPs, which are generally accepted to have been recycled by accretion.

5. Thermal evolution simulations for the observed NSs

5.1. Assumptions and setup

To identify and test the dominant heating mechanisms, we compared the predictions of thermal evolution models to the HST observations. The thermal evolution curves were calculated by integrating Eqs. (1), (8), and (9) with a modified version of the code of Petrovich & Reisenegger (2010), under the assumption of magnetic dipole spin-down with a constant magnetic field strength B = 3.2 × 10¹⁹(P[s]Ṗ)^1/2 G inferred from the observed values of the period P and its time-derivative Ṗ for the respective pulsar. The initial period was set to P₀ = 1 ms for MSPs and P₀ = 5 ms for CPs. For CPs we assumed a fully catalyzed envelope containing heavy elements (Gudmundsson et al. 1983), whereas for MSPs we assumed an accreted envelope with light elements (Potekhin et al. 1997), but we remind the reader that the late-time evolution of T_s^∞ does not depend on this choice.

For the CPs (B0950, J0108, and J2144), we assumed an initial core temperature T₀^∞ = 10¹¹ K for all curves. For the MSPs (J0437 and J2124) we set T₀^∞ = 10⁹ K, roughly matching the conditions expected at the end of the accretion phase. We verified that, for both cases, any choice of T₀ ≳ 10⁹ K does not affect our results, since the initial neutrino emission is so strong that it causes the temperature to drop to 10⁹ K in less than one year. When considering rotochemical heating, we set the initial values for the chemical imbalances equal to zero, which is irrelevant to our results as, for the magnetic field strengths considered here, evolution curves for different initial values rapidly converge (Reisenegger 1995; Fernández & Reisenegger 2005; Kantor & Gusakov 2021).

5.2. Results of simulations and comparison to observations

In Fig. 1, we show thermal evolution simulations for each of the observed pulsars, considering different internal heating mechanisms (or none) and the normal or Cooper-paired states of the nucleons in the core, as summarized in Table 3. The CPs B0950 and J0108 were combined in the same panel, because they have very similar magnetic fields, and therefore, as long as they had similar initial rotation periods and internal properties, they could be regarded as two different points along the same thermal evolution curves, for whose calculation we considered the magnetic field value inferred for B0950. We plotted the data point for the measured temperature or upper limit of each pulsar at the time t at which the computed spin-down parameters P(t) and Ṗ(t) match the observed values. This time is in general somewhat shorter than the characteristic age, τ ≡ P/(2Ṗ), which would agree in the limit of infinitely rapid initial rotation.

Fig. 1. Surface temperature evolution curves for each observed pulsar considering models with different reheating mechanisms. In each panel, all curves were calculated assuming magnetic dipole spin-down with a constant dipole moment inferred from the measured spin-down parameters P and Ṗ of the respective pulsars. B0950 and J0108 are plotted together, because they have nearly the same inferred dipole field. Their curves were generated with the magnetic field of B0950. For the MSPs (J0437 and J2124), the initial conditions were taken as P0 = 1 ms and T0∞ = 109 K, whereas for the CPs (B0950, J0108, and J2144) we assumed P0 = 5 ms and T0∞ = 1011 K. The (1σ) error bars or upper limits for the pulsars are placed at the times at which the present spin parameters are reached. The models used for Murca and Durca reactions are given in Table 2 as Model M and Model D, respectively. The labels correspond to the following cases (see Table 3 for a schematic summary): I Passive cooling with Murca reactions in normal (nonsuperfluid and nonsuperconducting) matter, solid black line. II Passive cooling with Durca reactions in normal matter, dashed gray line. III Rotochemical heating with Murca reactions in normal matter, solid maroon line. IV Rotochemical heating with Durca reactions in normal matter, dashed maroon line. V Vortex creep with Murca reactions setting J = 3 × 1043 erg s and assuming normal particles in the core, solid green line. VI Same as V, but with Durca reactions, dashed green line. VII Crustal heating (only present in MSPs) with Murca reactions and normal particles in the core, solid cyan line. VIII Same as VII, but with Durca reactions, dashed cyan line. IX Rotochemical heating with Murca reactions considering normal protons and superfluid neutrons with a uniform energy gap Δn = 1.5 MeV, with reaction rates reduced by a factor f = 10−2, solid orange line. X Same as IX, but with Durca reactions and reaction rates reduced by a factor f = 10−8, dashed orange line. XI Rotochemical heating with Murca reactions considering normal protons and superfluid neutrons with Δn = 1.5 MeV and reaction rates reduced by a factor f = 10−2, as well as vortex creep with J = 3 × 1043 erg s, magenta line with star markers.

One can see that the curves without reheating (labeled I and II) drop to very low temperatures, T_s^∞ ∼ 10⁴ K at early times t < 10⁷ yr, so they cannot match the relatively high temperatures T_s^∞ ∼ 10⁵ K, of the observed pulsars. This is shown here for matter that is not Cooper paired, but remains true for any plausible assumptions regarding Cooper pairing of neutrons and/or protons (Page et al. 2013). Thus, some heating mechanism is required to explain the observations.

For rotochemical heating in CPs containing normal matter, which is discussed and illustrated in some detail in Gonzalez & Reisenegger (2010) and labeled III and IV in the lower panels of Fig. 1, neutrino emission is the primary cooling mechanism from the NS’s birth until t ∼ 10⁴ − 10⁵ years. During this period, the chemical imbalances (shown in Fig. 3 of Gonzalez & Reisenegger 2010) increase as a result of the star’s spin-down, reaching a maximum that largely depends on the initial period and the strength of the external dipole field. Beyond the neutrino cooling age, the NS cools down through thermal photon emission from its surface. At this stage, the effect of spin-down becomes negligible, but the slow release of the previously accumulated chemical imbalances makes the temperature decrease much more slowly than in the absence of heating. We have verified that curves with different initial periods asymptotically converge, because the late evolution of the chemical imbalances in CPs is fully determined by the beta reactions, while the temperature is determined by the balance of the heat input through these reactions against the cooling through photons from the NS surface (Eq. (3)).

For rotochemical heating in MSPs with normal matter, previously studied by Reisenegger (1995), Fernández & Reisenegger (2005), and Gonzalez & Reisenegger (2010) and labeled III and IV in the upper panels of Fig. 1, the thermal evolution is slightly different. From the end of recycling (t = 0) until t ∼ 10⁴ − 10⁵ yr, neutrino emission is the dominant cooling mechanism, and the chemical imbalances increase due to spin-down. However, the latter happens more slowly than in CPs, because of the weaker dipole field, so the imbalances do not saturate during this period. As the star transitions to the photon cooling age, the chemical imbalances continue to grow until reaching a quasi-steady state in t ∼ 10⁶ − 10⁷ yr. At this stage, the spin-down remains important, but it is balanced by the effect of the reactions (Eq. (11)), so the NS reaches the quasi-steady state discussed in Sect. 4.1, and the chemical imbalances remain roughly constant. Similarly, the heat input from the reactions is again balanced by the cooling effect through thermal photon emission from the NS surface (Eq. (3)).

Rotochemical heating with Murca reactions in normal matter (III) can account for the observed temperatures of the CPs B0950 and J0108, and is consistent with the upper limits for the CP J2144 and the MSP J2124. However, as argued in Sect. 4.1, it is insufficient to explain the high temperature of MSP J0437. Rotochemical heating with Durca reactions (IV), whose physics in principle is quite similar to the Murca scenario, produces faster beta reactions in normal matter, which leads to the saturation of the chemical imbalances at lower values and, consequently, cooler temperatures compared to Murca reactions. Any choice of initial period, initial chemical imbalances, and NS model for Durca reactions in normal matter is found to underestimate the observed temperatures.

For rotochemical heating in a NS core with superfluid neutrons, superconducting protons (with no magnetic flux penetrating them), or both, illustrated by curves IX and X in all panels, beta reactions are suppressed as long as the chemical imbalances remain below the threshold Δ_thr given by Eq. (20). As explained in Sect. 4.1, this threshold must be rather high, Δ_thr ∼ 1.5 MeV, to explain the thermal photon emission of J0437. For the imbalances to reach this high threshold, the initial rotation of the pulsar must be fast, P₀ ≲ 2 ms, which is plausible for MSPs, but not for CPs. Since, in our simulations, we chose P₀ = 1 ms for the MSPs, P₀ = 5 ms for the CPs, and Δ_thr = 1.5 MeV for all, we find that the MSPs reach the threshold, so rotochemical heating is activated, explaining the thermal emission of J0437, whereas CPs do not reach the threshold and therefore are not reheated, cooling down to low temperatures⁶. Thus, it is not possible to explain all the observations with rotochemical heating as the only heating mechanism, unless a small amount of magnetic flux is present inside a core with superconducting protons, which we expect to study in a forthcoming paper. We note that in the simulations with Cooper pairing we assumed uniform and isotropic gaps, but reduced the reaction rates by constant factors f = 10⁻² for Murca and f = 10⁻⁸ for Durca reactions, which are roughly the maximum values that suppress the oscillations studied in Petrovich & Reisenegger (2011) and in Appendix A. Any smaller values do not affect the results significantly. These reductions somehow mimic the effect of nonuniform energy gaps, which would only allow reactions in some thin shell within the NS core where Δ_thr is smallest (González-Jiménez et al. 2015).

Another potentially important heating mechanism is vortex creep, represented as curves V and VI in all panels of Fig. 1. This mechanism leads to significant deviations from standard cooling models once a star reaches the photon cooling age in CPs and somewhat later in MSPs. The heat generated by vortex friction in the crust is proportional to the star’s angular deceleration $|\dot{\mathrm{\Omega }}|$, which is generally larger in CPs than in MSPs.

In our analysis of vortex creep, we treat the uncertain excess angular momentum, J, as a free parameter, which is constrained by the upper limit for J2144, corresponding to J ≤ 3 × 10⁴³ erg s. Adopting this upper bound for both of our vortex creep models (V and VI), it can account for the temperature of B0950 but is insufficient for J0437. Thus, if J is approximately universal for all NSs (as we expect since it is determined by the structure of the NS crust), vortex creep cannot by itself account for all the observations.

Crustal heating (labeled VII and VIII in the upper panels of Fig. 1) occurs only in recycled pulsars, that is, MSPs. Similarly to vortex creep, the quasi-steady temperature in this case does not depend on the type of Urca reactions that occur in the core, but is solely determined by the pycnonuclear reactions in the deeper layers of the crust and the current spin-down of the NS, quantified by $\mathrm{\Omega }\dot{\mathrm{\Omega }}$. We included all possible pycnonuclear reactions listed in Haensel & Zdunik (1990) to establish an upper limit for this mechanism. Both MSPs achieve very similar temperatures because they have comparable spin-down rates. However, this upper limit is too low to account for the observed temperature of J0437.

In both vortex creep and crustal heating, the late-time temperature is set by a balance between the respective heating effect and surface photon emission. Thus, it depends only on the current spin parameters of the pulsar and is unaffected by the type of Urca reactions and the Cooper pairing in the core.

Since rotochemical heating with Cooper-paired nucleons with a large gap can account for the surface temperature of J0437 and vortex creep for that of B0950, but not vice-versa, we also model the combination of both effects (each with the same parameters as before), which is labeled XI in all panels. It can be seen that, at any given time, the temperature is roughly the higher one among the two obtained by considering each effect independently. At late times, rotochemical heating dominates for MSPs, explaining the temperature of J0437, and vortex creep dominates for CPs, explaining that of B0950.

It is interesting to note that this model also predicts that the temperatures of the other pulsars should be close to the upper limits inferred from their HST observations. Based on this fact, we speculate that the emission detected for J2124 and possibly detected for J0108 should have an important thermal contribution, and the thermal emission of J2144 should not be much lower than the upper limit from its nondetection.

5.3. Sensitivity to initial conditions

As mentioned in Sect. 3, NSs are born with very high core temperatures, T_c ∼ 10¹¹ K. Although not quite as extreme, MSPs also emerge from their accretion phase with a high core temperature, T_c ∼ 10^8 − 9 K. In the early neutrino cooling stage, for T_c ≳ 10⁸ K, the luminosity is a strong function of temperature, ∝T_c⁶ for Durca and ∝T_c⁸ for Murca reactions (assuming normal matter), so the NS quickly loses any memory of its initial temperature. The pulsars included in the present study have surface temperatures T_s ≲ 2 × 10⁵ K, corresponding to core temperatures of at most a few times 10⁶ K, so their present state should be independent of their initial temperature.

Similarly, it has been pointed out that, at such high initial temperatures, any initial chemical imbalance is quickly erased (Reisenegger 1995; Fernández & Reisenegger 2005), justifying our assumption of vanishing initial values for η_e and η_μ, which we also verified through simulations with non-zero initial values.

The effect of the initial rotation period P₀ is less obvious. First, we note that

$$\begin{array}{c}\hfill \frac{P_0^2}{P^2}=1-\frac{t}{\tau }·\end{array}$$(25)

Thus, as long as t ≲ τ, we have P ∼ P₀. During this time, the rotational properties of the pulsar do not vary much, so the probability of detecting it as a radio pulsar is also roughly constant, and the probability to detect it at an age < t is ∼t/τ. For the characteristic ages τ ∼ 10^7 − 10 yr, it is very unlikely for any of them to be young enough to explain its temperature through passive cooling, which requires t ≲ 10⁵ yr.

If vortex friction or crustal heating dominates, the luminosity at t ≫ 10⁵ yr depends just on the present spin-down parameters, $\dot{\mathrm{\Omega }}$ or $\mathrm{\Omega }\dot{\mathrm{\Omega }}$, respectively, not on the initial conditions. This imposes a very weak constraint on the initial period, which only needs to be slightly shorter than the present one.

For rotochemical heating with Cooper-paired nucleons, the initial period can have a strong qualitative effect by determining whether the chemical imbalances can grow enough to reach the threshold Δ_thr and thus activating the internal heating. For this to happen, it requires

$$\begin{array}{c}\hfill \frac{1}{P_0^2}-\frac{1}{P^2}\gtrsim \frac{1}{{(1.8\mathrm{ms})}^2}\frac{{\mathrm{\Delta }}_{\mathrm{thr}}}{1.5\mathrm{MeV}}·\end{array}$$(26)

For MSPs, it appears possible to satisfy this condition, as there are a few known MSPs with periods P ≈ 1.5 ms (e.g., Backer et al. 1982; Hessels et al. 2006; Bassa et al. 2017; Ho et al. 2019), although it is not obvious that all MSPs have had such fast rotations. Among CPs, the fastest known has a period P = 16 ms (Marshall et al. 1998), and therefore it is very unlikely for any of them to satisfy this condition.

For J0437, as discussed in Sect. 2, its age inferred from the temperature of its white dwarf companion, t = 6.0 ± 0.5 Gyr (Durant et al. 2012), is quite close to its (intrinsic) characteristic age, τ = 6.7 ± 0.2 Gyr, suggesting that its “initial” period (as it emerged from the accretion phase and started spinning down), P₀, was substantially shorter than its current period P. This makes it plausible, but by no means certain, that P₀ was short enough for the chemical imbalances to reach Δ_thr, activating rotochemical heating and explaining the observed temperature. For J2124, the absence of a binary companion does not allow for such an assessment.

6. Conclusions

Passive cooling models of NSs predict surface temperatures T_s^∞ ≪ 10⁵ K at ages t > 10⁶ yr. In the present paper, we have considered deep UV observations of five pulsars with characteristic ages much larger than this (all with τ > 10⁷ yr), four of which have detectable emission. In two cases (MSP J0437 and CP B0950), this can be confidently interpreted as thermal emission from their whole surface, leading to measurements of T_s^∞ ∼ 10⁵ K, which imply the presence of some heating mechanism. For the other three cases (MSP J2124 and CPs J0108 and J2144), we considered upper limits on T_s^∞ based on their UV detection or nondetection.

We modeled the thermal evolution of NSs with three potentially important heating mechanisms, namely rotochemical heating, vortex creep, and crustal heating, with and without Cooper pairing of neutrons and/or protons, but ignoring the possible presence of magnetized structures within the superconducting protons. We found that none of the heating mechanisms by itself can match all the observations, indicating that a realistic model needs to include other processes (such as external heating by relativistic magnetospheric particles or the presence of a magnetic field in the superconductor) or at least two mechanisms acting together. In fact, it is possible to account for all observational data with a model that combines vortex creep with rotochemical heating in the presence of a large, uniform⁷ and isotropic Cooper pairing gap ≈1.5 MeV in one of the nucleon species (or in both of them, with an equivalent combination Δ_thr ≈ 1.5 MeV) as well as a fast initial rotation for J0437 (P₀ ≲ 2 ms, from Eq. (26)). In this model, rotochemical heating is the dominant heating mechanism in MSPs, getting activated as they spin down enough for their chemical imbalances to exceed Δ_thr, thus explaining the temperature of J0437. CPs, on the other hand, are unlikely to have a fast enough initial rotation for this to happen, but the temperature of B0950 can be explained by vortex creep. In this model, the predicted temperatures for the three other pulsars are close to their observational upper limits, making further observations of these and other nearby, old pulsars in multiple filters worthwhile.

Data availability

The neutron star thermal evolution models listed in Table 3 are publicly available at Zenodo: https://zenodo.org/records/18378150.

Acknowledgments

This work was supported by ANID-FONDECYT grant 1201582 and the Center for Astrophysics and Associated Technologies (CATA; ANID Basal grant FB210003).

References

Appendix A: Thermal oscillations

As mentioned in Sect. 4.1, Petrovich & Reisenegger (2011) found that the evolution of an MSP with rotochemical heating in the presence of Durca reactions and a large, uniform, and isotropic Cooper pairing gap leads to strong oscillations of the chemical imbalances and the temperature. On the other hand, González-Jiménez et al. (2015) did not find oscillations when considering the more realistic case of nonuniform and anisotropic gaps, in which reactions become activated only in a small fraction of the core volume. We explore this effect for the case of Murca reactions with Δ_thr = 1.5 MeV in Fig. A.1, with reaction rates multiplied by different factors f ≤ 1. We see that, for large f, oscillations occur, whereas for smaller values they are eliminated and the MSP reaches a quasi-steady state that is nearly independent of f (compare the curves with f = 10⁻² and 10⁻⁵), with the transition happening somewhere in the range 10⁻² < f < 10⁻¹. We found the same behavior for Durca reactions with Δ_thr = 1.5 MeV, in this case with a threshold value f ∼ 10⁻⁸.

Fig. A.1. Evolution of the chemical imbalances and the temperature for an MSP with initial conditions P0 = 1 ms, T0∞ = 109 K, and ηe, 0∞ = ημ, 0∞ = 0, considering the effect of rotochemical heating with Murca reactions, normal protons, and superfluid neutrons with Δn∞ = 1.5 MeV. Each panel has the reaction rates reduced by a different factor f = 1, 10−1, 10−2, 10−5, respectively.

All Tables

All Figures

Fig. 1. Surface temperature evolution curves for each observed pulsar considering models with different reheating mechanisms. In each panel, all curves were calculated assuming magnetic dipole spin-down with a constant dipole moment inferred from the measured spin-down parameters P and Ṗ of the respective pulsars. B0950 and J0108 are plotted together, because they have nearly the same inferred dipole field. Their curves were generated with the magnetic field of B0950. For the MSPs (J0437 and J2124), the initial conditions were taken as P0 = 1 ms and T0∞ = 109 K, whereas for the CPs (B0950, J0108, and J2144) we assumed P0 = 5 ms and T0∞ = 1011 K. The (1σ) error bars or upper limits for the pulsars are placed at the times at which the present spin parameters are reached. The models used for Murca and Durca reactions are given in Table 2 as Model M and Model D, respectively. The labels correspond to the following cases (see Table 3 for a schematic summary): I Passive cooling with Murca reactions in normal (nonsuperfluid and nonsuperconducting) matter, solid black line. II Passive cooling with Durca reactions in normal matter, dashed gray line. III Rotochemical heating with Murca reactions in normal matter, solid maroon line. IV Rotochemical heating with Durca reactions in normal matter, dashed maroon line. V Vortex creep with Murca reactions setting J = 3 × 1043 erg s and assuming normal particles in the core, solid green line. VI Same as V, but with Durca reactions, dashed green line. VII Crustal heating (only present in MSPs) with Murca reactions and normal particles in the core, solid cyan line. VIII Same as VII, but with Durca reactions, dashed cyan line. IX Rotochemical heating with Murca reactions considering normal protons and superfluid neutrons with a uniform energy gap Δn = 1.5 MeV, with reaction rates reduced by a factor f = 10−2, solid orange line. X Same as IX, but with Durca reactions and reaction rates reduced by a factor f = 10−8, dashed orange line. XI Rotochemical heating with Murca reactions considering normal protons and superfluid neutrons with Δn = 1.5 MeV and reaction rates reduced by a factor f = 10−2, as well as vortex creep with J = 3 × 1043 erg s, magenta line with star markers.

In the text

	Fig. A.1. Evolution of the chemical imbalances and the temperature for an MSP with initial conditions P0 = 1 ms, T0∞ = 109 K, and ηe, 0∞ = ημ, 0∞ = 0, considering the effect of rotochemical heating with Murca reactions, normal protons, and superfluid neutrons with Δn∞ = 1.5 MeV. Each panel has the reaction rates reduced by a different factor f = 1, 10−1, 10−2, 10−5, respectively.
In the text