Cooling Theory Faced with Old Warm Neutron Stars: Role of Non-Equilibrium Processes with Proton and Neutron Gaps

TL;DR

This paper investigates how non-equilibrium beta processes, influenced by proton and neutron pairing gaps, can explain the unexpectedly high surface temperatures of old warm neutron stars through rotochemical heating.

Contribution

It introduces a detailed model including both proton and neutron pairing gaps to explain old warm neutron stars' temperatures via rotochemical heating.

Findings

Simultaneous inclusion of proton and neutron gaps enhances heating effects.

Model explains observed temperatures of old warm millisecond pulsars.

Initial rotational period influences temperature predictions for old neutron stars.

Abstract

Recent observations have found several candidates for old warm neutron stars whose surface temperatures are above the prediction of the standard neutron star cooling scenario, and thus require some heating mechanism. Motivated by these observations, we study the non-equilibrium beta process in the minimal cooling scenario of neutron stars, which inevitably occurs in pulsars. This out-of-equilibrium process yields the late time heating in the core of a neutron star, called the rotochemical heating, and significantly changes the time evolution of the neutron star surface temperature. To perform a realistic analysis of this heating effect, we include the singlet proton and triplet neutron pairing gaps simultaneously in the calculation of the rate and emissivity of this process, where the dependence of these pairing gaps on the nucleon density is also taken into account. We then compare the predicted surface temperature of neutron stars with the latest observational data. We show the simultaneous inclusion of both proton and neutron gaps is advantageous for the explanation of the old warm neutron stars since it enhances the heating effect. It is then found that the observed surface temperatures of the old warm millisecond pulsars, J2124-3358 and J0437-4715, are explained for various choices of nucleon gap models. The same setup is compatible with the observed temperatures of ordinary pulsars including old warm ones, J0108-1431 and B0950+08, by choosing the initial rotational period of each neutron star accordingly. In particular, the upper limit on the surface temperature of J2144-3933 can be satisfied if its initial period is $\gtrsim 10\,\mathrm{ms}$.

Figures

Tables

Table 1. Table 1: Observational data used in this work. t sd subscript 𝑡 sd t_{\text{sd}} , t kin subscript 𝑡 kin t_{\text{kin}} , T s ∞ superscript subscript 𝑇 𝑠 T_{s}^{\infty} , and P 𝑃 P denote the spin-down age, kinematic age, effective surface temperature, and period of neutron stars, respectively. The sixth column shows the atmosphere model used in the estimation of the surface temperature, where H, BB, C, and PL indicate hydrogen, blackbody, carbon, and power-law, respectively, while M represents a magnetized NS hydrogen atmosphere model, such as NSA [ 55 ] and NSMAX [ 56 ] . Data are taken from ATNF Pulsar Catalogue [ 57 , 58 ] unless other references are shown explicitly.

Name	${\log }_{10}{t}_{\mathrm{sd}}$	${\log }_{10}{t}_{\mathrm{kin}}$	${\log }_{10}{T}_s^{\infty }$	$P$	Atmos.
	[yr]	[yr]	[K]	[s]	model
PSR J2124-3358	${10.0}_{-0.1}^{+0.2}$ [27]		$4.7-5.3$ [27]	$4.9\times {10}^{-3}$	BB+PL
PSR J0437-4715	${9.83}_{-0.02}^{+0.01}$ [26]		$5.1-5.5$ [26]	$5.8\times {10}^{-3}$	BB
PSR J2144-3933	$8.5$		[30]	$8.5$	BB
PSR J0108-1431	$8.3$ [59]		$5.0-5.7$ [28]	$0.81$	BB
PSR B0950+08	$7.2$		$5.0-5.5$ [29]	$0.25$	BB+PL
RX J2143.0+0654	$6.6$		$6.082\left(2\right)$ [60]	$9.4$	BB
RX J0806.4-4123	$6.5$		$6.005\left(5\right)$ [61]	$11.4$	BB
PSR B1929+10	$6.5$		[62]	$0.23$	BB+PL
RX J0420.0-5022	$6.3$		$5.742\left(3\right)$ [63]	$3.5$	BB
PSR J2043+2740	$6.1$		$5.64\left(8\right)$ [64]	$0.096$	H
RX J1605.3+3249	$4.5$ [65]	$5.6-6.6$ [66]	$5.88\left(1\right)$ [67]	$3.4$ [65]	BB
RX J0720.4-3125	$6.3$	$5.8-6.0$ [68]	$5.5-5.7$ [69]	$8.4$	BB
RX J1308.6+2127	$6.2$	$5.5-6.2$ [70]	$6.07\left(1\right)$ [71]	$10.3$	BB
PSR B1055-52	$5.7$		$5.88\left(8\right)$ [64]	$0.20$	BB
PSR J0357+3205	$5.7$		${5.62}_{-0.08}^{+0.09}$ [72]	$0.44$	M+PL
RX J1856.5-3754	$6.6$	$5.66\left(5\right)$ [68]	$5.6-5.7$ [73]	$7.1$	BB
PSR J1741-2054	$5.6$		${5.85}_{-0.02}^{+0.03}$ [74]	$0.41$	BB+PL
PSR J0633+1748	$5.5$		$5.71\left(1\right)$ [75]	$0.24$	BB+PL
PSR J1740+1000	$5.1$		$6.04\left(1\right)$ [76]	$0.15$	BB
PSR B0656+14	$5.0$		$5.81\left(1\right)$ [77]	$0.38$	BB+PL
PSR B2334+61	$4.6$		$5.5-5.9$ [78]	$0.50$	M
PSR J0538+2817	$5.8$	$4.3-4.8$ [79]	$6.02\left(2\right)$ [79]	$0.14$	H
XMMU J1732-344		$4.0-4.6$ [80, 81]	$6.25\left(1\right)$ [81]		C
PSR B1706-44	$4.2$		${5.8}_{-0.13}^{+0.13}$ [82]	$0.10$	M+PL

Table 2. Table 2: The values of Z n ​ p subscript 𝑍 𝑛 𝑝 Z_{np} , Z n ​ p ​ ℓ subscript 𝑍 𝑛 𝑝 ℓ Z_{np\ell} , and W n ​ p ​ ℓ subscript 𝑊 𝑛 𝑝 ℓ W_{np\ell} in Eqs. ( 13 ) and ( 14 ), which are taken from Ref.  [ 24 ] .

$M$	$Z_{npe}$	$Z_{np\mu }$	$Z_{np}$	$W_{npe}$	$W_{np\mu }$
[$M_{\odot }$]	[${10}^{-61}\mathrm{erg}$]	[${10}^{-61}\mathrm{erg}$]	[${10}^{-61}\mathrm{erg}$]	[${10}^{-13}\mathrm{erg}{s}^2$]	[${10}^{-13}\mathrm{erg}{s}^2$]
$1.4$	$10$	$12$	$4$	$-1.5$	$-2$
$1.8$	$6$	$7$	$2$	$-1.4$	$-1.8$

Equations

$$C\frac{d{T}^{\infty }}{dt}$$

$$n+N_1\to p+N_2+\ell +{\bar{\nu }}_{\ell },$$

$$p+N_2+\ell \to n+N_1+{\nu }_{\ell },$$

$$Q_{M,N\ell }$$

$$\times \left[f_1{f}_2(1-f_3)(1-f_4)(1-f_{\ell })+(1-f_1)(1-f_2)f_3{f}_4{f}_{\ell }\right],$$

$$L_{\nu ,M}^{\infty }$$

$${ϵ}_N(\mathbf{p})\simeq {\mu }_N+\mathrm{sign}(p-p_{F,N})\sqrt{{\Delta }_N^2+v_{F,N}^2(p-p_{F,N}{)}^2},$$

$$\tilde{N}\tilde{N}\to [NN]+{\nu }_{\ell }+{\bar{\nu }}_{\ell },$$

$${\eta }_{\ell }\equiv {\mu }_n-{\mu }_p-{\mu }_{\ell },$$

$$Q_{M,N\ell }=Q_{M,N\ell }^{(0)}{I}_{M,ϵ}^N,$$

$$L_H^{\infty }$$

$$\Delta {\Gamma }_{M,N\ell }$$

$$\times \left[f_1{f}_2(1-f_3)(1-f_4)(1-f_{\ell })-(1-f_1)(1-f_2)f_3{f}_4{f}_{\ell }\right].$$

$$\Delta {\Gamma }_{M,N\ell }=\frac{Q_{M,N\ell }^{(0)}}{T(r)}{I}_{M,\Gamma }^N.$$

$$\frac{d{\eta }_e^{\infty }}{dt}$$

$$\frac{d{\eta }_{\mu }^{\infty }}{dt}$$

$$\dot{\Omega }(t)=-k\Omega (t{)}^n,$$

$$\Omega (t)=\frac{2\pi }{\sqrt{P_0^2+2P\dot{P}t}},$$

$$I_{0}\equiv 2\int_{0}^{\infty}dx_{\nu}\,x_{\nu}^{3}\,\biggl{[}\prod_{i}^{5}\int_{-\infty}^{\infty}dx_{j}\,f_{j}\biggr{]}\,\delta\biggl{(}\sum_{j=1}^{5}x_{j}-x_{\nu}\biggr{)}=\frac{11513\pi^{8}}{60480}~{}.$$

$$I_{M,ϵ}^N$$

$$\displaystyle\times\left[f(x_{\nu}-\xi_{\ell}-z_{n}-z_{p}-z_{N_{1}}-z_{N_{2}})+f(x_{\nu}+\xi_{\ell}-z_{n}-z_{p}-z_{N_{1}}-z_{N_{2}})\right]\delta^{3}\biggl{(}\sum_{j=1}^{5}\bm{p}_{j}\biggr{)}\,,$$

$$I_{M,\Gamma }^N$$

$$\displaystyle\times\left[f(x_{\nu}-\xi_{\ell}-z_{n}-z_{p}-z_{N_{1}}-z_{N_{2}})-f(x_{\nu}+\xi_{\ell}-z_{n}-z_{p}-z_{N_{1}}-z_{N_{2}})\right]\delta^{3}\biggl{(}\sum_{j=1}^{5}\bm{p}_{j}\biggr{)}\,,$$

$$x_N\equiv \frac{{ϵ}_N-{\mu }_N}{T},x_{\nu }\equiv \frac{{ϵ}_{\nu }}{T},v_N\equiv \frac{{\Delta }_N({\mathbf{p}}_{F,N},T)}{T},z_N\equiv \mathrm{sign}(x_N)\sqrt{x_N^2+v_N^2},{\xi }_{\ell }\equiv \frac{{\eta }_{\ell }}{T},$$

$$A_0^N$$

$$T_{c}(k_{F})=T_{0}\exp\biggl{[}-\frac{(k_{F}-k_{0})^{2}}{\Delta k^{2}}\biggr{]}~{},$$

$${\tilde{I}}_{M,ϵ}^N$$

$$\times \left[f(x_{\nu }-{\xi }_{\ell }-z_n-z_p-z_{N_1}-z_{N_2})+f(x_{\nu }+{\xi }_{\ell }-z_n-z_p-z_{N_1}-z_{N_2})\right],$$

$${\tilde{I}}_{M,\Gamma }^N$$

$$\times \left[f(x_{\nu }-{\xi }_{\ell }-z_n-z_p-z_{N_1}-z_{N_2})-f(x_{\nu }+{\xi }_{\ell }-z_n-z_p-z_{N_1}-z_{N_2})\right].$$

$${\tilde{I}}_{M,ϵ}^N$$

$$\times {\int }_{r_n}^{1-r_{N_1}-r_{N_2}-r_p}d{u}_n{\int }_{r_{N_1}}^{1-u_n-r_{N_2}-r_p}d{u}_{N_1}{\int }_{r_{N_2}}^{1-u_n-u_{N_1}-r_p}d{u}_{N_2}$$

$$\times \frac{u_n}{\sqrt{u_n^2-r_n^2}}\frac{u_{N_1}}{\sqrt{u_{N_1}^2-r_{N_1}^2}}\frac{u_{N_2}}{\sqrt{u_{N_2}^2-r_{N_2}^2}}{K}_{ϵ}(u_n+u_{N_1}+u_{N_2},r_p),$$

$${\tilde{I}}_{M,\Gamma }^N$$

$$\times {\int }_{r_n}^{1-r_{N_1}-r_{N_2}-r_p}d{u}_n{\int }_{r_{N_1}}^{1-u_n-r_{N_2}-r_p}d{u}_{N_1}{\int }_{r_{N_2}}^{1-u_n-u_{N_1}-r_p}d{u}_{N_2}$$

$$\times \frac{u_n}{\sqrt{u_n^2-r_n^2}}\frac{u_{N_1}}{\sqrt{u_{N_1}^2-r_{N_1}^2}}\frac{u_{N_2}}{\sqrt{u_{N_2}^2-r_{N_2}^2}}{K}_{\Gamma }(u_n+u_{N_1}+u_{N_2},r_p),$$

$$K_{ϵ}(u,r)$$

$$+\frac{1}{30}\left(16r^4+83r^2(1-u{)}^2+6(1-u{)}^4\right)\sqrt{(1-u{)}^2-r^2},$$

Full text

1 Introduction

The cooling theory of a neutron star (NS) has been studied for over the years and tested against observations of the thermal emission from NSs [1, 2, 3]. In particular, the minimal cooling paradigm [4, 5, 6] is successful in explaining these observations on the same footing. It consists of the early time ($t\lesssim 10^{5}\,\mathrm{yr}$) neutrino emission from the core and late time ($t\gtrsim 10^{5}\,\mathrm{yr}$) photon emission from the surface. Fast neutrino emission processes such as the direct Urca process are not included—this assumption is valid unless NSs are very heavy [7]. The minimal cooling takes account of the effect of the nucleon superfluidity111Throughout this paper, neutron superfluidity and proton superconductivity are collectively referred to as nucleon superfluidity. as it plays important roles in NS cooling (see, e.g., Refs. [8, 9, 10]). In the NS core, the spin-singlet pairing of protons and spin-triplet pairing of neutrons are formed once the internal temperature gets lower than the corresponding critical temperatures, giving rise to energy gaps in the spectra of elementary excitations of nucleons. Such an energy gap suppresses the slow neutrino emission [11, 12, 13]. In the meantime, the onset of nucleon pairings triggers the rapid neutrino emission from the Cooper pair breaking and formation (PBF) process [14, 15, 16, 17, 18, 19]. The formation of nucleon pairings also changes the contribution of nucleons to the heat capacity of the NS [20]. In any case, the neutrino emission rate gets highly suppressed as the NS cools down, and eventually the surface photon emission dominates the neutrino emission, which then results in a rapid decrease in the NS surface temperature to $T_{s}^{\infty}\lesssim 10^{4}\,\mathrm{K}$ at $t\gtrsim 1~{}\,\mathrm{Myr}$. This is a generic consequence of the cooling theory, which can in principle be tested via the observation of the temperature of old NSs.

This minimal cooling theory is build upon the assumption that nucleons and leptons (electrons and muons) in the NS core are in chemical (or beta) equilibrium through the Urca reactions. This assumption, however, turns out to be invalid for (especially old) NSs. As the rotation rate of a pulsar decreases, the centrifugal force decreases [21]. This reduction makes the NS continuously contract, which perturbs the local number density of each particle species away from the equilibrium value. On the other hand, the timescale of the Urca reaction is typically much longer than that of the NS contraction (especially for old NSs), and thus the beta equilibrium cannot be maintained. This has a significant impact on the cooling of NSs, since the imbalance in the chemical potentials of nucleons and leptons, which quantifies the degree of the departure from beta equilibrium, is partially converted to the heat inside the NS [21, 22, 23]. This heating mechanism due to the non-equilibrium Urca process is sometimes called the rotochemical heating. In the presence of the rotochemical heating, the NS surface temperature reaches $T_{s}^{\infty}\sim 10^{5}\,\mathrm{K}$ at $t\gtrsim\,\mathrm{Gyr}$ [24], which is in contrast to the prediction of the minimal cooling with the beta equilibrium.

Intriguingly, several recent observations of old pulsars suggest that the late time heating operates in NSs. It is reported in Refs. [25, 26] that the (far) UV spectrum of the millisecond pulsar (MSP), J0437-4715, is consistent with the thermal emission with $T_{s}^{\infty}=(1.25-3.5)\times 10^{5}\,\mathrm{K}$. The age of J0437-4715 is estimated independently from the pulsar spin-down age and the age of the companion white dwarf (WD), with both being compatible with the age of $t\simeq(6-7)\times 10^{9}\,\mathrm{yr}$. The surface temperature of $\sim 10^{5}\,\mathrm{K}$ for such an old NS cannot be explained by the cooling theory without heating. Further evidences of old “warm” NSs are provided by the MSP J2124-3358 [27] and ordinary pulsars J0108-1431 [28] and B0950+08 [29]. On the other hand, the observation of the ordinary old pulsar J2144-3933 imposes an upper limit on its surface temperature, $T_{s}^{\infty}<4.2\times 10^{4}\,\mathrm{K}$ [30], giving an evidence for the presence of an old “cold” NS. It is quite important to study if we can explain these observations by means of the minimal cooling with the non-equilibrium beta processes, i.e., with the rotochemical heating.222 There are several other heating mechanisms to explain these warm NSs, including superfluid vortex creep [31, 32, 33, 34, 35, 36] and rotationally-induced deep crustal heating [37]. See also Ref. [38] for comparison of different internal heating mechanisms.

To that end, it is necessary to include the non-equilibrium effect into the minimal cooling paradigm in a consistent manner. As mentioned above, the nucleon superfluidity plays a crucial role in the minimal cooling; this is also true for the rotochemical heating. The nucleon pairing gaps provide the threshold of the Urca reactions [39] because when the chemical imbalance is below this threshold, the Urca reactions are so suppressed that the heating is ineffective. In addition, the size of the threshold determines how much imbalance is to be generated, with a larger imbalance enhancing the heating effect. The effect of nucleon superfluidity on the non-equilibrium beta processes has been considered in the literature [40, 41, 42, 43]. In Refs. [40, 42], the rates and emissivities of the non-equilibrium beta processes were computed in the presence of either proton superconductivity or neutron superfluidity. In Ref. [41], the rotochemical heating was discussed for isotropic and uniform nucleon pairings. The anisotropic and density-dependent neutron pairing gaps were considered in Ref. [43], but the effect of proton superconductivity was not included in the computation. In reality, however, both the proton superconductivity and nucleon superfluidity are expected to occur in the NS core. The proton and neutron parings have the position- and temperature-dependence, and the neutron paring has also the angular-dependence. To perform a realistic simulation of the rotochemical heating, therefore, we need to include all of these effects simultaneously. Indeed, it is known in the calculation of the modified Urca processes that the simultaneous inclusion of all of these effects can considerably change the result [13]—this should also be the case for the non-equilibrium counterpart of the processes.

In this paper, we include both the proton singlet and neutron triplet pairings into the calculation of the NS thermal evolution in the presence of the non-equilibrium beta processes, with their density- and temperature-dependence fully taken into account. We then compare the results to the latest observational data of the NS surface temperature. We find that the heating with both nucleon pairings can be stronger than that only with neutron pairing [43], which turns out to be advantageous in explaining the old warm MSPs.

Meanwhile, it is found that the same setup can also account for the temperature data of the ordinary classical pulsars if the diversity of their initial spin periods is taken into account. We also discuss the compatibility of our scenario with the so-called X-ray dim Isolated Neutron Stars (XDINSs).

This paper is organized as follows. We review the minimal cooling and rotochemical heating in Sec. 2. In Sec. 3, we summarize the observational data of NS surface temperatures. We then give our numerical analysis in Sec. 4 and devote Sec. 5 to conclusions. In appendix A, we summarize the formulas for phase space factors. In Appendix B, we study the effect of varying magnetic field for the rotochemical heating.

2 Thermal evolution of neutron stars

In this section, we review the theory of NS thermal evolution [1, 2, 3]. Let us denote the local temperature in the core of a NS by $T(r)$, where $r$ denotes the distance from the NS center. This local temperature in general depends on the position, especially for a very young NS. It is however known [1, 2, 3] that the typical timescale of thermal relaxation in a NS is $\sim 10^{2-3}\,\mathrm{yr}$, and thus a NS with the age $t\gtrsim 10^{4}\,\mathrm{yr}$ can safely be regarded as isothermal. Since our main focus is on old NSs, in this work, we assume that NSs have already reached an isothermal state. In this case, the red-shifted internal temperature defined by $T^{\infty}\equiv T(r)e^{\Phi(r)}$ is constant throughout the NS core, with $e^{2\Phi(r)}=-g_{tt}(r)$ the time component of the metric. The evolution of this red-shifted temperature is then governed by

[TABLE]

where $L_{\nu}^{\infty}$ and $L_{\gamma}^{\infty}$ are the red-shifted luminosities of neutrino and photon emissions, respectively, $C$ denotes the total heat capacity of the NS, and $L_{H}^{\infty}$ is the source of heating if exists. In Sec. 2.1, we review the minimal cooling theory, assuming chemical equilibrium among nucleons and leptons. In this case, NSs just cool down due to the first two terms in the right-handed side of Eq. (1), with $L_{H}^{\infty}=0$. We then discuss in Sec. 2.2 the internal heating effect caused by the non-equilibrium beta processes in rotating pulsars, which yields a non-zero $L_{H}^{\infty}$. In Sec. 2.3, we summarize the pairing gap models for nucleon superfluidity used in this work, which affect the evaluation of $L_{\nu}^{\infty}$ and $L_{H}^{\infty}$ significantly.

2.1 Minimal cooling

The minimal cooling [4, 5, 6] is a successful paradigm that can explain many NS surface temperatures. In this scenario, the energy loss of a NS for $t\lesssim 10^{5}\,\mathrm{yr}$ is caused by the neutrino emission from the core, whose dominant processes are the modified Urca and PBF: $L_{\nu}^{\infty}\simeq L_{\nu,M}^{\infty}+L_{\nu,\text{PBF}}^{\infty}$. The fast neutrino emission processes such as the direct Urca process are not included. At later times, the surface photon emission becomes the dominant source for the NS cooling. In this subsection, we give a brief review on these processes.

The modified Urca process consists of the reactions

[TABLE]

where $N_{1}=N_{2}=n$ (neutron branch) or $N_{1}=N_{2}=p$ (proton branch) and $\ell=e,\mu$. The emissivity, the energy loss rate per unit time and volume, of this process is given by

[TABLE]

where $j=1,2,3,4$ denote the nucleons $n,N_{1},p,N_{2}$, respectively, $\delta^{4}(P_{f}-P_{i})$ the energy-momentum conserving delta function, $1/2\times\sum_{\mathrm{spin}}|\mathcal{M}_{M,N\ell}|^{2}$ the matrix element summed over all the particles’ spins with the symmetry factor, and $f$’s the Fermi-Dirac distribution functions. In the NS core, nucleons and charged leptons are highly degenerate, and thus only those in a thermal shell near the Fermi surface can participate in the reaction, making the emitted neutrino energy $\epsilon_{\nu}$ of $\mathcal{O}(T)$. This is much smaller than the Fermi momenta of nucleons and charged leptons, and thus we can neglect the neutrino momentum in the momentum delta function as well as in the matrix element. The matrix element is calculated in Ref. [44] based on the free one-pion exchange, with the angular dependence of the matrix element on the momenta of nucleons and leptons neglected. We also adopt this approximation in the following analysis.333 Recent development in this calculation suggests that there may be a significant change in the result if, e.g., the in-medium effects are included [45, 46, 47, 48]. We will study the impact of such effects on the calculation of the non-equilibrium beta processes on another occasion. We then obtain the total luminosity of the modified Urca process as

[TABLE]

where $dV=4\pi r^{2}e^{\Lambda(r)}dr$ with $e^{2\Lambda(r)}=g_{rr}(r)$ being the radial component of the metric.

The effect of the nucleon superfluidity enters through the distribution functions. In the presence of nucleon superfluidity, the energy of a quasi-particle near the Fermi surface has the form

[TABLE]

where $\mu_{N}$ is the chemical potential, $p_{F,N}$ the Fermi momentum, and $v_{F,N}$ the Fermi velocity. The energy gap $\Delta_{N}$ depends on the Fermi momentum $\bm{p}_{F,N}$ and temperature $T$: $\Delta_{N}=\Delta_{N}(\bm{p}_{F,N},T)$. The gap develops a non-zero value below the critical temperature, which results in suppression in the emissivity roughly by a factor of $e^{-\Delta_{N}/T}$. More detailed estimation of the suppression factor is given in Refs. [12, 13], where the beta equilibrium (equilibrium between the reactions (2) and (3)) is assumed. This assumption does not hold for spin-down pulsars, and the deviation from the beta equilibrium gives rise to heating in the NS core, as we see in Sec. 2.2.

As mentioned above, the minimal cooling does not include the fast cooling processes such as the direct Urca process, which can occur only in the very dense region of a heavy NS [7]. For instance, the Akmal-Pandharipande-Ravenhall (APR) equation of state (EOS) [49] allows the direct Urca to operate for $M\gtrsim 1.97\,M_{\odot}$. We use this EOS in the following analysis. If such a fast cooling process occurs, it dominates the other neutrino emission processes and causes a rapid cooling of the NS. At present, the masses of most of the isolated NSs whose surface temperature is measured are unknown. Since our main focus is on NSs with relatively high surface temperatures, in this work we just assume that the direct Urca process does not occur in the NSs we consider—this is equivalent to the assumption that the NSs in question have masses $\lesssim 1.97\,M_{\odot}$. We however note in passing that an observation of NS surface temperature that is lower than our prediction may be explained just by including the direct Urca process. Of course, the validity of this assumption should be confirmed by measuring the NS masses in the future observations.

In addition to the suppression of the modified Urca process, the nucleon superfluidity triggers another source of neutrino emission called the PBF process:

[TABLE]

where $\tilde{N}$ represents a quasi-particle of the nucleon Fermi liquid and $[NN]$ indicates the nucleon pairing. This PBF process can occur only in the presence of nucleon superfluidity, and is active when the temperature is just below the critical temperature [6]; for lower temperatures, this process is again suppressed by a Boltzmann factor. It turns out that the PBF process tends to be dominant in the neutrino emission after the onset of nucleon superfluidity [13]. In particular, other neutrino emission processes such as bremsstrahlung are usually sub-dominant. In our analysis, we include only the dominant neutrino emission processes, i.e., the modified Urca and the PBF, and neglect the other processes such as bremsstrahlung, for simplicity.

After the PBF process becomes ineffective (typically for $t\gtrsim 10^{5}\,\mathrm{yr}$), the NS cools via the surface photon emission. The photon emission obeys the Stefan-Boltzmann law: $L_{\gamma}=4\pi R^{2}\sigma_{B}T_{s}^{4}$, where $T_{s}$ is the surface temperature and $\sigma_{B}$ is the Stefan-Boltzmann constant. This surface temperature is different from the internal temperature as the envelope shields the surface from the interior. The relation between $T^{\infty}$ and $T_{s}$ depends on the constituent of the envelop, which is often parametrized by the ratio between the mass of the light elements in the envelop, $\Delta M$, and the NS mass, $M$: $\Delta M/M$.444The quantity $\eta\equiv g_{14}^{2}\Delta M/M$ is also used in the literature, where $g_{14}$ is the surface gravity in units of $10^{14}~{}\mathrm{cm}\cdot\mathrm{s}^{-2}$. We use the relation given in Ref. [50] in the following analysis.555For more recent discussions on the envelope models, see Refs. [51, 52].

For the heat capacity $C$ in Eq. (1), we consider the contribution of neutrons, protons, electrons, and muons. The contribution of nucleons is modified after the nucleon superfluidity occurs, and is significantly suppressed when the temperature is much lower than its critical temperature [1]. In this case, the contribution of leptons dominates the total heat capacity of the NS.

All in all, by solving Eq. (1) with $L_{\nu}^{\infty}$, $L_{\gamma}^{\infty}$, and $C$ given above, we obtain the temperature evolution of a NS. It is then found [4, 5, 6] that the resultant temperature curves are compatible with the observed surface temperatures of NSs, except for several old warm NSs mentioned in the introduction. We will discuss in the subsequent sections that these exceptions can also be explained in the minimal setup, once the effect of the non-equilibrium beta processes is taken into account.

2.2 Rotochemical Heating

The rotochemical heating occurs in the very same setup as in the minimal cooling. The only difference is whether the beta equilibrium of the Urca processes is assumed or not. In the case of an actual pulsar, its rotational rate keeps decreasing, which results in a continuous reduction of the centrifugal force. Consequently, NSs contract continuously and the equilibrium number densities of nucleons and charged leptons change at all times. The number densities of these particles in a NS follow the equilibrium values if the Urca reactions are fast enough. It however turns out that the typical timescale of these reactions is much longer than that of the NS contraction, especially for an old NS. As a result, the Urca reactions are no longer in beta equilibrium. The departure from the beta equilibrium leads to an imbalance in the chemical potentials of nucleons and leptons, which keeps increasing over the time. The energy stored in the chemical imbalance is released partly via the neutrino emission and partly as a heat, generating a non-zero $L_{H}^{\infty}$ in Eq. (1). We emphasize that this heating effect due to the non-equilibrium beta reactions is an inevitable consequence for actual rotating NSs and thus needs to be included into the calculation of the temperature evolution of NSs. We review this heating mechanism in Sec. 2.2.1, and discuss the effect of superfluidity on the rotochemical heating in Sec. 2.2.2.

2.2.1 Heating rate

As in the minimal cooling, we consider a NS which is comprised of nucleons, electrons, and muons. We again assume that the fast neutrino emission processes do not operate in the NSs in question. As we discussed above, the Urca processes in spin-down pulsars are out of equilibrium, which leads to imbalance in the chemical potentials of nucleons and leptons. We denote this imbalance by

[TABLE]

for $\ell=e,\mu$. This parameter quantifies the departure from the beta equilibrium and is equal to the amount of the energy released in each reaction.

The neutrino emissivity in the out-of-equilibrium modified Urca process is again calculated by using Eq. (2.1), with relaxing the assumption of the beta equilibrium. Following Refs. [41, 43], we factorize the emissivity as

[TABLE]

where $Q_{M,N\ell}^{(0)}$ denotes the emissivity computed for a non-superfluid NS in beta equilibrium. The dimensionless function $I^{N}_{M,\epsilon}$ represents the phase space integral, which incorporates the effect of the departure from beta equilibrium as well as that of nucleon superfluidity.

The departure from beta equilibrium generates heat due to the entropy production [21, 22, 23]. The heating luminosity is written as

[TABLE]

where $\Delta\Gamma_{M,N\ell}$ is the difference between the reaction rates of the processes (2) and (3):

[TABLE]

We again factorize $\Delta\Gamma_{M,N\ell}$ by using the phase space factor $I^{N}_{M,\Gamma}$ as follows [41, 43]:

[TABLE]

We study the phase space factors $I^{N}_{M,\epsilon}$ and $I^{N}_{M,\Gamma}$ in more detail in Sec. 2.2.2.

The time evolution of $T^{\infty}$ is still described by Eq. (1), now with a heating source $L_{H}^{\infty}$. The chemical imbalance parameters $\eta_{\ell}$ are also time-dependent. It is discussed in Ref. [39] that the diffusion timescale of the chemical imbalance is shorter than the time scale of its evolution, and thus we can safely assume that the red-shifted imbalance parameters $\eta_{\ell}^{\infty}\equiv\eta_{\ell}e^{\Phi(r)}$ are constant throughout the NS core. The evolution of $\eta_{\ell}^{\infty}$ is then determined by [24]

[TABLE]

where $\Omega(t)$ is the angular velocity of the pulsar. $Z_{np}$, $Z_{np\ell}$ and $W_{np\ell}$ are the constants defined in Ref. [24], which depend on the stellar structure of the NS. The first terms in the right-hand side of these equations describe the equilibration by the modified Urca process, while the second terms show the effect of spin-down, which makes the NS out of beta equilibrium.

For the pulsar spin-down, we assume the power-law deceleration:

[TABLE]

where $k$ and $n$ are positive constants. Furthermore, we take $n=3$ in the following analysis; this braking index corresponds to the case where the loss of the rotational energy is caused solely by the magnetic dipole radiation. In this case, we can solve Eq. (15) as

[TABLE]

where $P$ and $\dot{P}$ are the present values of the pulsar period and its derivative, respectively, and $P_{0}$ is the initial period $P(t=0)$. Note that the combination $P\dot{P}=-4\pi^{2}\dot{\Omega}/\Omega^{3}=4\pi^{2}k$ is time-independent for $n=3$ and related to the dipole magnetic filed; for instance, $B\sim 3.2\times 10^{19}\,(P\dot{P}/\mathrm{s})^{1/2}\,\mathrm{G}$ for a NS of radius $R=10\,\mathrm{km}$ and moment of inertia $I=10^{45}\,\mathrm{g}\,\mathrm{cm^{2}}$. We also find from Eq. (16) that if $P_{0}\ll P$, then $t=P/(2\dot{P})$. This is called the spin-down age of a pulsar, which is widely used for the estimation of the pulsar age.

2.2.2 Phase space factors

Now let us give more concrete expressions for the phase space factors $I^{N}_{M,\epsilon}$ and $I^{N}_{M,\Gamma}$ defined in Eqs. (9) and (12), respectively:666The prefactor is equal to the inverse of

I_{0}\equiv 2\int_{0}^{\infty}dx_{\nu}\,x_{\nu}^{3}\,\biggl{[}\prod_{i}^{5}\int_{-\infty}^{\infty}dx_{j}\,f_{j}\biggr{]}\,\delta\biggl{(}\sum_{j=1}^{5}x_{j}-x_{\nu}\biggr{)}=\frac{11513\pi^{8}}{60480}~{}.

(17)

We also note that $I^{N}_{M,\Gamma}$ in Eq. (19) has the opposite sign to the corresponding phase space integrals given in Refs. [41] and [43], while it is consistent with those in Ref. [40].

[TABLE]

where we have defined

[TABLE]

and $f(x)=1/(e^{x}+1)$. The constant factor $A_{0}^{N}$ is given by

[TABLE]

where $\Theta_{p}=1$ for $p_{F,n}<3p_{F,p}+p_{F,\ell}$ and [math] for otherwise, and $j=1,2,\dots,5$ correspond to $n,p,\ell,N_{1},N_{2}$, respectively. As can be seen from the definition in Eqs. (9) and (12), we have $I^{N}_{M,\epsilon}=1$ and $I^{N}_{M,\Gamma}=0$ in the limit of $v_{n}=v_{p}=\xi_{\ell}=0$, i.e., for a non-superfluid NS in beta equilibrium. The angular integral in the above equations is non-trivial in the presence of a neutron triplet pairing since its gap amplitude depends on the direction of $\bm{p}_{F,n}$ (see Sec. 2.3); otherwise it is just reduced to the factor $A_{0}^{N}$.

For superfluid matter, the phase space factors provide a threshold for the modified Urca process: $\eta_{\ell}\gtrsim 3\Delta_{n}+\Delta_{p}$ for the neutron branch and $\eta_{\ell}\gtrsim\Delta_{n}+3\Delta_{p}$ for the proton branch [39]. Hence, $\Delta_{\mathrm{th}}=\mathrm{min}\{3\Delta_{n}+\Delta_{p},\Delta_{n}+3\Delta_{p}\}$ is the threshold of the rotochemical heating—for $\eta_{\ell}\lesssim\Delta_{\mathrm{th}}$, heating does not occur because the modified Urca reaction is suppressed ($\Delta\Gamma_{M,N\ell}\simeq 0$). For a very young NS, the Urca reaction is fast enough so that the chemical equilibrium is maintained, i.e., $\eta_{\ell}=0$. Later, the NS departs from beta equilibrium due to the spin-down and $\eta_{\ell}$ monotonically increases until it exceeds $\Delta_{\mathrm{th}}$, after which the accumulated $\eta_{\ell}$ is converted to heat.777We however note that if $\Delta_{\mathrm{th}}$ is large and/or the increase rate of $\eta_{\ell}$ is small, $\eta_{\ell}$ may never exceed the threshold and thus rotochemical heating is always ineffective. We also see such situations in Sec. 4.3. Therefore, the rotochemical heating is efficient usually at late times, when $T\ll\Delta_{N}$ and $T\ll\eta_{\ell}$. In such a situation, we can safely exploit the zero temperature approximation in the calculation of the phase space factors [41].

A numerical calculation in Ref. [24] shows that the late-time heating indeed occurs in a non-superfluid NS. In this case, we have analytical expressions for the phase space factors [21]. Such analytical expressions for the neutron (proton) branch can also be obtained for the case where only protons (neutrons) form a constant paring gap in the limit of $T\to 0$ [41]. As for the numerical evaluation of the phase space factors, Refs. [40, 42] give the results for the case in which either proton or neutron has a non-zero gap. Reference [41] also performs the numerical computation of the phase space factors in the presence of the neutron and proton singlet uniform pairings. In Ref. [43], neutron triplet pairings whose gap has density and temperature dependence are considered, but the effect of proton superfluidity is neglected. In this work, we include the effect of both the singlet proton and triplet neutron pairing gaps with taking account of their density and temperature dependence. For the calculation of the phase space factors, we use the zero temperature approximation as in Ref. [41]; see Appendix A for more details.

2.3 The nucleon gap models

The strength of neutrino emission and the rotochemical heating depends on the size of nucleon gaps through the phase space factors. In the NS core, singlet pairing of protons and triplet pairing of neutrons are formed once the internal temperature has fallen below their critical temperatures. For proton singlet pairing, the gap is independent of the orientation of the proton momentum $\bm{p}_{F,p}$. On the other hand, neutron triplet gap is dependent on the angle $\theta$ between the neutron momentum $\bm{p}_{F,n}$ and the quantization axis. This angular dependence is different among states which have different projections $m_{J}$ of the total angular momentum on the quantization axis. Currently, it is unclear which of the states is actually formed in the NS core. In this work, we just consider the $m_{J}=0$ state, which has the angular dependence of $\Delta_{N}(\bm{p}_{F,N},T)\propto\sqrt{1+3\cos^{2}\theta}$. The application of our analysis to the case of the $|m_{J}|=2$ state is straightforward.888See Ref. [43] for a comparison of the rotochemical heating in the cases of $m_{J}=0$ and $|m_{J}|=2$.

Calculation of gap profile suffers from uncertainty due to the nucleon interactions and medium effects. There are many theoretical calculations/models of the zero temperature gap, $\Delta_{N}(\bm{p}_{F,N},T=0)$, which is proportional to the critical temperature. Given the zero temperature gap, we can obtain the gap amplitude for $T\neq 0$ by using the fitting formula in Ref. [1]. For the proton gap, several theoretical models are shown in Ref. [8], among which we use the AO model [53] as a small gap and the CCDK model [54] as a large gap. The calculation of neutron triplet gap is more uncertain, so we just use the phenomenological models named “a”, “b” and “c” in Ref. [4], for which the critical temperature is given by a Gaussian-shape function with respect to the Fermi wave number $k_{F}\equiv p_{F,N}/\hbar$:

[TABLE]

with $T_{0}=1.0\times 10^{9}\,\mathrm{K}$, $k_{0}=1.8\,\mathrm{fm}^{-1}$, and $\Delta k=0.5\,\mathrm{fm}^{-1}$ for “a”, $T_{0}=3.0\times 10^{9}\,\mathrm{K}$, $k_{0}=2.0\,\mathrm{fm}^{-1}$, and $\Delta k=0.5\,\mathrm{fm}^{-1}$ for “b”, and $T_{0}=1.0\times 10^{10}\,\mathrm{K}$, $k_{0}=2.5\,\mathrm{fm}^{-1}$, and $\Delta k=0.7\,\mathrm{fm}^{-1}$ for “c”. In Fig. 1, we show the critical temperature/zero temperature gap of the models used in our analysis as a function of the Fermi wave number $k_{F}$. As for the neutron triplet pairing gap, we show the amplitude for $\cos\theta=0$. The dashed line indicates the core-crust boundary and the dotted lines correspond to the NS center for $M=1.4\,M_{\odot}$ and $1.8\,M_{\odot}$. As can be seen from the plots, even at low temperatures, protons can be in the normal fluid state near the NS center depending on the NS mass and gap models, while neutrons are in the superfluid state in the entire core for any gap model.

3 Observations of neutron star temperatures

Before going to our numerical study, we summarize the current status of the NS temperature observations. We focus on isolated NSs with the age $t>10^{4}\,\mathrm{yr}$, for which we can safely assume that the thermal and diffusion relaxation in the NS core has already been completed. In Tab. 1, we list the NSs whose surface temperature is measured, together with two NSs (PSR J2144-3933 and PSR B1929+10) for which only the upper bound on the surface temperature is obtained. For most of the NSs in the table, only the spin-down age $t_{\mathrm{sd}}=P/(2\dot{P})$ can be used for the estimation of their age, while in some cases, the kinematic age $t_{\mathrm{kin}}$, which is derived from the motion of the supernova remnant, is also available. In this work, we use the kinematic age if available, and the spin-down age otherwise.999We note that once we fix the initial period and the value of $P\dot{P}$ (hence the dipole magnetic field), the NS age is unambiguiously determined from the observed value of the NS period through Eq. (16). In particular, the spin-down age should agree to the real pulsar age without uncertainty if $P_{0}\ll P$. The values of $t_{\mathrm{sd}}$ and $P$ are taken from ATNF Pulsar Catalogue [57, 58] unless other references are shown explicitly.

The surface temperature is estimated from the fit of the observed spectrum with atmosphere models. Generically speaking, the hydrogen/carbon atmosphere model fits the thermal emission of young pulsars ($t\lesssim 10^{5}\,\mathrm{yr}$) well, while the heavy element model or blackbody offers a good fit for a middle-aged or old pulsar. In addition, if a NS has a relatively large magnetic field, a magnetized NS hydrogen atmosphere model such as NSA [55] and NSMAX [56] may improve its spectrum fit. In the sixth column in Tab. 1, we show the atmosphere model used in the evaluation of the surface temperature of each NS; H, BB, C, PL, and M represent the hydrogen, blackbody, carbon, power-law, and magnetized NS hydrogen atmosphere models, respectively. For some NSs in the table, there are several atmosphere models that can fit the observed spectrum. In such cases, we choose the model which is considered to give the best fit and/or yields a NS radius of a plausible size ($\sim 10$ km). Moreover, there are some cases where the use of two or more blackbody components improves the fit due to the presence of hot spots, though we do not show this explicitly in the table. In these cases, we use the temperature of the component with a NS radius that is consistent with the typical NS radius. Note, however, that the fits with atmosphere models are often performed with the mass, radius, or distance of the NS being fixed, though these quantities are not precisely known for most of the NSs—if this is the case, additional systematic uncertainty might be present. For more details, see the references cited in the fourth column of the table.

The first two pulsars in Tab. 1, PSR J2124-3358 and J0437-4715, are classified into the MSPs. They have small $P$ and $\dot{P}$, and hence a small dipole magnetic field, compared to the ordinary (classical) pulsars. PSR J0437-4715 is the closest millisecond pulsar at present. Its rotational period, mass, and distance are estimated to be $5.76$ ms, $1.44\pm 0.07~{}M_{\odot}$, and $d=156.79\pm 0.25$ pc, respectively [83]. This pulsar is in a binary system accompanied with a white dwarf. The spin-down age of PSR J0437-4715 is estimated in Ref. [26] to be $t_{\text{sd}}=(6.7\pm 0.2)\times 10^{9}$ years with the Shklovskii correction [84] included. This is in a good agreement with the estimated age of the white dwarf, $t_{\text{WD}}=(6.0\pm 0.5)\times 10^{9}$ years. In Ref. [26], it is found that a fit with the Rayleigh-Jeans law in the far UV range is consistent with a blackbody emission from the whole NS surface with a temperature of $T_{s}^{\infty}=(1.25-3.5)\times 10^{5}$ K. As argued in Ref. [26], it is unlikely for the observed surface temperature to be due to heat flow coming from the magnetosphere regions. Since the minimal cooling theory predicts $T_{s}^{\infty}\ll 10^{3}\,\mathrm{K}$ for $t\sim 10^{10}\,\mathrm{yr}$, the observed surface temperature requires late time heating. PSR J2124-3358 is an isolated millisecond pulsar with a period of $4.93$ ms [83]. Its spin-down age, after corrected by the Shklovskii effect, is $11^{+6}_{-3}\times 10^{9}$ years [27] for the distance $d=410^{+90}_{-70}$ pc. Its surface temperature is obtained with a blackbody plus power-law fit to be $(0.5-2.1)\times 10^{5}$ K [27], with the radius fixed to be 12 km. This is also well above the cooling theory prediction.

We also have examples of old warm ordinary pulsars: PSR J0108-1431 and B0950+08. PSR J0108-1431 is an old NS with the spin-down age of $2.0\times 10^{8}$ years [59], where the Shklovskii correction is taken into account. The analysis in Ref. [28] with a Rayleigh-Jeans spectrum fit shows that its surface temperature is $T_{s}^{\infty}=(7-10)\times 10^{4}\,(d_{130}/R_{13})^{2}$ K, where $d_{130}$ is the distance in units of 130 pc and $R_{13}$ is the apparent radius in units of 13 km. Within the error of the distance, $d=210^{+90}_{-50}$ pc [85], the maximum (minimum) temperature is estimated as $T_{s}^{\infty}=5.3\times 10^{5}$ K ($1.1\times 10^{5}$ K) for a radius of 13 km. PSR B0950+08 has the spin-down age of $1.75\times 10^{7}$ years. Its surface temperature is obtained with a power-low plus blackbody spectrum fit in Ref. [29] as $(1-3)\times 10^{5}$ K, with other parameters such as the pulsar radius varied in a plausible range.

Seven middle aged pulsars, RX J2143.0+0654, J0806.4-4123, J0420.0-5022, J1308.6+2127, J0720.4-3125, J1856.5-3754, and J1605.3+3249 are classified into the X-ray Dim Isolated Neutron Stars (XDINSs), which are also dubbed as the Magnificent Seven. They exhibit thermal X-ray emission without any signature of magnetospheric activity, and have a rather long spin period. See Refs. [86, 87, 88] for reviews of XDINSs. Their surface temperatures are found to be $\sim 10^{6}\,\mathrm{K}$, which are again higher than the prediction of the cooling theory. The inferred dipole magnetic fields of these NSs are relatively large: $B\sim 10^{13}-10^{14}\,\mathrm{G}$. As we will discuss in the following sections, in these NSs, a different type of heating mechanism due to the magnetic field decay [89, 90] may operate.

Contrary to the above examples, PSR J2144-3933 is an old “cool” NS. This is one of the slowest pulsars, having $P=8.51$ s, and its spin-down age is $t_{\mathrm{sd}}=333\,\mathrm{Myr}$ with the Shklovskii correction taken into account. Assuming the blackbody spectrum, the authors in Ref. [30] obtained an upper limit on the surface temperature of J2144-3933: $T_{s}^{\infty}<4.2\times 10^{4}\,\mathrm{K}$. This is the lowest limit on the surface temperature of NSs for the moment.

In the next section, we discuss if the minimal cooling setup with the non-equilibrium beta processes is compatible with the observed surface temperatures in Tab. 1.

4 Results

Now we show the results of our numerical analysis, where we follow the thermal evolution of NSs with the effect of the non-equilibrium beta reactions included. We then compare the results with the observed surface temperatures given in Tab. 1.

4.1 Physical input

We perform the numerical study in the framework of the minimal cooling with the non-equilibrium beta process discussed in Sec. 2. The following inputs are common to all of the analyses:

• •

APR EOS [49].

• •

Initial condition: $T^{\infty}=10^{10}\,\mathrm{K}$ and $\eta_{e}^{\infty}=\eta_{\mu}^{\infty}=0$.

• •

Protons and neutrons form singlet (${}^{1}S_{0}$) and triplet (${}^{3}P_{2}\,(m_{J}=0)$) pairings in the core, respectively.

• •

The pulsar braking index $n=3$, i.e., $\Omega(t)$ obeys Eq. (16).

• •

We use the formula given in Ref. [50] for the relation between $T^{\infty}$ and $T_{s}$.

For the superfluid gap models, we use the CCDK and AO models for proton and the “a”, “b” and “c” models for neutron, which are shown in Fig. 1. The values of $Z_{np}$, $Z_{np\ell}$, and $W_{np\ell}$ in Eqs. (13) and (14) are read from Fig. 3 in Ref. [24], which are summarized in Tab. 2. The numerical values of EOS and the solution of Tolman- Oppenheimer-Volkoff (TOV) equation are taken from NSCool [91].

We divide the NSs listed in Tab. 1 into two categories: MSPs and the others. The latter contains ordinary pulsars and XDINSs. We exploit a representative parameter set for each category as follows:

Millisecond pulsars

MSPs have much smaller $P$ and $\dot{P}$ than ordinary pulsars. With MSP J0437-4715 in mind, we use the following parameters for this category:

• •

$M=1.4\,M_{\odot}$.

• •

$P=5.8\,\mathrm{ms}$.

• •

$\dot{P}=5.7\times 10^{-20}$.

• •

$\Delta M/M=10^{-7}$.

We also note that the values of $P$ and $\dot{P}$ of J2124-3358, $P=4.9\,\mathrm{ms}$ and $\dot{P}=2.1\times 10^{-20}$, are fairly close to those of J0437-4715, while its mass is unknown. We have fixed the amount of the light elements in the envelope, $\Delta M/M=10^{-7}$, as it turns out that the result is almost independent of this choice for old NSs such as J0437-4715 and J2124-3358.

Ordinary pulsars and XDINSs

For ordinary pulsars and XDINSs, we use

• •

$M=1.4\,M_{\odot}$ or $1.8\,M_{\odot}$.

• •

$P=1\,\mathrm{s}$.

• •

$\dot{P}=1\times 10^{-15}$.

• •

$\Delta M/M=10^{-7}$ or $10^{-15}$.

Note that $P$ and $\dot{P}$ affect the rotochemical heating only through Eq. (16), and thus the result depends only on the combination $P\dot{P}$. Ordinary pulsars have $P\dot{P}\sim 10^{-17}-10^{-13}$, corresponding to $B\sim 10^{11}-10^{13}\,\mathrm{G}$. The dependence of the thermal evolution on $P\dot{P}$ is weaker than that on gap models and $P_{0}$, and thus we fix it to be $P\dot{P}=1\times 10^{-15}\,\mathrm{s}$ in the following analysis.

Once we fix the NS parameters as above, the time evolution of the NS surface temperature depends only on the nucleon gap models and the initial period $P_{0}$. As we see in Sec. 2.2.2, the heating rate depends on the nucleon pairing gaps via the phase space factors. On the other hand, the initial period $P_{0}$ affects the time evolution of the NS angular velocity $\Omega(t)$ in Eq. (16), through which the accumulation rate of the chemical imbalance is modified. We will study these effect in the following subsections.

4.2 Millisecond pulsars

We first compute the evolution of the redshifted surface temperature $T_{s}^{\infty}$ for the MSP category. The resultant temperature evolution is shown in Fig. 2, where the initial period is taken to be $P_{0}=1\,\mathrm{ms}$. We use the CCDK (AO) model for proton pairing in the left (right) panel. The red, blue, and green lines correspond to the “a”, “b”, and “c” models for neutron pairing, respectively. The solid (dashed) lines are for the case with (without) proton superfluidity. The observed surface temperatures of J0437-4715 and J2124-3358 are also shown by the black points with black solid lines indicating the uncertainty.

In all the cases given in Fig. 2, the chemical imbalance $\eta_{\ell}^{\infty}$ overcomes the threshold of rotochemical heating $\Delta_{\mathrm{th}}$ at $t=10^{6-7}\,\mathrm{yr}$, at which the surface temperatures quickly rise to $T_{s}^{\infty}\sim 10^{6}\,\mathrm{K}$. After that the system reaches a quasi-steady state, where the increase in $\eta_{\ell}^{\infty}$ due to spin-down is compensated by the consumption via the Urca processes and the heating rate balances with the photon cooling rate. During this stage, $T_{s}^{\infty}$ gradually decreases due to the decline of the term $|\Omega\dot{\Omega}|$.

As we discussed in Sec. 2, in this work we include the effect of both proton and neutron superfluidity simultaneously. We can see the consequence of this simultaneous inclusion by comparing the solid and dashed lines for each case. Let us first study the case with the proton CCDK for $P_{0}=1\,\mathrm{ms}$, i.e., the left panel in Fig. 2. In this case, for all neutron gap models, the heating effect starts to be visible at a later time in the presence of proton superconductivity. This is because the additional contribution from the proton pairing gap increases the threshold of rotochemical heating $\Delta_{\mathrm{th}}$, and thus an extra amount of $\eta_{\ell}^{\infty}$ needs to be accumulated. This results in a delay in the onset of rotochemical heating. Moreover, a larger value of $\Delta_{\mathrm{th}}$ leads to a larger value of $\eta_{\ell}^{\infty}$ eventually. This then results in a higher $T_{s}^{\infty}$ at late times, since the heating power is proportional to $\eta_{\ell}$ as in Eq. (10). This feature can also be seen for every neutron gap model in the left panel in Fig. 2.

Next, we examine the cases with the AO proton gap model shown in the right panel in Fig. 2. For this proton paring gap, we do not see enhancement in $T_{s}^{\infty}$ due to the proton superfluidity for the neutron gaps “a” and “b”. As we see in the left panel in Fig. 1, the size of the AO proton gap is smaller than the CCDK gap, and even vanishes deep inside the NS core. For this reason, the rotochemical threshold $\Delta_{\mathrm{th}}$ is determined almost solely by the neutron gap, which makes the effect of proton superconductivity invisible. For the neutron “c” gap, on the other hand, the gap amplitude is very large near the NS center, and hence heating is ineffective there. Instead, the rotochemical heating mainly occurs in the intermediate region where the AO proton gap is sizable, which makes the effect of proton gap manifest. This observation indicates that it is crucial to take account of the density dependence of the nucleon pairing gaps for the evaluation of the rotochemical heating effect.

The results in Fig. 2 show that the observed surface temperatures of J0437-4715 and J2124-3358 can be explained by the heating effect of non-equilibrium beta reactions, especially for moderate/large nucleon gaps. In particular, for the neutron “c” gap model, the simultaneous inclusion of proton superconductivity improves the fit considerably such that the predicted thermal evolution is totally consistent with the observed temperatures.

4.3 Ordinary pulsars and XDINSs

Next we discuss the second category, which is comprised of ordinary pulsars and XDINSs. In Fig. 3, we show the time evolution of the surface temperature $T_{s}^{\infty}$ for this category with $P_{0}=1\,\mathrm{ms}$. We use the neutron gap “a” (“c”) in the upper (lower) panels and the CCDK (AO) proton gap in the left (right) panels. The red and blue lines correspond to $M=1.4\,M_{\odot}$ and $1.8\,M_{\odot}$, while the solid and dashed lines represent the cases for the heavy ($\Delta M/M=10^{-15}$) and light ($\Delta M/M=10^{-7}$) element envelope models, respectively. It is found that another choice of the envelope parameter $\Delta M/M$ just falls in between these two cases. We also plot the observed surface temperatures of ordinary pulsars and XDINSs with black and green points, respectively, with the horizontal (vertical) lines indicating the uncertainty in the kinematic age ($T_{s}^{\infty}$). The bars with down arrows correspond to the upper limits on $T_{s}^{\infty}$ of J2144-3933 and B1929+10.

For the cases shown in Fig. 3, the rotochemical heating begins earlier than $10^{4}\,\mathrm{yr}$. The difference from the MSP category is due to the larger $P\dot{P}$ of the ordinary pulsars and XDINSs. The predicted temperatures tend to be higher for the light element envelope than the heavy element one at earlier times, but this difference disappears at later times. In addition, the evolution curves depend on the NS mass. This dependence is, however, rather non-trivial compared with that observed in Ref. [43], where only the neutron triplet paring is taken into account and the temperature always gets higher for a lighter NS mass. This complexity is caused by the non-trivial dependence of the rotochemical threshold on the matter density, as both proton and neutron pairing gaps contribute to the threshold.

It is found from the top two panels that the small neutron gap “a” offers the predictions consistent with most of the surface temperatures of the ordinary pulsars. We also find that three hot XDINSs, having relatively high temperatures $T_{s}^{\infty}\sim 10^{6}\,\mathrm{K}$ at $t\sim 10^{6}\,\mathrm{yr}$, are located above the predictions. The bottom two panels, on the other hand, show that the large neutron gap “c” gives higher temperatures than the “a” model, and can be consistent only with a part of the NS temperatures shown in the figure. The hot XDINSs are still above the prediction, while several young pulsars ($t\sim 10^{4}\,\mathrm{yr}$ and $T_{s}^{\infty}<10^{6}\,\mathrm{K}$) are below the thermal evolution curves.

A different choice of $P_{0}$ may change the evolution of $T_{s}^{\infty}$ via the modification of the increase rate in $\eta_{\ell}^{\infty}$. In general, a larger initial period makes the departure from the beta equilibrium smaller and hence the heating effect milder. To see this, in Fig. 4, we show the time evolution of the surface temperature for $P_{0}=10\,\mathrm{ms}$ for the CCDK (left) and AO (right) proton gap models and the “a” (top) or “c” (bottom) neutron gap models. It is found that this slightly larger initial period strongly suppresses the heating effect. This can be seen in the bottom panels, where the heating does not occur for both of the NS masses due to the large neutron gap of “c”. For the neutron gap “a”, on the other hand, the internal heating can occur when the nucleon pairing gaps are sufficiently suppressed in the NS center so that even a small amount of the chemical imbalance can trigger the rotochemical heating. The predictions in all of these cases are consistent with the surface temperatures of relatively young pulsars ($t\sim 10^{4}-10^{6}\,\mathrm{yr}$) of $10^{5}\,\mathrm{K}<T_{s}^{\infty}<10^{6}\,\mathrm{K}$. The surface temperature of B0950+08 may be explained for some cases. The upper limit on the surface temperature of J2144-3933 can also be satisfied.

We note that, for most of the NS parameter choices, at later times, the intensity of the rotochemical heating has little dependence on $P\dot{P}$ and thus on its dipole magnetic field, as discussed in detail in App. B. This is because for a sufficiently old pulsar, the spin-down effect is very weak and the evolution of $\eta_{\ell}^{\infty}$ is solely determined by the equilibration caused by the modified Urca process. This observation, therefore, indicates that we cannot explain the low surface temperature of J2144-3933 by changing the magnetic field while keeping $P_{0}$ small.

In summary, the observed surface temperatures of ordinary pulsars are compatible with the rotochemical heating for both small and large neutron gaps. For a small neutron gap, the predicted curves of the time evolution of $T_{s}^{\infty}$ tend to be in good agreement with most of the observations for $P_{0}=1\,\mathrm{ms}$. In the case of a large neutron gap, on the other hand, the middle-aged and old pulsars can still be explained with $P_{0}=1\,\mathrm{ms}$, while it is required to assume a larger initial period for young pulsars; for an initial period of $1\,\mathrm{ms}<P_{0}<10\,\mathrm{ms}$, the heating can still occur but predicted temperatures tend to be lower than those for $P_{0}=1\,\mathrm{ms}$, with which we can explain all of the data below the curves shown in Fig. 3. In particular, the J2144-3933 limit can be satisfied for any nucleon gap models if we take a sufficiently large initial period. It is intriguing to note that the J2144-3933 has one of the longest periods observed so far, as mentioned in Sec. 3. This observation may suggest that the initial period of this NS is also relatively long so that the rotochemical heating has never operated and the NS has been cooled down below the limit $T_{s}^{\infty}<4.2\times 10^{4}\,\mathrm{K}$.

On the other hand, the observed temperatures of some of the XDINSs may be higher than the prediction of the rotochemical heating. This discrepancy could be merely due to the errors of their age and temperature; in particular, the temperatures shown in the Tab. 1 may suffer from large systematic uncertainty because of our ignorance of their masses, radii, and distances. We also note that the temperature data of the three especially hot XDINS—J2143+0654, RX J0806.4-4123, and RX J1308.6+2127—used in this analysis are obtained using a spectrum fit with one blackbody component plus Gaussian absorption lines. If there is a hot spot in these NSs, however, such a fit tends to give a higher temperature and a smaller NS radius than the actual ones. In fact, the inferred radii of J2143+0654 and RX J0806.4-4123, $R=3.10(4)\,\mathrm{km}$ for $d=500\,\mathrm{pc}$ [60] and $R=2.39(15)\,\mathrm{km}$ for $d=250\,\mathrm{pc}$ [61], respectively, are considerably smaller than the typical size of the NS radius $\sim 10-15\,\mathrm{km}$. In addition, these NSs exhibit X-ray pulsations [65], which indicate that the temperature distribution on their surface is inhomogeneous. It is indeed found that a spectrum analysis based on a two-temperature blackbody model gives a larger radius and a lower surface temperature for these XDINSs [92], with which the discrepancy between the prediction and observation is significantly reduced. Another potential explanation for the discrepancy is that the XDINSs have undergone the magnetic field decay, which makes the spin-down age differ from the actual age. Moreover, the magnetic field decay itself can be another source of heating [93, 90]. If this is the case, the inclusion of only the rotochemical heating effect may be insufficient to explain the temperatures of the XDINSs.

5 Conclusions

We have studied the non-equilibrium beta process in the minimal cooling scenario, which gives rise to the late time heating in NSs. Extending the previous works, we have included the singlet proton and triplet neutron pairing gaps simultaneously in the calculation of the rate and emissivity of the process, with their density dependence taken into account. We then compare the time evolution of the NS surface temperature predicted in this framework with the latest observations of the NS temperatures, especially with the recent data of the old warm NSs. It is found that the simultaneous inclusion of both proton and neutron gaps is advantageous for the explanation of the old warm NSs, since it increases the threshold of rotochemical heating and thus enhances the heating effect. We find that the observed surface temperatures of warm MSPs, J2124-3358 and J0437-4715, are explained for various choices of nucleon gap models. The same setup can also explain the temperatures of ordinary pulsars by choosing the initial rotational period of each NS accordingly. In particular, with $P_{0}=10\,\mathrm{ms}$ or larger, the upper limit on the surface temperature of J2144-3933 can be satisfied.

To explain the observation, we require $P_{0}\lesssim 10\,\mathrm{ms}$ for old warm pulsars and $P_{0}\gtrsim 10\,\mathrm{ms}$ for old cold one. This assumption is reasonable for hot MSPs since their current periods are also as small as $\mathop{\mathcal{O}}(1)\,\mathrm{ms}$. On the other hand, for ordinary pulsars, several recent studies suggest that they are born with $P_{0}=\mathop{\mathcal{O}}(10-100)\,\mathrm{ms}$ [94, 95, 96, 97, 98, 99, 100, 101], which is apparently in tension with the requirement of $P_{0}\lesssim 10\,\mathrm{ms}$. Nevertheless, the initial period of a NS highly depends on the detail of the supernova process where the NS was created. Given that a fully satisfactory simulation for supernova explosion process has not yet been available, we regard the issue of the NS initial period as an open question, expecting future simulations and astrophysical observations to answer this problem.

As we have seen in Sec. 4.3, the surface temperatures of the three XDINSs J2143+0654, RX J0806.4-4123, and RX J1308.6+2127 are higher than the prediction of the rotochemical heating. This discrepancy may be due to large magnetic fields of these XDINSs. In fact, XDINSs commonly have absorption feature in the X-ray spectrum [102, 103], which is interpreted as the proton cyclotron resonance or atomic transition. Both require the magnetic field larger than about $10^{13}\,\mathrm{G}$, and such a strong magnetic field can affect the thermal evolution of XDINSs. If this is the case, we need a more involved analysis for these XDINSs including the magnetic field evolution.

Regarding the ambiguity in the choice of pairing gaps, there is a hint from the observation. The recently observed rapid cooling of the NS in the supernova remnant Cassiopeia A (Cas A) may be the first observational signature of the neutron superfluidity in the NS core [104, 105, 106, 107, 52]. The measured decline in temperature is explained with a large proton gap such as the CCDK model and a small neutron gap with the critical temperature of $T_{c}\simeq 5\times 10^{8}\,\mathrm{K}$, which is about a factor of 2 smaller than that for the “a” model. The time evolution of $T_{s}^{\infty}$ for these gaps is close to that for the proton CCDK + neutron “a” model, and thus is consistent with the observations of ordinary pulsars as shown in Fig. 3. This setup can also explain the surface temperature of the MSP J2124-3358, while the predicted value of $T_{s}^{\infty}$ for J0437-4715 is slightly below the observed one. Although a large neutron gap such as “c” model is favored by the observation of these MSPs, it is not consistent with the Cas A NS observation. Further exploration of a unified explanation for the temperatures of the Cas A NS and the old warm NSs will be given on another occasion.

Finally, we emphasize that the non-equilibrium beta process discussed in this paper is not an adhoc assumption but an inevitable consequence of rotating NSs, and therefore the heating mechanism based on this process should always be taken into account. In this sense, the minimal cooling plus rotochemical heating is the minimal scenario for the NS thermal evolution. Intriguingly, we have found that this “minimal” setup is compatible with the observed surface temperatures of NSs for the moment, without relying on exotic physics. Further developments in the evaluation of nucleon pairing gaps, as well as additional data of NS surface temperatures, allow us to test this minimal scenario in the future.

Acknowledgements

We thank Teruaki Enoto, Kenji Fukushima, and Hideyuki Umeda for valuable discussions and suggestions. We also thank Mikhail E. Gusakov for useful correspondence. This work is supported in part by the Grant-in-Aid for Scientific Research A (No.16H02189 [KH]), Young Scientists B (No.17K14270 [NN]), Innovative Areas (No.26104001 [KH], No.26104009 [KH], No.18H05542 [NN]). The work of KY was supported by JSPS KAKENHI Grant Number JP18J10202.

Appendix

Appendix A Phase space factors

In this appendix, we give a detailed discussion on the phase space integrals in Eqs. (18) and (19). While the nucleon pairing is negligible, we use the formulas in Ref. [21]. In the numerical calculation, we neglect the pairing if $3v_{n}+v_{p}<1$ for the neutron branch and if $v_{n}+3v_{p}<1$ for the proton branch. If the pairing operates but $\xi_{\ell}<1$, we neglect non-equilibrium effects. We use the results in Ref. [13] for $I^{N}_{M,\epsilon}$, and set $I^{N}_{M,\Gamma}=0$. If the pairing operates and $\xi_{\ell}\geq 1$, we numerically perform the integrals using the zero temperature approximation; we show the detail of this calculation in this appendix. We have checked that the result scarcely depends on the choice of the threshold value of pairing and $\xi_{\ell}$.

A.1 Energy integral

We define the energy integral part of the phase space factors by101010$\tilde{I}^{N}_{M,\Gamma}$ in Eq. (24) has the opposite sign to $I^{N}_{M,\Gamma}$ in Ref. [41].

[TABLE]

Following the argument in Ref. [41], we use the zero temperature approximation in this calculation; namely, we replace the Fermi-Dirac distribution $f(x)$ by the step function $\Theta(-x)$. The phase space factors are then written in terms of the following integral expressions:

[TABLE]

where $r_{N}=v_{N}/\xi$, $u_{N}=|z_{N}|/\xi$, and

[TABLE]

Notice that the integral region of Eqs. (25) and (26) is different from that in Eqs. (B.12) and (B.13) in Ref. [41]. We also find that the sign of the first term in Eq. (27) is opposite to that in Eq. (B.11) in Ref. [41].

A.2 Angular integral

Next, we perform the angular integration. Due to the different angular dependence of proton and neutron pairing gaps, we need to separately treat the integrals for the proton and neutron branches.

Proton branch

For the proton branch, we can carry out the angular integration of the momenta of the protons and leptons trivially since the proton singlet gap is isotropic. The integration with respect to the neutron momentum direction is then reduced to a simple average as

[TABLE]

where only $v_{n}$ is dependent on $\Omega_{n}$.

Neutron branch

The neutron branch involves three neutrons, $n$, $n_{1}$ and $n_{2}$, so the angular integral is generically complicated. Here we neglect the proton and lepton momenta in the momentum conserving delta function. Then, the three neutron momenta form an equilateral triangle, and we can put the angular integral into

[TABLE]

with $\cos\theta_{n}$ the polar angle of $n$ around the quantization axis and $\varphi_{n_{1}}$ the azimuthal angle of $n_{1}$ around $\bm{p}_{n}$. The relative angles among the momenta of $n$, $n_{1}$ and $n_{2}$ are fixed because of the momentum conservation. The polar angles of $\bm{p}_{n_{1}}$ and $\bm{p}_{n_{2}}$ with respect to the quantization axis are written as

[TABLE]

respectively. The integrands in Eq. (30) depend on $\cos\theta_{n}$, $\cos\theta_{n_{1}}$, and $\cos\theta_{n_{2}}$ through the triplet neutron pairing gap $\Delta_{n}\propto\sqrt{1+3\cos^{2}\theta}$.

Appendix B Rotochemical heating with different magnetic field for ordinary pulsars

In Sec. 4, we fix the spin-down rate $P\dot{P}=10^{-15}\,\mathrm{s}$, i.e., the dipole magnetic field $B\sim 10^{12}\,\mathrm{G}$, for ordinary pulsars, though it generically takes a value in a rather broad range: $P\dot{P}=10^{-17}-10^{-13}\,\mathrm{s}$, corresponding to $B\sim 10^{11-13}\,\mathrm{G}$. Since this value controls the evolution of angular velocity (see Eq. (16)), a change in $P\dot{P}$ may affect the intensity of the rotochemical heating. In this appendix, we study the effects of varying $P\dot{P}$ on the temperature evolution of ordinary pulsars. We fix the other NS parameters to be $M=1.4\,M_{\odot}$ and $\Delta M/M=10^{-15}$ in the following analysis.

Figure 5 shows the thermal evolution for $P\dot{P}=5.8\times 10^{-17}\,\mathrm{s}$ (blue), $4.2\times 10^{-15}\,\mathrm{s}$ (red), and $1\times 10^{-13}\,\mathrm{s}$ (purple), corresponding to $B\sim 10^{11}\,\mathrm{G}$, $10^{12}\,\mathrm{G}$, and $10^{13}\,\mathrm{G}$, respectively. Note that $5.8\times 10^{-17}\,\mathrm{s}$ is the observed value of $P\dot{P}$ for B0950+08 and $4.2\times 10^{-15}\,\mathrm{s}$ is that for J2144-3933. We use the neutron “a” and proton CCDK (left panels) or AO (right panels) gap models. The top, middle, and bottom panels show the cases of $P_{0}=1\,\mathrm{ms}$, $5\,\mathrm{ms}$, and $10\,\mathrm{ms}$, respectively. For $P_{0}=1\,\mathrm{ms}$ and $5\,\mathrm{ms}$, we see that the surface temperature is lower for a larger $P\dot{P}$ at middle age. This is because with a larger $P\dot{P}$, rotochemical heating begins earlier, resulting in a smaller $\eta_{\ell}$ at middle age. For $t\gtrsim 10^{8}\,\mathrm{yr}$, however, the surface temperatures of different $P\dot{P}$ converge to each other since the spin-down effect becomes so small that the evolution of $\eta_{\ell}$ is mainly determined by the modified Urca process. For $P_{0}=10\,\mathrm{ms}$, the rotochemical heating never occurs for the proton CCDK gap model (bottom left). In the case of the proton AO model, on the other hand, the heating effect is small but still visible (bottom right). Moreover, the difference of surface temperatures at $t\gtrsim 10^{6}\,\mathrm{yr}$ is larger than that for $P_{0}=1\,\mathrm{ms}$ or $P_{0}=5\,\mathrm{ms}$. In this case, the accumulated $\eta_{\ell}$ is close to the rotochemical threshold, and hence the heating luminosity is very sensitive to the difference in the $\eta_{\ell}$, making the temperature evolution highly dependent on other NS parameters such as $P\dot{P}$.

Notice that although a large $P\dot{P}$ (and thus a large dipole magnetic field) tends to predict low surface temperature, it does not help explain the low temperature of J2144-3933 compared to that of J0108-1431 and B0950+08, as indicated by the red and blue curves in Fig. 5, which are always close to each other. Consequently, the different temperatures of these NSs should be attributed to the difference in $P_{0}$.