Vortex Creep Heating in Neutron Stars

Motoko Fujiwaraa,b*** E-mail address<EMAIL_ADDRESS>Koichi
Hamaguchia,c††† E-mail address<EMAIL_ADDRESS>Natsumi
Nagataa‡‡‡ E-mail address<EMAIL_ADDRESS>
and Maura E. Ramirez-Quezadaa,d§§§ E-mail address: me.quezada@hep-
th.phys.s.u-tokyo.ac.jp,

aDepartment of Physics, University of Tokyo, Bunkyo-ku, Tokyo 113–0033, Japan

b Physik-Department, Technische Universität, München, James-Franck-Straße,
85748 Garching, Germany
cKavli IPMU (WPI), University of Tokyo, Kashiwa, Chiba 277–8583, Japan

dDual CP Institute of High Energy Physics, C.P. 28045, Colima, México

Recent observations of old warm neutron stars suggest the presence of a
heating source in these stars, requiring a paradigm beyond the standard
neutron-star cooling theory. In this work, we study the scenario where this
heating is caused by the friction associated with the creep motion of neutron
superfluid vortex lines in the crust. As it turns out, the heating luminosity
in this scenario is proportional to the time derivative of the angular
velocity of the pulsar rotation, and the proportional constant $J$ has an
approximately universal value for all neutron stars. This $J$ parameter can be
determined from the temperature observation of old neutron stars because the
heating luminosity is balanced with the photon emission at late times. We
study the latest data of neutron star temperature observation and find that
these data indeed give similar values of $J$, in favor of the assumption that
the frictional motion of vortex lines heats these neutron stars. These values
turn out to be consistent with the theoretical calculations of the vortex-
nuclear interaction.

##  1 Introduction

A neutron star (NS) is incredibly dense and exists under extreme conditions of
pressure and temperature that cannot be found in other places in the universe.
While the internal structure of NSs remains elusive, indirect evidence
suggests the existence of a neutron superfluid in the inner crust. The first
indication of superfluidity in NSs came from the observation of non-zero
pairing energy associated with attractive forces, leading to the formation of
an energy gap and hence superfluidity [1]. This scenario was predicted to
occur in stars with neutron cores [2]. After the discovery of pulsars, the
energy gap in NS matter has been further studied; see, e.g., Refs. [3, 4, 5]
for a recent review on superfluidity in NSs.

In a rotating NS, the irrotational property of a superfluid requires the
formation of vortex lines, whose distribution determines the angular velocity
of the superfluid component. In the inner crust region, these vortex lines are
fixed at certain positions by the interactions with nuclei and cannot move
freely. This preserves the rotational speed of the superfluid component and
prevents it from following the slowdown of the pulsar rotation, giving rise to
the deviation in the rotational speed between the superfluid and other
components. This deviation increases until the vortex lines are forced to move
by the Magnus force, which increases as the difference in the rotational speed
increases. This vortex-line dynamics leads to some observational consequences.
A well-known example is the pulsar glitch phenomenon, namely, sudden changes
in the rotational frequency of NSs,111See Refs. [6, 7, 8, 9] for a recent
review on pulsar glitches.  which could be attributed to an avalanche of
unpinning of superfluid vortex lines [10, 11]. Another phenomenon is the
heating effect caused by the friction associated with the creep motion of
vortex lines [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22], which is the
subject of the present paper.

Our motivation to revisit this heating mechanism is provided by the recent
observations of old warm NSs [23, 24, 25, 26, 27, 28, 29], whose observed
temperature is considerably higher than that predicted in the standard NS
cooling scenario [30, 31, 32, 33, 34, 35]. This work aims to study whether
these observations can be explained by the vortex-creep heating effect. With
this objective in mind, we focus on the following characteristic property of
the vortex-creep heating. As we see below in detail, the heating luminosity in
this heating mechanism is proportional to the time derivative of the angular
velocity of the pulsar rotation, and the proportional constant is determined
only by the NS structure and the vortex-nuclear interactions. As a result, the
value of this proportional constant, denoted by $J$ in this paper, is almost
universal over NSs. In addition, we can obtain the value of $J$ for an old NS
by observing its temperature and its pulsar motion since, at late times, the
heating luminosity balances with the luminosity of photon emission, which is
determined by the surface temperature of the NS. As it turns out, the present
data of old warm NSs indeed show similar values of $J$ for these stars, in
agreement with the prediction of the vortex-creep heating scenario. We also
find that these values are compatible with the $J$ parameter evaluated from
the calculations of the nuclear pinning force available in the literature.

The remainder of this paper is structured as follows. In Section 2, we provide
a brief overview of the thermal evolution of NSs, focusing on the isothermal
phase. In Section 3, we describe the vortex creep heating mechanism for NSs
and explain how we can relate a universal parameter with the late-time
temperature prediction. In Section 4, we summarize the numerical evaluation of
the pinning force, from which we calculate the parameter $J$. In Section 5, we
study the recent observations of old and warm NSs to assess the current status
of the vortex creep heating hypothesis. Finally, we summarize our findings and
conclude our discussion in Section 6.

##  2 Thermal evolution of neutron star

This section reviews the surface temperature prediction of NSs based on their
thermal evolution and sources of heating and cooling. The temperature
distribution within NSs is characterized by a local core temperature $T$,
which exhibits a temperature gradient only during the early stages of the NS,
typically within the first 10–100 years of its existence [36, 37]. The high
thermal conductivity of the highly degenerate electron gas causes the core to
become isothermal over time, reaching thermal equilibrium. Hence, the red-
shifted temperature of NSs,
$T^{\infty}(\bar{r},t)=T(\bar{r},t)e^{\phi(\bar{r})}$, where $\phi(\bar{r})$
specifies the gravitational redshift, reaches a constant value
$T^{\infty}(\bar{r},t)\simeq T^{\infty}(t)$ and only the outermost layers
exhibit an appreciable temperature gradient. The relativistic equation
describing the thermal evolution of NSs after thermal relaxation is given by
[38, 39]

$\displaystyle
C(T^{\infty})\frac{dT^{\infty}}{dt}=-L_{\nu}^{\infty}(T^{\infty})-L_{\gamma}^{\infty}(T^{\infty})+L_{\rm
H}^{\infty}.$ (2.1)

Here, $C$ represents the total heat capacity of the NS and is temperature-
dependent. The right-hand side of the equation expresses the red-shifted
luminosity for three different processes, namely, neutrino cooling
$L_{\nu}^{\infty}$, photon cooling $L_{\gamma}^{\infty}$, and heating source
$L_{\rm H}^{\infty}$.

At temperatures below a few $\times~{}10^{9}~{}\mathrm{K}$, the NS becomes
transparent to neutrinos, allowing them to escape without interacting with the
stellar matter and carrying away energy. Therefore, the cooling process during
the early stage of the star’s life is dominated by neutrino emission. At later
times, typically for $t\gtrsim 10^{5}~{}\mathrm{yrs}$,222The dominance of
photon emission might be delayed for massive NSs, in which the rapid neutrino
emission via the direct Urca process could occur; we do not consider this
possibility in what follows.  photon emission dominates neutrino emission. The
thermal photon emission follows a blackbody spectrum. The photon luminosity
$L_{\gamma}$ can be described by the Stefan-Boltzmann law and related to
surface temperature $T_{\rm s}$ as

$\displaystyle L_{\gamma}=4\pi R_{\rm NS}^{2}\sigma_{\rm SB}T_{\rm s}^{4},$
(2.2)

in the local reference frame of the NS, where $\sigma_{\rm SB}$ is the
Stephan-Boltzmann constant and $R_{\rm NS}$ the NS radius. To relate the
internal temperature $T$ of a NS to its surface temperature $T_{\rm s}$, we
use the heat envelope model proposed by Potekhin et al. in 1997 [40].
According to this model, the observed thermal emission from old isolated NSs
can be explained by the heat trapped in a thin envelope surrounding the star’s
crust.

If we have an internal heating source, the heating luminosity will balance
with the photon luminosity at a sufficiently late time,

$\displaystyle L^{\infty}_{\rm H}$ $\displaystyle\simeq L^{\infty}_{\gamma},$
(2.3)

which determines the surface temperature. It is pointed out that NSs have some
internal heating mechanisms, such as vortex creep heating, rotochemical
heating, and magnetic field decay. These heating effects may become visible at
late times and can operate even for isolated NSs. These internal mechanisms
are comprehensively compared with the observed surface temperature in Refs.
[20, 22], which is recently revisited by Ref. [41].

We note in passing that the balance equation (2.3) generically holds with good
accuracy for old NSs that are older than $10^{5}~{}\mathrm{yrs}$. This can be
seen by estimating the typical timescale $\tau_{\rm eq}$ for the relaxation
into the equilibrium state:

$\displaystyle\tau_{\rm eq}$ $\displaystyle\simeq\frac{CT^{\infty}}{4\pi
R_{\rm NS}^{2}\sigma_{\rm SB}T_{\rm s}^{4}}$ (2.4) $\displaystyle\sim 3\times
10^{4}~{}\mathrm{yrs}\left(\frac{C}{10^{35}~{}\mathrm{erg/K}}\right)\left(\frac{R_{\rm
NS}}{11.43~{}{\rm km}}\right)^{-2}\left(\frac{T_{\rm
s}}{10^{5}~{}\mathrm{K}}\right)^{-4}\left(\frac{T^{\infty}}{10^{6}~{}\mathrm{K}}\right)~{},$
(2.5)

where we have used $T^{\infty}\sim 10^{6}~{}\mathrm{K}$ corresponding to the
surface temperature, at the equilibrium phase, of $T_{\rm s}\sim
10^{5}~{}\mathrm{K}$.333To derive the typical scale of $\tau_{\rm eq}$, we
fitted the relation between the NS internal temperature $T$ and the NS surface
temperature $T_{\rm s}$ for $t_{\rm age}\lesssim 10^{6}~{}\mathrm{yrs}$ [40]
by assuming $T\propto T_{\rm s}^{2}$. It is found that this timescale is
shorter than the age of old NSs, assuring the equilibrium condition (2.3).

##  3 Review of vortex creep heating

In this section, we will review vortex creep heating, where the presence of a
superfluid in the inner crust of a NS plays a key role. In this region, vortex
lines are thought to exist as a consequence of the NS rotation. In Sec. 3.1,
we will derive the equation of motion for this rotational motion by
introducing two different angular velocities for the inner crust superfluid
and the other part of the star. In Sec. 3.2, we will consider the dynamics of
a vortex line and evaluate its radial velocity. Finally, in Sec. 3.3, we will
assess the effect of vortex creep on the late-time temperature prediction of
NSs.

###  3.1 Equations of motion for crust and superfluid

To describe a rotating NS, let us divide it into two components depending on
how external torque exerts on it, the crust component and the superfluid
component [42]. The crust component comprises a rigid crust, a lattice of
nuclei, and charged particles tightly coupled to electromagnetic field lines
[42]. This component is directly affected by external torque provided by
pulsar magnetic radiation. The superfluid component refers to ${}^{1}S_{0}$
superfluid of neutrons.444The neutron triplet (${}^{3}P_{2}$) superfluid in
the core region is classified into the crust component in the two-component
model [42]. This is because the neutron superfluid in this region is expected
to coexist with the proton superconductor and be tightly coupled to the crust
component [43, 44]. This superfluid phase is believed to appear in the inner
crust region based on pairing gap evaluations [3, 4, 5]. This component is
just indirectly affected by the external torque through the interaction with
the crust component.

Figure 1: The structure of a NS and vortex lines. Left: The grey, blue, and
brown regions represent the outer crust, inner crust, and core regions,
respectively. The red lines represent vortex lines. Right: A single vortex
line in the inner crust. The vortex line in the neutron superfluid is attached
to the outer crust and has two boundaries.

On the left side of Fig. 1, we show a schematic diagram of a NS. The grey
region is the outer crust, composed of ions and electrons. The blue layer
represents the inner crust region and has an approximately $\sim
1~{}\mathrm{km}$ thickness. In this region, nuclei, electrons, and neutrons
exist, and neutrons are expected to be in the form of the neutron singlet
superfluid. The brown region is the NS core, whose internal structure remains
uncertain and a subject of ongoing research and debate. In the two-component
model, the neutron superfluid in the blue layer is classified into the
superfluid component, while all the other parts are classified into the crust
component.

Based on this two-component model, we derive the equations of motion and
describe the rotation. For later convenience, we divide the system into a thin
disc and introduce the cylindrical coordinate $(r,\varphi,z)$. Since we treat
the crust component as a rigid body, its angular velocity $\Omega_{\rm c}(t)$
is independent of $r$. On the other hand, the superfluid angular velocity
varies with $r$, and we denote it by $\Omega_{\rm s}(t,r)$. The equation of
motion for the crust component is

$\displaystyle I_{\rm c}\dot{\Omega}_{\rm c}(t)=N_{\rm ext}(t)+N_{\rm
int}(t),$ (3.6)

where $I_{\rm c}$ represents the moment of inertia of the crust component, and
the dot represents the derivative with respect to time $t$. On the right-hand
side, we divide the torque into two parts: The first term $N_{\rm ext}(t)$
represents the external torque acting on the system. The second term $N_{\rm
int}(t)$ corresponds to the effects of internal torques in the system

$\displaystyle N_{\rm int}(t)$ $\displaystyle=-\int
dI_{p}(r)\frac{\partial\Omega_{\rm s}}{\partial t}(t,r),$ (3.7)

where $dI_{p}$ represents the differential inertial momenta. The integral is
taken over the region where the crust component interacts with the superfluid
component, called the pinning region. These two components connect through the
pinning force between the vortex line and the nuclei-like object in the inner
crust, as we will discuss in Sec. 3.2.

The equations of the superfluid component are obtained by introducing two
fundamental measures of rotation in the fluid: vorticity and circulation. The
vorticity vector $\bm{\omega}$ characterizes the microscopic features of
rotation and is a locally determined value defined as the curl of the fluid
velocity $\bm{v}$,

$\displaystyle\bm{\omega}\equiv\nabla\times\bm{v}.$ (3.8)

In contrast, the circulation $\Gamma$ measures the macroscopic rotation and is
defined over a finite region. It is given by the line integral of the fluid
velocity $\bm{v}$ around a closed path $C$, which can be expressed as a
surface integral over any surface $S$ with the boundary $C$,

$\displaystyle\Gamma\equiv\oint_{C}\bm{v}\cdot
d\bm{\ell}=\iint_{S}\bm{\omega}\cdot d\bm{S},$ (3.9)

where $d\bm{\ell}$ and $d\bm{S}$ denote the line and surface elements,
respectively. We used Stokes’ theorem to obtain the last expression.

Superfluid motion obeys the potential flow condition,

$\displaystyle\nabla\times\bm{v}{{}_{\rm s}}=\bm{0},$ (3.10)

where $\bm{v}{{}_{\rm s}}$ denotes the superfluid velocity. This condition,
implying the absence of the vorticity in the superfluid, holds because the
superfluid velocity is proportional to the gradient of the phase of the
condensate wave-function of the superfluid. Nevertheless, we still have a
nonzero circulation if there exists a singular object known as a vortex line.
In Fig. 1, we show a schematic picture of a single vortex line in the NS inner
crust. The vortex line is a string-like configuration with a thickness of the
order of femtometers (red curve). The circulation for each vortex line is
quantized in units of

$\displaystyle\kappa$ $\displaystyle\equiv\frac{h}{2m_{n}},$ (3.11)

where $h$ is the Planck constant and $m_{n}$ is the neutron mass. This
quantization follows from the condition that the wave-function of the
condensate is single-valued, and thus the change in its phase must be $2\pi k$
with $k$ as an integer. Since a vortex line with $k=1$ is energetically
favored and stable, the system with larger angular velocity contains a larger
number of vortex lines [45]. The number of vortex lines will be saturated if
the total circulation reaches that of the rigid rotation as a whole system,
which is expected in NSs. The vortex lines in the inner crust have boundaries
corresponding to the normal matter in the outer crust. Under this
circumstance, the circulation for the contour $C$ in Fig. 1 is uniquely
determined and, because of the potential flow condition (3.10), can be
regarded as topological as it remains unchanged under deformations of $C$
(unless it passes through another vortex line).

We may express the superfluid velocity on average in the same form as the
normal fluid,555This relation is confirmed by observing the shape of the free
surface for rotating superfluid in liquid He system [46, 47, 48].

$\displaystyle\Braket{\bm{v}_{\rm s}}=\bm{\Omega}_{\rm s}\times\bm{r}~{},$
(3.12)

where $\bm{r}$ denotes the position vector from the center of the NS. The
total circulation of a superfluid system is equal to the sum of the
circulation of each vortex line. This means that the number of vortex lines is
directly related to the superfluid angular velocity $\bm{\Omega}_{\rm s}$. By
substituting Eq. (3.12) in Eq. (3.9), we obtain

$\displaystyle\Gamma_{\rm superfluid}$
$\displaystyle=\int_{C}d\bm{\ell}\cdot\bigl{(}\bm{\Omega}_{\rm
s}(t,r)\times\bm{r}\bigr{)}=\int_{0}^{r}dr^{\prime}~{}2\pi r^{\prime}\kappa
n(t,r^{\prime}),$ (3.13)

where we choose the integral path $C$ around the edge of the disc with radius
$r$, and $n(t,r)$ is the number density of the vortex lines per unit area.
Noting that only the radial motion of vortex lines changes $n$ due to the
axial symmetry around the rotation axis, we obtain the following conservation
law,

$\displaystyle\frac{\partial n}{\partial t}+\nabla\cdot(nv_{r}\bm{e}_{r})=0,$
(3.14)

where $v_{r}$ is the vortex velocity in the radial direction $\bm{e}_{r}$,
which we call the creep rate. Combining Eqs. (3.13) and (3.14), we obtain the
equation of motion for the superfluid component:

$\displaystyle\frac{\partial\Omega_{\rm s}}{\partial t}$
$\displaystyle=-\left(2\Omega_{\rm s}+r\frac{\partial\Omega_{\rm s}}{\partial
r}\right)\frac{v_{r}}{r}.$ (3.15)

The NS rotation is described by Eqs. (3.6) and (3.15) coupled through a
nonzero $v_{r}$. In other words, by switching off the radial motion of the
vortex line, we have $\partial\Omega_{\rm s}/\partial t=0$. In this case,
$N_{\rm int}$ turns out to be zero, and the two equations of motion are
decoupled.

###  3.2 Dynamics of a vortex line

In the inner crust, a vortex line feels two forces, the pinning force and the
Magnus force. The pinning force arises from the interaction between a vortex
line and a nuclear-like object within the inner crust, resulting in the
pinning of the vortex line to the lattice of nuclei where the energy is
minimized [10]. As long as the pinning force is dominant, the vortex lines
remain attached to the crust and move at the same velocity as the crust
component:

$\displaystyle\bm{v}_{\rm VL}(t,r)=\bm{\Omega}_{\rm c}(t)\times\bm{r},$ (3.16)

where $\bm{v}_{\rm VL}(t,r)$ denotes the velocity of a vortex line.

One way to quantify the pinning force is to compare the energies associated
with different configurations of a vortex line.

Figure 2:  The nuclear pinning and interstitial pinning for the vortex line
pinning configurations. Red cylinders and black spheres represent vortex lines
and nuclei, respectively.

In Fig. 2, we show two possibilities for the pinning configurations, the
nuclear pinning and the interstitial pinning. The difference in the energies
between these two configurations defines the pinning energy,

$\displaystyle E_{\rm pin}\equiv E_{\rm NP}-E_{\rm IP},$ (3.17)

where $E_{\rm NP}$ ($E_{\rm IP}$) denotes the energy of the nuclear
(interstitial) pinning configuration. A nuclear pinning configuration occurs
when the pinning energy is negative, and the vortex line is directly attached
to the nuclei lattice. Conversely, when the pinning force is repulsive, the
vortex line is pinned in the interstitial regions. For a nuclear pinning
configuration, a rather crude estimate of the pinning force per unit length is
then given by

$\displaystyle f_{\rm pin}|_{\rm NP}$ $\displaystyle\simeq\frac{|E_{\rm
pin}|}{\Delta r\Delta L},$ (3.18)

where $\Delta r$ is the distance between the nuclear and interstitial pinning
positions and $\Delta L$ is the distance between the successive pinning sites
along a vortex line. These two quantities are expected to be of the order of
the Wigner-Seitz radius $R_{\rm WS}$, the radius of an imaginary sphere whose
volume is equal to the average volume per nucleus in each region.

Due to this pinning effect, the relative velocity between the superfluid and
the vortex lines is developed,

$\displaystyle\delta\bm{v}\equiv\bm{v}_{\rm s}-\bm{v}_{\rm
VL}=\delta\bm{\Omega}\times\bm{r},$ (3.19)

where we use Eqs. (3.12) and (3.16) to obtain the last expression and
introduce the relative angular velocity,

$\displaystyle\delta\bm{\Omega}$ $\displaystyle\equiv\bm{\Omega}_{\rm
s}-\bm{\Omega}_{\rm c}.$ (3.20)

This velocity difference induces the Magnus force per unit length of a vortex
line,

$\displaystyle\bm{f}_{\rm Mag}$
$\displaystyle=\rho\,(\delta\bm{v})\times\bm{\kappa},$ (3.21)

where $\rho$ is the superfluid density and $\bm{\kappa}$ is the vorticity
vector, which is parallel to the vortex line (hence, parallel to the
rotational axis of the NS) and has the absolute value
$|\bm{\kappa}|\equiv\kappa$.

Figure 3: The Magnus force acting on a vortex line. We define the direction of
$\bm{\kappa}$ as the direction of the right-hand screw.

We see that the Magnus force always acts in the outward direction, as
illustrated in Fig. 3, since the crust component rotates slower than the
superfluid component due to the deceleration by the external torque.

Vortex lines overcome the trapping through thermal fluctuations or quantum
tunneling [49, 50], and start to creep outward. If this creep rate is rapid
enough, the superfluid rotation can smoothly follow the change in the crust
rotation, and the system can reach a steady state. To see if this is the case
for the NSs of our interest, let us briefly review the evaluation of the
vortex creep rate $v_{r}$ following Ref. [50]. In this analysis, the pinning
force is modeled by a periodic potential with a period equal to the span of
the nuclei lattice $(\sim R_{\rm WS})$ and a height equal to the pinning
energy $(\sim|E_{\rm pin}|)$. The Magnus force is considered as a bias that
tilts the periodic potential. Let us introduce the transition rate for a
vortex line to move from a local minimum into the next local minimum as ${\cal
R}_{\rm VC}$. The zero-point frequency in the vicinity of local minima of the
potential $\omega_{0}$ controls ${\cal R}_{\rm VC}$ and is obtained through
the quantization of the vortex system [49, 50]. If the vortex tension is
negligible compared to the pinning force, we obtain666In Ref. [50], the case
where the vortex tension dominates the pinning force is also studied, and the
conclusion of the steady state turns out to remain unchanged.

$\displaystyle\omega_{0}\simeq\frac{\pi\kappa\Lambda}{4R_{\rm WS}^{2}}\simeq
1.2\times 10^{20}\,{\rm s}^{-1}\left(\frac{R_{\rm WS}}{50~{}{\rm
fm}}\right)^{-2}\left(\frac{\Lambda}{2}\right),$ (3.22)

where $\Lambda$ characterizes the vortex tension,
$T_{v}=\rho\kappa^{2}\Lambda/(4\pi)$, and is evaluated as
$2\lesssim\Lambda\lesssim 10$ in the inner crust [51]. The transition rate is
then given by [49, 50]

$\displaystyle{\cal R}_{\rm
VC}\simeq\frac{\omega_{0}}{2\pi}\,e^{-\frac{|E_{\rm pin}|}{k_{\rm B}T_{\rm
eff}}},$ (3.23)

where $T_{\rm eff}$ is defined by

$\displaystyle k_{\rm B}T_{\rm
eff}\equiv\frac{\hbar\omega_{0}}{2}\coth\left(\frac{T_{\rm
q}}{T}\right)\sim\begin{cases}k_{\rm B}T&(T\gg T_{\rm q})\\\
\frac{\hbar\omega_{0}}{2}&(T\ll T_{\rm q})\end{cases},$ (3.24)

with $T$ being the temperature of the inner crust and

$\displaystyle T_{\rm q}\equiv\frac{\hbar\omega_{0}}{2k_{\rm B}}\simeq
3.8\times
10^{8}~{}\mathrm{K}~{}\left(\frac{\omega_{0}}{10^{20}~{}\mathrm{s}^{-1}}\right).$
(3.25)

As can be seen from these expressions, the unpinning occurs predominantly
through thermal activation (quantum tunneling) for $T\gg T_{\rm q}$ ($T\ll
T_{\rm q}$). In particular, for NSs as old as those considered in the
following analysis, their internal temperature is $\lesssim 10^{8}$ K.
Therefore, the vortex creep motion in these NSs is triggered by quantum
tunneling. By using the transition rate in Eq. (3.23), we evaluate the creep
rate as $v_{r}\simeq{\cal R}_{\rm VC}\cdot R_{\rm WS}$.

With the vortex creep rate obtained above, we can estimate the typical
distance at which a vortex line travels through the creep motion during the
lifetime of the NS:

$\displaystyle v_{r}\times t_{\rm age}$ $\displaystyle\simeq R_{\rm
WS}\cdot\frac{\omega_{0}}{2\pi}\cdot e^{-\frac{2|E_{\rm
pin}|}{\hbar\omega_{0}}}\times t_{\rm age}$ $\displaystyle\simeq 160~{}{\rm
km}\times\left(\frac{\omega_{0}}{10^{20}~{}\rm
s^{-1}}\right)\left(\frac{R_{\rm WS}}{50~{}\rm fm}\right)\left(\frac{t_{\rm
age}}{10^{5}~{}\rm yr}\right),$ (3.26)

where we set $|E_{\rm pin}|=1$ MeV as a representative value. We find that
this distance is considerably larger than the crust thickness ($\simeq 1$ km),
implying that the vortex creep motion would reach a steady state in old NSs.

Once the system enters the steady phase, the crust and superfluid components
decelerate at the same rate, i.e.,

$\displaystyle\frac{\partial\Omega_{\rm s}(t,r)}{\partial t}=\dot{\Omega}_{\rm
c}\equiv\dot{\Omega}_{\infty}~{}.$ (3.27)

Note that $\Omega_{\rm c}$ and $|\dot{\Omega}_{\infty}|$ are identified as the
current observed values of the angular velocity and deceleration rate of NSs,
respectively. As a consequence, the relative angular velocity between the
crust and superfluid components becomes independent of time, and this is found
to be fairly close to the critical angular velocity determined by the
condition $f_{\rm pin}=f_{\rm Mag}$ [50]:

$\displaystyle\delta\Omega_{\infty}\simeq\delta\Omega_{\rm cr}~{},$ (3.28)

with

$\displaystyle\delta\Omega_{\rm cr}$
$\displaystyle\equiv|\delta\bm{\Omega}|_{f_{\rm pin}=f_{\rm Mag}}=\frac{f_{\rm
pin}}{\rho\kappa r}.$ (3.29)

###  3.3 Prediction of surface temperature

As vortices move outward, the rotational energy of the superfluid is
dissipated through their frictional interaction with the normal components of
the inner crust, which heats the NS. This heating luminosity is computed as
[12, 15]

$\displaystyle L_{\mathrm{H}}$ $\displaystyle=N_{\rm ext}\Omega_{\rm
c}(t)-\frac{d}{dt}\left[\frac{1}{2}I_{\rm c}\Omega_{\rm
c}^{2}(t)+\frac{1}{2}\int dI_{\rm p}\Omega_{\rm s}^{2}(t,r)\right]$
$\displaystyle=\int dI_{\rm p}\left[\Omega_{\rm c}(t)-\Omega_{\rm
s}(t,r)\right]\frac{\partial\Omega_{\rm s}(t,r)}{\partial t}$
$\displaystyle=\int dI_{\rm p}\,\delta\Omega\left|\frac{\partial\Omega_{\rm
s}(t,r)}{\partial t}\right|~{},$ (3.30)

where we use Eqs. (3.6) and (3.7), and the inertial momenta $I_{\rm p}$ is
integrated over the region where the pinning process efficiently occurs. This
expression is further simplified in the steady state with the condition
(3.27),

$\displaystyle L_{\mathrm{H}}$ $\displaystyle=J|\dot{\Omega}_{\infty}|~{},$
(3.31)

where we define,

$\displaystyle J\equiv\int dI_{\rm p}\delta\Omega_{\infty}~{},$ (3.32)

and $\delta\Omega_{\infty}$ denotes the steady-state value of the relative
angular velocity. As we see, the heating luminosity is proportional to the
current deceleration rate of the pulsar. Note that the value of $J$ can be
estimated from Eqs. (3.27), (3.29), and (3.32) by specifying the value of
$f_{\rm pin}$ and the region of pinning. We will evaluate the proportional
coefficient $J$ in Sec. 4.2.

In particular, if this heating luminosity balances with the photon luminosity,
which we expect to occur for old NSs as we discussed in Sec. 2, the surface
temperature $T_{\rm s}$ can be estimated using Eq. (2.3) as

$\displaystyle T_{\rm s}^{\rm eq}$
$\displaystyle\equiv\left(\frac{J|\dot{\Omega}_{\infty}|}{4\pi R_{\rm
NS}^{2}\sigma_{\rm SB}}\right)^{\frac{1}{4}}$ (3.33) $\displaystyle\simeq
1.0\times
10^{5}~{}\mathrm{K}\left(\frac{J}{10^{43}~{}\mathrm{erg~{}s}}\right)^{\frac{1}{4}}\left(\frac{|\dot{\Omega}_{\infty}|}{10^{-14}~{}\mathrm{s}^{-2}}\right)^{\frac{1}{4}}\left(\frac{R_{\rm
NS}}{11.43~{}\mathrm{km}}\right)^{-\frac{1}{2}}.$ (3.34)

In the steady vortex creep scenario, therefore, the surface temperature is
predicted as a function of $J$, $|\dot{\Omega}_{\infty}|$, and $R_{\rm NS}$,
free from the uncertainty of the initial condition and the subsequent
temperature evolution. We will examine this prediction against observation in
Sec. 5.

##  4 Theoretical approaches for the vortex pinning

To compute the energy dissipation due to vortex creep, we evaluate the
parameter $J$ in Eq. (3.32). We first review the calculation of the pinning
force $f_{\rm pin}$ available in the literature in Sec. 4.1; the values of
$f_{\rm pin}$ which our analysis is based on are summarized in Appendix A.
Then, in Sec. 4.2, we estimate possible ranges of $J$ using the results in
Sec. 4.1, which will be compared with observation in the subsequent section.

###  4.1 Evaluation of pinning force

To evaluate the pinning force, we need to analyze a nucleon many-body system
at high densities, which generically suffers from technical difficulties due
to less-known properties of nuclear interactions. A traditional method of
treating nuclear interactions is to model a form of the interaction and fit it
to experimental data, such as nucleon-nucleon scattering [52, 53, 54]. For
example, the Argonne interaction [54], which we consider in the following
analysis, is a two-body nucleon potential fitted to nucleon scattering data
and deuteron properties. This sort of bare interaction does not include in-
medium effects. The many-body calculation based on bare interaction is
necessary to obtain the properties of nucleon systems. At the same time, it is
still challenging to perform it in general due to, e.g., the strong repulsive
core. An alternative method is to use an effective interaction incorporating
in-medium effect phenomenologically. Skyrme-type interactions [55, 56] are
well-known examples, which consist of contact (zero-range) interactions with
momentum-dependent coefficients. The parameters of the interactions are
determined by fitting to the experimental data of binding energies and radii
of several nuclei. There are many sets of fitting parameters used in the
literature [57, 58, 59], such as SLy4 [60] and SkM* [61]. There are also
finite-range interactions, such as the Gogny interaction [62].

There are several approaches to analyzing the nuclear matter, and the
following two are often used for the calculation of the pinning energy:

  * •

Quantum approach

A standard method to analyze a quantum multi-body system is to calculate the
energy levels of a single particle in the mean field of a self-consistent
potential by solving the corresponding Schrödinger equation. A neutron pairing
interaction is then considered to determine the pairing field as in the BCS
theory. These solutions are obtained via an iterative process such that they
satisfy the self-consistent conditions. This method is called the Hartree-
Fock-Bogoliubov (HFB) method and adopted in Refs. [63, 64, 65, 66].

  * •

Semi-classical approach

The quantum approach based on the HFB method often requires a high
computational cost. To evade this, a semi-classical approach based on the
Thomas-Fermi approximation is also frequently used, where nuclear matter is
regarded as a many-body system of nucleons subject to the Pauli exclusion
principle and moving independently from each other in a mean-field potential.
The energy of the system for a given chemical potential is obtained by the
variational principle.777Strictly speaking, we minimize the modified
Hamiltonian defined by $H^{\prime}=H-\mu N$, where $\mu$ and $N$ denote the
chemical potential and the number of particles, respectively. We also note
that for NSs, the temperature $T$ can be regarded as zero; thus, the free and
internal energy are equivalent. This approach is used in Refs [67, 68, 69, 70,
71, 72].

We can then estimate the pinning force from the pinning energy obtained above.
There are two approaches for this calculation:

  * •

Microscopic calculation

We may estimate the pinning force through Eq. (3.18) for the nuclear pinning
configuration. We refer to this estimation as the microscopic approach since,
as illustrated in the most right window in Fig. 4, it focuses on the
microscopic scale of $\mathcal{O}(R_{\rm WS})$. This approach considers the
interaction between a vortex and the single nucleus in the Wigner-Seitz cell,
and thus the interaction of the vortex with other distant nuclei is neglected.
We obtain the pinning force per unit length by just multiplying the pinning
force per nucleus with the number of nuclei in the unit length along the
vortex line ($\simeq 1/R_{\mathrm{WS}}$).

  * •

Mesoscopic calculation

Vortex lines are much longer than the lattice spacing; therefore, each vortex
line pins onto a large number of nuclei in reality. Such a vortex line does
not align to the crystal axis over its total length in general. The mesoscopic
approach considers this realistic configuration by taking the average of the
force exerted on a vortex over the possible directions of the vortex line with
respect to the crystal lattice. This calculation focuses on the length-scale
$L\sim(10^{2}$–$10^{3})\times R_{\rm WS}$, for which the vortex line can be
regarded as a straight line, as illustrated in the middle window in Fig. 4—we
call this scale the mesoscopic scale. The derived pinning force thus tends to
be smaller than those obtained with the microscopic calculation.

Figure 4: The landscape of a vortex-line configuration at different length-scales.  | Semi-classical | Quantum
---|---|---
Microscopic | Mean field: | Pairing: | Hartree-Fock: | Pairing:
Woods-Saxon | Argonne  | SLy4 | contact force
$f_{\mathrm{pin}}=(1-7)\times 10^{-3}$ [71] |  $f_{\mathrm{pin}}=(3-4)\times 10^{-4}$ [65]
Mesoscopic | Mean field: | Pairing: | Hartree-Fock: | Pairing:
Woods-Saxon | Argonne | SLy4 | contact force
$f_{\mathrm{pin}}=8\times 10^{-7}-4\times 10^{-4}$ [72] |  $f_{\mathrm{pin}}=5\times 10^{-7}-8\times 10^{-5}$ [66]
Table 1:  Calculations of the vortex pinning force considered in this work.
The first row in each column lists the potentials for the mean-field and
pairing interactions used in the evaluations; the second row shows the range
of the pinning force in the inner crust in units of
$\mathrm{MeV}\cdot\mathrm{fm}^{-2}$.

All in all, we have $2\times 2=4$ combinations for the prescription of the
pinning force calculation. In the following discussion, we consider a
representative calculation with a specific choice of nuclear interactions for
each combination, as summarized in Table 1.

For the microscopic semi-classical approach, we consider the calculation in
Ref. [71] with a Woods-Saxon potential for the mean-field potential and the
Argonne interaction for the neutron-neutron pairing interaction. It is found
that the nuclear pinning configuration (see Fig. 2) occurs only in high-
density regions. The pinning force per nuclear pinning site is estimated by
$E_{\mathrm{pin}}/R_{\mathrm{WS}}$ for this configuration in Ref. [71]. To
convert this into the pinning force per unit length $f_{\mathrm{pin}}$, we
multiply this by a factor of $1/(2R_{\mathrm{WS}})$, as cubic cells with the
side length $2R_{\mathrm{WS}}$ are used in the calculation of Ref. [71]. As a
result, we obtain the values of $f_{\mathrm{pin}}=\text{($1$--$7$)}\times
10^{-3}~{}\mathrm{MeV}\cdot\mathrm{fm}^{-2}$, depending on the position in the
inner crust. We list quantities relevant to this calculation in Table 3 in
Appendix A.888To get some ideas about the dependence of this calculation on
the choice of potentials, we note that Ref. [71] also considers the case where
the Gogny interaction is used instead of the Argonne interaction. In this
case, the nuclear pinning occurs in higher density regions, and the pinning
forces tend to be larger by a factor of a few, compared with the calculation
with the Argonne interaction; see Table 4 in Appendix A for this calculation.
Notice that the values of the pinning force shown here should be regarded as a
ballpark estimate. The lower limit ($f_{\mathrm{pin}}=1\times
10^{-3}~{}\mathrm{MeV}\cdot\mathrm{fm}^{-2}$) could have been overestimated
because of the discretization of the positions at which the pinning energy is
estimated; as seen in Table 3, there is a density region of $\rho\sim
10^{13}~{}\mathrm{g}\cdot\mathrm{cm}^{-3}$ around which
$E_{\mathrm{pin}}\simeq 0$, leading to a very small pinning force if we use
Eq. (3.18). The upper limit could also be underestimated since the pinning
force obtained by this equation corresponds to the average taken over the
distance between the nuclear center and the interstitial position.

For the microscopic quantum approach, we consider the results given in Ref.
[65], where the SLy4 Skyrme interaction is used for the Hartree-Fock
calculation and a density-dependent contact interaction is used for the
neutron-neutron pairing interaction. The parameters of the contact interaction
are discussed in Ref. [73]. Table 5 summarises the relevant quantities for
this calculation. In contrast to the semi-classical analysis [71], nuclear
pinning occurs in lower-density regions in this case. However, it also occurs
in the highest-density regions (see Table 5 in Appendix A).999The qualitative
feature described here does depend on the choice of the nuclear interactions.
As shown in Ref. [65], if we use the SkM* interaction instead of SLy4, the
nuclear pinning never occurs at high densities.  The values of the pinning
force estimated from the pinning energy and the Wigner-Seitz radius as in the
semi-classical calculation are $f_{\mathrm{pin}}=\text{($3$--$4$)}\times
10^{-4}~{}\mathrm{MeV}\cdot\mathrm{fm}^{-2}$. See Ref. [65] for detailed
discussions regarding the difference between the semi-classical and quantum
results.

The semi-classical mesoscopic calculation is given in Ref. [72], where the
nuclear potentials are taken to be the same as in the semi-classical
microscopic calculation in Table 1. The resultant values of the pinning force
are found to be $f_{\mathrm{pin}}=8\times 10^{-7}$ – $4\times
10^{-4}~{}\mathrm{MeV}\cdot\mathrm{fm}^{-2}$, which are summarized in Table 6.
We see that these values are considerably smaller than those for the semi-
classical microscopic calculation in Ref. [71] due to the averaging over the
vortex-line directions.

For the quantum mesoscopic calculation, we consider the result given in Ref.
[66], where the SLy4 interaction and a contact interaction are used for the
Hartree-Fock calculation and pairing interactions, respectively, as in the
quantum microscopic calculation in Table 1. As summarized in Table 7, the
pinning force is found to be in the range $f_{\mathrm{pin}}=5\times
10^{-7}-8\times 10^{-5}~{}\mathrm{MeV}\cdot\mathrm{fm}^{-2}$.101010If we use
the SkM* potential instead of SLy4, we obtain slightly smaller values of
$f_{\mathrm{pin}}$, as also shown in Table 7.  We again find that these values
are much smaller than those obtained with the quantum microscopic approach.

Before concluding this section, we note that we can also calculate the pinning
force with a three-dimensional dynamical simulation of a vortex [74, 75, 76,
77].111111See also Ref. [78] for a kinetic approach. Such a calculation tends
to be costly; thus, a certain degree of simplification is usually required for
the moment. Besides, the current evaluation is limited to a few benchmark
values of density. The estimated values of the pinning force are consistent
with the above estimates.

###  4.2 Theoretical evaluation of $J$

For the evaluation of $J$ in Eq. (3.32), it is convenient to change the
coordinate from cylindrical coordinates $(r,\varphi,z)$ to spherical
coordinates $(R,\theta,\phi)$:

$\displaystyle J$ $\displaystyle=\int_{\rm pin}dr\,dz\,d\varphi\,\rho
r^{3}\cdot\frac{f_{\mathrm{pin}}}{\rho\kappa r}\simeq\int_{R_{\rm in}}^{R_{\rm
out}}dR\,d\theta\,d\phi\,R^{3}\sin^{2}\theta\cdot\frac{f_{\rm
pin}}{\kappa}\,,$ (4.35)

where we approximate $\delta\Omega_{\infty}$ in Eq. (3.32) by
$\delta\Omega_{\mathrm{cr}}$ in Eq. (3.29), which holds with good accuracy in
the situation of our interest [50] as we mentioned in Sec. 3.2. We perform the
integral over the range $[R_{\rm in},R_{\rm out}]$ where the pinning force is
evaluated. We use the Akmal-Pandharipande-Ravenhall (APR) [79] equation of
state to determine the NS core size and the equation of state tabulated in
Crust_EOS_Cat_HZD-NV.dat in NSCool [80] based on Refs. [81, 82] to determine
the density distribution in the crust.

For the evaluation method of pinning force, we focus on the mesoscopic
approach shown in Table 1, since for the microscopic calculation, the pinning
force is obtained only in a limited region in the crust, as can be seen in
Table 3–5. Considering that the evaluation of the pinning force suffers from
large uncertainty depending on calculation methods, we make the following
crude approximation in the calculation of the above integral—we neglect the
density dependence of $f_{\mathrm{pin}}$ and fix it to a value in the range
shown in Table 1. We thus obtain a range of $J$ accordingly, which we regard
as the uncertainty of this pinning force estimation. As a result, we obtain
$J=3.9\times 10^{40}-1.9\times 10^{43}~{}\mathrm{erg}\cdot\mathrm{s}$ for the
semi-classic mesoscopic calculation and $J=1.7\times 10^{40}-2.7\times
10^{42}~{}\mathrm{erg}\cdot\mathrm{s}$ for the quantum mesoscopic calculation.

##  5 Vortex creep heating vs. observation

No. | Type | Name | $\log_{10}t_{\rm sd}$ | $\log_{10}t_{\rm kin}$ | $\log_{10}|\dot{\Omega}|$ | $\log_{10}T_{\rm s}$ | $\log_{10}J_{\rm obs}$
---|---|---|---|---|---|---|---
|  |  | [yr] | [yr] |  [s-2] | [K] | [erg s]
1. | [Y] | PSR B1706-44 | $4.2$ | — | $-10.3$ |  $5.68$–$6.34$ [83] |  $41.9$–$44.6$
2. | [O] | PSR J1740+1000 | $5.1$ | — | $-11.2$ |  $5.89$ [84] | $43.8$
3. | [Y] | PSR B2334+61 | $4.6$ | — | $-11.3$ |  $5.76$ [85] | $43.3$
4. | [O] | PSR B0656+14 | $5.0$ | — | $-11.6$ |  $5.87$ [86] | $44.1$
5. | [O] | PSR J0633+1746 | $5.5$ | — | $-11.9$ |  $5.71$ [87] | $43.7$
6. | [Y] | PSR J0538+2817 | $5.8$ |  $4.60^{+0.18}_{-0.30}$[88] | $-11.9$ |  $6.02$ [88] | $45.0$
7. | [O] | PSR B1055-52 | $5.7$ | — | $-12.0$ |  $5.81$ [89] | $45.1$
8. | [X] | RX J1605.3+3249 | $4.5$ |  $5.66^{+0.04}_{-0.07}$ [90] | $-12.1$ |  $5.86$ [91] | $44.5$
9. | [O] | PSR J2043+2740 | $6.1$ | — | $-12.1$ |  $<5.95$ [92] | $<44.8$
10. | [O] | PSR J1741-2054 | $5.6$ | — | $-12.2$ |  $5.85$ [93] | $44.6$
11. | [O] | PSR J0357+3205 | $5.7$ | — | $-12.4$ |  $5.62$ [94] | $43.8$
12. | [O] | PSR B0950+08 | $7.2$ | — | $-13.6$ |  $4.78$–$5.08$ [29] |  $41.7$–$42.9$
13. | [X] | RX J0420.0-5022 | $6.3$ | — | $-13.8$ |  $5.74$ [95] | $45.8$
14. | [M] | PSR J0437-4715 | $9.20$ | — | $-14.0$ |  $5.54$ [25] | $45.1$
15. | [X] | RX J1308.6+2127 | $6.2$ |  $5.74^{+0.16}_{-0.26}$ [96] | $-14.2$ |  $6.08$ [97] | $47.5$
16. | [X] | RX J0720.4-3125 | $6.3$ |  $5.93^{+0.07}_{-0.26}$ [98] | $-14.2$ |  $6.02$ [99] | $47.3$
17. | [M] | PSR J2124-3358 | $9.58$ | — | $-14.3$ |  $4.70$–$5.32$ [26] |  $42.0$–$44.5$
18. | [X] | RX J1856.5-3754 | $6.6$ |  $5.66^{+0.04}_{-0.05}$[98] | $-14.4$ |  $5.65$ [100] | $46.0$
19. | [X] | RX J2143.0+0654 | $6.6$ | — | $-14.5$ |  $5.67$–$6.06$ [101, 102] |  $46.1$–$47.8$
20. | [X] | RX J0806.4-4123 | $6.5$ | — | $-14.6$ |  $6.01$ [103] | $47.6$
21. | [O] | PSR J0108-1431 | $8.3$ | — | $-15.1$ |  $4.43$–$4.74$ [28] |  $41.8$–$43.1$
22. | [O] | PSR J2144-3933 | $8.4$ | — | $-16.4$ |  $<4.62$ [104] | $<43.8$

Table 2:  The data of the NSs considered in this paper. We classify them into
four types—ordinary pulsars younger than $10^{5}~{}\mathrm{yrs}$ [Y], ordinary
pulsars older than $10^{5}~{}\mathrm{yrs}$ [O], XDINSs [X], and millisecond
pulsars [M]. The values of $t_{\rm sd}$ and $\dot{\Omega}$ without references
are computed from the data given in the ATNF pulsar catalogue [105, 106]. The
value of $J_{\rm obs}$ is evaluated as in Eq. (5.36) with $R_{\rm
NS}=11.43~{}\mathrm{km}$.

We now compare the prediction of the vortex-creep heating mechanism with
observation. For this purpose, it is useful to calculate the following
quantity for each NS:121212We neglect the gravitational redshift factor since
its effect is within the $\mathcal{O}(1)$ uncertainty of $J$ discussed below.

$\displaystyle J_{\rm obs}$ $\displaystyle\equiv\frac{4\pi R_{\rm
NS}^{2}\sigma_{\rm SB}T_{\rm s}^{4}}{|\dot{\Omega}|}~{}.$ (5.36)

As evident from Eqs.(2.2), (2.3), and (3.31), this corresponds to the $J$
parameter inferred from the observation of each NS. This inference assumes the
steady creeping of vortices (discussed in Sec.3.2) and the balance between
vortex-creep heating luminosity and photon cooling luminosity, which we expect
to hold if the NS is older than $\sim 10^{5}$ years. Since NSs are comparable
in size and mass, we expect that $J_{\rm obs}$ is also roughly equal (up to a
factor of $\mathcal{O}(1)$) for every NS. We test this expectation by using
the data of $J_{\rm obs}$ for old NSs. We also compare the values of
$J_{\mathrm{obs}}$ with the theoretical computations given in Sec. 4.2.

In Table 2, we list the values of $J_{\mathrm{obs}}$ for the NSs we consider
in this paper. We select isolated NSs older than $10^{4}$ yrs. In evaluating
$J_{\mathrm{obs}}$, we have just assumed $R_{\rm NS}=11.43~{}\mathrm{km}$ for
all NSs, as the radius is poorly known for most NSs; we keep in mind that this
may introduce an $\mathcal{O}(1)$ error in the determination of
$J_{\mathrm{obs}}$. We also show the age, surface temperature, and
$\dot{\Omega}=2\pi\dot{P}/P^{2}$ of the NSs, ($P$ and $\dot{P}$ are the period
and its time derivative, respectively). Regarding the NS age, we use the
kinetic age $t_{\rm kin}$ if available. Otherwise, we use the spin-down age
$t_{\rm sd}=P/(2\dot{P})$. We calculate $t_{\rm sd}$ and $\dot{\Omega}$ from
$P$ and $\dot{P}$ given in the Australia Telescope National Facility (ATNF)
pulsar catalogue [105, 106]. Notice that the surface temperatures of some of
the old NSs in this table are much higher than the predicted temperature in
the standard NS cooling scenario [30, 31, 32, 33, 34, 35].

It is important to note that not all NSs listed in Table 2 are useful for
testing the vortex-creep heating. As we have discussed in Sec. 2, photon
emission becomes the dominant cooling source for NSs older than $\sim 10^{5}$
years. For younger NSs, this may not be the case, so
$L^{\infty}_{\gamma}\lesssim L^{\infty}_{\rm H}$ instead of (2.3), for which
the values of $J_{\mathrm{obs}}$ in Table 2 may be underestimated. To
distinguish such young NSs from others, we indicate them by the type [Y] in
the table. Another class of NSs that are inappropriate for our test is the
X-ray Dim Isolated NSs (XDINSs). These NSs are considered to be descendants of
magnetars [103, 95] that experienced the decay of magnetic fields before. This
may make these NSs hotter than ordinary NSs of the same age [107, 95, 108,
109], resulting in an overestimate of $J_{\mathrm{obs}}$. We denote these NSs
by the type [X]. The rest of the NSs, which we use for the test of the vortex-
creep heating, are classified into old ordinary pulsars [O] and millisecond
pulsars [M].

The uncertainty in the determination of $J_{\mathrm{obs}}$ stems mainly from
that in the surface temperature, which is significant due to its quartic
dependence on $T_{\rm s}$. Generically speaking, it is very difficult to
identify all of the sources of uncertainties in the measurement of the NS
surface temperature, and it is often the case that the error shown in the
literature is only a part of them, such as those from the spectrum fitting,
the determination of the distance and/or radius of the NS, and so on. At
present, it is fair to say that the NS temperature measurement typically
suffers from $\mathcal{O}(1)$ uncertainty, as can be seen in, e.g., Ref.
[110]. Motivated by this, we include a factor of two uncertainty in $T_{\rm
s}$ for the stars in Table 2 for which only the central value is presented.
For the other stars, we describe our prescription for the error estimation in
Appendix B. We have checked that the errors thus obtained are similar to or
more conservative than those adopted in Ref. [110].

Figure 5:  The values of the $J$ parameter obtained from the observation. The
grey triangles, blue circles, green inverse triangles, and orange stars
correspond to the young ordinary pulsars ([Y]), the old ordinary pulsars
([O]), the XDINSs ([X]), and the millisecond pulsars ([M]), respectively. The
points with an arrow indicate upper limits. The red shaded region shows
observationally favored range, $J\simeq
10^{42.9\text{--}43.8}~{}\mathrm{erg}\cdot\mathrm{s}$. For comparison, we also
show the values of $J$ estimated with the mesoscopic calculations by black
bars.

In Fig. 5, we show the range of $J_{\mathrm{obs}}$ estimated as described
above for each NS listed in Table 2. The grey triangles, green inverse
triangles, blue circles, and orange stars represent the young ordinary pulsars
([Y]), the XDINSs ([X]), the old ordinary pulsars ([O]), and the millisecond
pulsars ([M]), respectively. The points with an arrow indicate that we only
have an upper limit on $J_{\rm obs}$ for those NSs. Recall that we are
concerned only with the NSs represented by the blue [O] and orange [M] points.
It is found that the estimated values of $J_{\mathrm{obs}}$ for these NSs are
in the same ballpark, $J\sim 10^{43}~{}\mathrm{erg}\cdot\mathrm{s}$, even
though their $|\dot{\Omega}|$’s distribute over orders of magnitude. This is
in good agreement with the prediction of the vortex-creep heating mechanism.
On the other hand, $J_{\mathrm{obs}}$ for the green points (XDINSs [X]) tend
to be larger than this, as expected.

We also show the theoretical estimations given in Sec. 4.2 in the upper panel
of Fig. 5. We see that the semi-classical mesoscopic calculation is consistent
with the observation, given that this theoretical estimation suffers from a
NS-dependent uncertainty of $\mathcal{O}(1)$ coming from the integration in
Eq. (4.35), in addition to that from the estimation of $f_{\mathrm{pin}}$. The
quantum mesoscopic calculation can explain some of the points with a small
$J_{\mathrm{obs}}$, but they are not large enough to explain, e.g., that of
J0437-4715. However, we note that this theoretical estimation is still allowed
by the observations since it just results in a lower heating luminosity than
the observed one. If this is the case, the vortex-creep heating may operate
but there exists another heating mechanism that dominates the vortex-creep
heating, such as the rotochemical heating [111, 112, 113, 114, 115, 116, 117,
118, 119].

It is premature to establish the existence of the vortex-creep heating, as
well as to conclude if an extra heating mechanism is required to be present.
To that end, we need to accumulate more data on the surface temperature of old
NSs with high accuracy, which we anticipate to be provided by future optical,
UV, and X-ray observations.131313See, for instance, Ref. [120]. Nevertheless,
obtaining a current compilation of the value of $J$ suggested by the
observation is intriguing. Considering intrinsic $\mathcal{O}(1)$ uncertainty
in $J$, we determine its rough range by requiring that it covers the range
suggested by B0950+08, which favors the smallest value, and satisfies the
upper limit set by J2144-3944 based on non-observation of thermal flux. This
yields

$\displaystyle J\simeq 10^{42.9-43.8}~{}\mathrm{erg}\cdot\mathrm{s}~{},$
(5.37)

which we show as the red band in Fig. 5.

(a) Ordinary pulsars  (b) Millisecond pulsars

Figure 6: The evolution of NS surface temperature with (without) the vortex
creep heating effect in the red band (black dashed line). The band corresponds
to the range of $J$ in Eq. (5.37). The dots with error bars show the observed
temperatures, presented in Table 2, with the same colors as in Fig. 5.

Finally, in Fig. 6, we show the evolution of NS surface temperature with
(without) the vortex creep heating effect in the red band (black dashed
line).141414In these plots, we use the APR equation of state [79] for a NS
mass of $1.4M_{\odot}$ to calculate the NS structure. For Cooper pairing gap
models, we use the SFB model [121] for the neutron singlet pairing, the model
“b” in Ref. [33] for the neutron triplet pairing, and the CCDK model [122] for
the proton singlet pairing. The temperature evolution at late times scarcely
depends on the choice of these models.  The band corresponds to the range of
$J$ in Eq. (5.37). The dots with error bars show the observed temperatures in
Table 2, with the same colors as in Fig. 5. In Fig. 6(a), we take
$P\dot{P}=10^{-15}$ s and the initial period $P_{0}=10~{}\mathrm{ms}$ to
calculate $|\dot{\Omega}(t)|$, where we assume that the external torque is
dominated by magnetic dipole radiation.151515In this case, we have
$\dot{\Omega}\propto-\Omega^{3}$, i.e., $P\dot{P}=\text{constant}$, and by
solving this we obtain
$\displaystyle|\dot{\Omega}|(t)=\frac{\pi}{\sqrt{2P\dot{P}}}\left[t+t_{\mathrm{sd},0}\right]^{-3/2},$
with $t_{\mathrm{sd},0}\equiv P_{0}^{2}/(2P\dot{P})$. For the choice of
parameters in Fig. 6(a) and Fig. 6(b), $t_{\mathrm{sd},0}\simeq 2\times
10^{3}$ and $5\times 10^{7}$ years, respectively. The value of $P\dot{P}$ is
related to the surface magnetic flux density $B_{s}$. In the ATNF pulsar
catalogue [105, 106], $B_{s}=3.2\times 10^{19}(P\dot{P})^{1/2}$ G is used for
this relation, with which we have $B_{s}\simeq 1.0\times 10^{12}$ G and
$5.8\times 10^{8}$ G in Fig. 6(a) and Fig. 6(b), respectively. These values
are typical for ordinary pulsars. Nevertheless, we note that
$|\dot{\Omega}(t)|$ obtained with these parameters do not exactly agree with
the observed values of $|\dot{\Omega}|$ in Table 2, so the data points shown
in this figure should be regarded as just an eye guide. In Fig. 6(b), we set
$P\dot{P}=3.3\times 10^{-22}$ s, which is the observed value for PSR
J2124-3358, and $P_{0}=1~{}\mathrm{ms}$. As we see in these plots, the
predicted temperature with the vortex creep heating starts to deviate from
that in the standard cooling scenario at $t\sim 10^{5}$ years and remains high
enough at later times to be compatible with the observed data.

##  6 Conclusion and discussion

We have revisited the vortex-creep heating mechanism in light of recent
observations of old warm NSs. As we have seen, this heating mechanism gives a
characteristic prediction that the heating luminosity is proportional to
$|\dot{\Omega}|$, with the proportional constant $J$ having an almost
universal value over NSs since the NS structure and the vortex-nuclear
interactions determine it. We have found that this prediction agrees with the
observational data of old NSs, with the favored range of $J$ in the same
ballpark as the theoretical calculations.

Notice that the scenario where vortex creep heating dominates all NSs can
readily be overturned if we discover a NS having $J$ much smaller than those
presented in Fig. 5. On the other hand, if we find a NS with a larger $J$, we
can disfavor our scenario only after excluding the existence of other heating
sources specific to this NS, such as accretion from its environment.

It is possible that other heating mechanisms also work in old NSs. Indeed, we
have already considered potential heating caused by the decay of magnetic
fields in XDINSs (see Sec. 5), and we have not used these NSs in our test of
the vortex-creep heating mechanism for this reason. Another heating mechanism
that may operate without relying on exotic phenomena is provided by the out-
of-equilibrium beta processes, which is dubbed rotochemical heating [111, 112,
113, 114, 115, 116, 117, 118, 119]. It is known that this rotochemical heating
mechanism can increase the surface temperature of old NSs up to $\sim 10^{6}$
K. Thus, its heating luminosity can be comparable to or even dominate the
vortex-creep heating. It would be worthwhile to study the vortex creep heating
in the presence of rotochemical heating and compare its prediction with the
temperature observations of old warm NSs. The NS heating caused by the
accretion of dark matter particles is also widely discussed in the literature
[123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137,
138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152,
153, 154, 155]. In this case, the surface temperature at late times is
predicted to be $\text{a few}\times 10^{3}$ K, and thus this heating effect
could be concealed by the vortex heating mechanism, making it improbable to
observe the dark matter signature through the temperature observation of old
NSs. A detailed study of this issue will be given in the forthcoming paper
[156].

## Acknowledgments

MF thanks Kazuyuki Sekizawa and Tomoya Naito for the fruitful discussion about
the current situation and the possible future directions in the pinning force
evaluation. This work is supported in part by the Collaborative Research
Center SFB1258 and by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy - EXC-2094 -
390783311 [MF], JSPS Core-to-Core Program (No. JPJSCCA 20200002 [MF]), the
Grant-in-Aid for Innovative Areas (No.19H05810 [KH and MRQ], No.19H05802 [KH],
No.18H05542 [MF and NN]), Scientific Research B (No.20H01897 [KH and NN]),
Young Scientists (No.21K13916 [NN]).

##  Appendix A Pinning force

In this appendix, we show the density dependence of the quantities relevant to
the pinning force calculations discussed in Sec. 4.1. Tables 3 and 4 are for
the microscopic semi-classical approach in Ref. [71], where the Argonne and
Gogny interactions are used for the nuclear pairing interaction, respectively.
The element corresponding to the cell nuclear composition and Wigner-Seitz
radius $R_{\rm WS}$ are derived for each baryon density $\rho$ in Ref. [81].
In Ref. [71], the pinning force is evaluated only for the nuclear pinning
configuration, and thus we show the values of $E_{\mathrm{pin}}$ and
$f_{\mathrm{pin}}=|E_{\mathrm{pin}}|/(2R_{\mathrm{WS}}^{2})$ only for this
case. The labels of NP and IP in the last column indicate the nuclear and
interstitial pinnings, respectively.

zone | Element | $\rho$ | $R_{\rm WS}$ | $\xi$ | $E_{\rm pin}$ | $f_{\rm pin}$ | config
---|---|---|---|---|---|---|---
|  | [g cm-3] | [fm] | [fm] | [MeV] | [MeV fm-2] |
1 | ${}^{320}_{40}\mathrm{Zr}$ | $1.5\times 10^{12}$ | $44.0$ | $7.02$ | — | — | IP
2 | ${}^{1100}_{50}\mathrm{Sn}$ | $9.6\times 10^{12}$ | $35.5$ | $4.34$ | — | — | IP
3 | ${}^{1800}_{50}\mathrm{Sn}$ | $3.4\times 10^{13}$ | $27.0$ | $8.54$ | $-5.2$ | $0.0036$ | NP
4 | ${}^{1500}_{40}\mathrm{Zr}$ | $7.8\times 10^{13}$ | $19.4$ | $11.71$ | $-5.1$ | $0.0068$ | NP
5 | ${}^{982}_{32}\mathrm{Ge}$ | $1.3\times 10^{14}$ | $13.8$ | $8.62$ | $-0.4$ | $0.0011$ | NP
Table 3: Quantities relevant to the pinning force calculation obtained with the microscopic semi-classical approach in Ref. [71], where the Argonne potential is used for the nuclear pairing interaction. $\rho$, $R_{\mathrm{WS}}$, and $\xi$ are the mass density, Wigner-Seitz radius, and coherence length, respectively.  zone | Element | $\rho$ | $R_{\rm WS}$ | $\xi$ | $E_{\rm pin}$ | $f_{\rm pin}$ | config
---|---|---|---|---|---|---|---
|  | [g cm-3] | [fm] | [fm] | [MeV] | [MeV fm-2] |
1 | ${}^{320}_{40}\mathrm{Zr}$ | $1.5\times 10^{12}$ | $44.0$ | $7.76$ | — | — | IP
2 | ${}^{1100}_{50}\mathrm{Sn}$ | $9.6\times 10^{12}$ | $35.5$ | $4.07$ | — | — | IP
3 | ${}^{1800}_{50}\mathrm{Sn}$ | $3.4\times 10^{13}$ | $27.0$ | $3.93$ | — | — | IP
4 | ${}^{1500}_{40}\mathrm{Zr}$ | $7.8\times 10^{13}$ | $19.4$ | $7.78$ | $-7.5$ | $0.010$ | NP
5 | ${}^{982}_{32}\mathrm{Ge}$ | $1.3\times 10^{14}$ | $13.8$ | $8.62$ | $-5.9$ | $0.015$ | NP
Table 4: Quantities relevant to the pinning force calculation obtained with the microscopic semi-classical approach in Ref. [71], where the Gogny potential is used for the nuclear pairing interaction.  zone | Element | $\rho$ | $n$ | $R_{\rm WS}$ | $\xi$ | $E_{\rm pin}$ | $f_{\rm pin}$ | config
---|---|---|---|---|---|---|---|---
|  | [g cm-3] | [fm-3] | [fm] | [fm] | [MeV] | [MeV fm-2] |
1a | ${}^{320}_{40}\mathrm{Zr}$ | $1.7\times 10^{12}$ | $0.001$ | $43.3$ | $4.43$ | $-1.08$ | $0.00029$ | NP
1b | ${}^{320}_{40}\mathrm{Zr}$ | $3.4\times 10^{12}$ | $0.002$ | $40.0$ | $4.21$ | $-1.20$ | $0.00038$ | NP
1c | ${}^{320}_{40}\mathrm{Zr}$ | $6.7\times 10^{12}$ | $0.004$ | $36.9$ | $3.93$ | — | — | IP
2a | ${}^{1100}_{50}\mathrm{Sn}$ | $1.3\times 10^{13}$ | $0.008$ | $33.0$ | $4.04$ | — | — | IP
2b | ${}^{1100}_{50}\mathrm{Sn}$ | $1.8\times 10^{13}$ | $0.011$ | $31.0$ | $4.12$ | — | — | IP
2c | ${}^{1100}_{50}\mathrm{Sn}$ | $2.8\times 10^{13}$ | $0.017$ | $28.0$ | $4.70$ | — | — | IP
3a | ${}^{1800}_{50}\mathrm{Sn}$ | $4.3\times 10^{13}$ | $0.026$ | $24.5$ | $6.05$ | — | — | IP
3b | ${}^{1800}_{50}\mathrm{Sn}$ | $6.2\times 10^{13}$ | $0.037$ | $21.4$ | $8.75$ | $-0.41$ | $0.00045$ | NP
Table 5:  Relevant quantities for the pinning force calculation obtained with
the microscopic quantum approach in Ref. [65].

The results for the microscopic quantum approach in Ref. [65] are summarized
in Table 5, where the SLy4 Skyrme interaction is used for the mean-field
interaction. The Wigner-Seitz radius $R_{\rm WS}$ shown in this table is
interpolated from the plot in Ref. [81]. For the evaluation of
$f_{\mathrm{pin}}$, we again use the formula
$f_{\mathrm{pin}}=|E_{\mathrm{pin}}|/(2R_{\mathrm{WS}}^{2})$ and show the
values only for the nuclear pinning, just for easy comparison with the semi-
classical calculations.

zone |  |  |  $f_{\rm pin}$ [MeV fm-2] |  |
---|---|---|---|---|---
$L=100R_{\rm WS}$ | $L=500R_{\rm WS}$ | $L=1000R_{\rm WS}$ | $L=2500R_{\rm WS}$ | $L=5000R_{\rm WS}$
1 | $7.63\times 10^{-6}$ | $2.29\times 10^{-6}$ | $1.49\times 10^{-6}$ | $9.30\times 10^{-7}$ | $7.68\times 10^{-7}$
2 | $2.12\times 10^{-5}$ | $6.34\times 10^{-6}$ | $4.06\times 10^{-6}$ | $2.62\times 10^{-6}$ | $2.12\times 10^{-6}$
3 | $1.43\times 10^{-4}$ | $4.49\times 10^{-5}$ | $2.74\times 10^{-5}$ | $1.54\times 10^{-5}$ | $1.14\times 10^{-5}$
4 | $3.84\times 10^{-4}$ | $1.10\times 10^{-4}$ | $6.89\times 10^{-5}$ | $4.23\times 10^{-5}$ | $3.32\times 10^{-5}$
5 | $7.85\times 10^{-5}$ | $2.71\times 10^{-5}$ | $1.89\times 10^{-5}$ | $1.36\times 10^{-5}$ | $1.12\times 10^{-5}$
Table 6:  The pinning force obtained in the semi-classical mesoscopic approach
for different values of $L$ over which the forces exerted on a vortex are
integrated [72]. The zones correspond to those in Table 3.

In Table 6, we show the pinning force obtained in the semi-classical
mesoscopic approach for different values of $L$ over which the forces exerted
on a vortex are integrated [72]. The zone numbers in the first column
correspond to those in Table 3. We show the calculation in which the reduction
of the pairing gap due to the polarization effects in the nuclear matter is
not included, corresponding to the choice of the reduction factor $\beta=1$
introduced in Ref. [72]. Because of the averaging procedure, we find that a
larger $L$ results in a smaller value of $f_{\mathrm{pin}}$.

model | zone | $E_{\rm pin}$ | config |  |  $f_{\rm pin}$ [MeV fm-2] |
---|---|---|---|---|---|---
|  | [MeV] |  | $L=1000R_{\rm WS}$ | $L=2500R_{\rm WS}$ | $L=5000R_{\rm WS}$
SLy4 | 1a | $-0.72$ | NP | $1.39\times 10^{-6}$ | $7.38\times 10^{-7}$ | $5.40\times 10^{-7}$
1b | $-0.91$ | NP | $1.97\times 10^{-6}$ | $1.08\times 10^{-6}$ | $8.00\times 10^{-7}$
1c | $-0.89$ | NP | $2.20\times 10^{-6}$ | $1.19\times 10^{-6}$ | $8.74\times 10^{-7}$
2a | $2.73$ | IP | $5.61\times 10^{-6}$ | $3.72\times 10^{-6}$ | $2.95\times 10^{-6}$
2b | $3.01$ | IP | $7.52\times 10^{-6}$ | $5.03\times 10^{-6}$ | $4.00\times 10^{-6}$
2c | $10.00$ | IP | $1.47\times 10^{-5}$ | $1.01\times 10^{-5}$ | $8.08\times 10^{-6}$
3a | $11.78$ | IP | $3.25\times 10^{-5}$ | $2.31\times 10^{-5}$ | $1.88\times 10^{-5}$
3b | $9.85$ | IP | $8.47\times 10^{-5}$ | $6.41\times 10^{-5}$ | $5.31\times 10^{-5}$
SkM* | 1a | $-0.72$ | NP | $3.61\times 10^{-7}$ | $1.60\times 10^{-7}$ | $9.50\times 10^{-8}$
1b | $-0.91$ | NP | $2.20\times 10^{-7}$ | $8.72\times 10^{-8}$ | $4.08\times 10^{-8}$
1c | $-0.89$ | NP | $2.83\times 10^{-6}$ | $1.82\times 10^{-6}$ | $1.44\times 10^{-6}$
2a | $2.73$ | IP | $5.68\times 10^{-6}$ | $3.73\times 10^{-6}$ | $2.93\times 10^{-6}$
2b | $3.01$ | IP | $7.58\times 10^{-6}$ | $5.01\times 10^{-6}$ | $4.00\times 10^{-6}$
2c | $10.00$ | IP | $1.25\times 10^{-5}$ | $8.50\times 10^{-6}$ | $6.76\times 10^{-6}$
3a | $11.78$ | IP | $2.54\times 10^{-5}$ | $1.80\times 10^{-5}$ | $1.45\times 10^{-5}$
3b | $9.85$ | IP | $8.00\times 10^{-5}$ | $5.93\times 10^{-5}$ | $4.79\times 10^{-5}$
Table 7:  The pinning force obtained in the quantum mesoscopic calculation
where the Skytme interactions, SLy4 and SkM*, are used for the mean-field
potential [66]. The zones correspond to those in Table 5.

Table 7 shows the pinning force obtained in the quantum mesoscopic calculation
where the SLy4 and SkM* Skytme interactions are used for the Hartree-Fock
calculation [66]. The zones correspond to those in Table 5. We again show the
results obtained without including the polarization effect, i.e., $\beta=1$ as
in the previous case.

Figure 7: The values of $f_{\rm pin}$ given in the tables in this appendix
against the density $\rho$. The filled and opened markers correspond to the
nuclear and interstitial pinnings.

Finally, we plot the values of the pinning force for each density region in
Fig. 7. The filled and opened markers correspond to the nuclear and
interstitial pinnings. As we see, the values of $f_{\rm pin}$ are distributed
in the range $10^{-8}$–$10^{-2}$ MeV $\cdot$ fm-2, depending on the evaluation
scheme and the selected nuclear potential.

##  Appendix B Selection criteria of NS data

We explain how we choose the range of uncertainty of $T_{\rm s}$ for each NS
shown in Table 2.

  * •

No. 1, PSR B1706-4: In Ref. [83], the X-ray data of PSR B1706-44 obtained in
XMM-Newton is fitted by the BB (blackbody), BB+PL (power law), and
atmosphere+PL models, and only BB+PL and atmosphere+PL models result in
acceptable $\chi^{2}$ values. The atmosphere model includes the light-element
NS atmosphere (e.g. dominated by Hydrogen) and shows large Wien excesses in
the high-energy region. Therefore, atmosphere+PL models tend to favor lower
temperature and larger radius for the emitting area. We selected the minimum
and maximum among the BB+PL and atmosphere+PL models,
$T^{\infty}=(0.48-2.2)\times 10^{6}~{}\mathrm{K}$, to include the uncertainty
coming from the choice of fitting models.

  * •

No. 9, PSR J2043+2740: Ref. [157] studied the XMM-Newton data of PSR
J2043+2740. Using the BB + PL model, the upper bound is derived as
$T^{\infty}_{s}<6.27\times 10^{5}~{}\mathrm{K}$, where $R_{\rm
NS}=10~{}\mathrm{km}$ is assumed for the emission radius. On the other hand,
Ref. [158] also fitted the X-ray data of the XMM-Newton and obtained an even
higher BB temperature, $T^{\infty}_{s}\simeq 9\times 10^{5}~{}\mathrm{K}$ with
the radiation radius $R^{\infty}\simeq
2~{}\mathrm{km}~{}\text{\cite[cite]{[\@@bibref{Number}{Zavlin:2004wt}{}{}]}}$.
Although the fitted radius is smaller than the expected NS radius, it is too
large to be interpreted as the magnetic cap radius. Thus, we can not exclude
the possibility that this BB temperature corresponds to the emission from the
NS surface. To evaluate $T_{\rm s}$ conservatively, we chose the highest value
of the BB temperature as an upper bound.

  * •

No. 12, PSR B0950+08: In Ref. [29], the optical-UV flux of PSR B0950+08
obtained in the Hubble Space Telescope (HST) far-UV (FUV) detector is
analyzed. The best-fit temperature is obtained as $T_{\rm s}=(6-12)\times
10^{4}~{}\mathrm{K}$, and we decided to use the proposed value. Note that the
conservative upper bound is also derived as $T_{\rm s}<1.7\times
10^{5}~{}\mathrm{K}$ by varying the parameter, such as the ratio of NS radius
and distance.

  * •

No. 17, PSR J2124-3358: Ref. [26] analyzed the optical data from the
J2124-3358. The BB+PL model gives the following possible range $T_{\rm
s}\in[0.5,2.1]\times 10^{5}~{}\mathrm{K}$ with the uncertainty from the
distance. We decided to select the original range, which is almost the same
uncertainty added by hand in the way we described in Sec. 5.

  * •

No. 19, RX J2143.0+0654: In Ref. [101], the X-ray data in XMM-Newton is fitted
using the BB absorption model, and the authors obtain the BB temperature
$k_{\rm B}T^{\infty}=104.0\pm 0.4~{}\mathrm{eV}$ with the BB radius as
$R^{\infty}=(3.10\pm 0.04)~{}\mathrm{km}$, where distance is fixed to be
$d=500~{}\mathrm{pc}$. This value is smaller than the typical NS radius, which
implies that this BB temperature is not from the surface but from the small
areas around the magnetic caps. In Ref. [102], the authors fitted the data
from Large Binocular Telescope (optical) by combining with the X-ray data in
XMM-Newton. Using the BB absorption model, $k_{\rm B}T^{\infty}=105.1\pm
0.9~{}\mathrm{eV}$ is obtained, which is consistent with the result in Ref.
[101]. They also perform fitting using the two-component BB model and the
hotter (cooler) component is obtained as $k_{\rm B}T=104~{}\mathrm{eV}$
($k_{\rm B}T=40~{}\mathrm{eV}$). It is impossible to eliminate uncertainty
from the model selection to fit the data from this situation, and thus we
choose all the possible ranges of the temperature, $k_{\rm B}T_{\rm
s}=40$–$106~{}\mathrm{eV}$.

  * •

No. 20, PSR J0108-1431: The X-ray data from the direction of J1080-1431
observed in XMM-Newton is fitted by the BB+PL model [159] with $k_{\rm
B}T=110^{+30}_{-10}~{}\mathrm{eV}$ and $R_{\rm
NS}=43^{+16}_{-9}~{}\mathrm{m}$. This small emission radius implies that this
BB component is not the cool surface temperature but the hot magnetic pole
component. We can interpret this result as the surface temperature is much
cooler than the magnetic pole, and thus, the hot component dominates the
observed flux. The latest analysis [28] analyses both the XMM-Newton and
optical data (HST, VLT). In particular, HST F140LP detected thermal emission,
and they put the conservative upper bound on the surface temperature as
$T_{\rm s}<5.9\times 10^{4}~{}\mathrm{K}$. To derive this conservative bound,
they included uncertainty from the parallax distance [160]. Furthermore, they
obtain the value of $T_{\rm s}$ by assuming the FUV flux is dominated by a
thermal component, $T_{\rm s}=27000-55000~{}\mathrm{K}$. We selected this
range to represent the uncertainty.

  * •

No. 22, PSR J2144-3933: The upper bound is obtained for the surface
temperature of J2144-3933 using XMM-Newton data (combining with the optical
data of Very Large Telescope (VLT)) [161] as $T_{\rm s}<2.3\times
10^{5}~{}\mathrm{K}$. The latest analysis [104] used deep optical and FUV
observation data by HST and derived the upper bound on the surface temperature
of J2144-3933. The conservative upper bound on the surface temperature is
derived based on the non-detection,

$\displaystyle T_{\rm s}<4.2\times 10^{4}~{}\mathrm{K},$ (B.38)

where a range of NS radius $R_{\rm NS}=[11,13]~{}\mathrm{km}$ and parallax
distance $d=172^{+20}_{-15}~{}\mathrm{pc}$ is considered to estimate the
uncertainty [160]. In this analysis, NS mass is fixed as $M_{\rm
NS}=1.4~{}M_{\odot}$.

Let us also comment on the rejected observational data from our list.

  * $-$

PSR B1929+10: The BB+PL fit is performed for the X-ray data [162]. However,
the magnetic pole component is reported to dominate the temperature because
the fitted radiation radius is much smaller than the NS radius. We conclude
this data is not appropriate to test the vortex creep heating.

  * $-$

XMMU J1732-344: We also omit XMMU J1732-344 from our list because the observed
value of $|\dot{\Omega}|$ is not determined. Once its pulsation data is fixed,
it is worth studying whether the vortex creep heating can explain this data;
its thermal emission is expected to exceed the value expected from minimal
cooling [163, 164] with its kinetic time information [163].

## References

  * [1] A. Bohr, B. R. Mottelson, and D. Pines, Possible Analogy between the Excitation Spectra of Nuclei and Those of the Superconducting Metallic State, Phys. Rev. 110 (1958) 936–938.
  * [2] A. B. Migdal, Superfluidity and the moments of inertia of nuclei, Nuclear Physics 13 (1959) 655–674.
  * [3] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, Stellar Superfluids, arXiv:1302.6626 (2013).
  * [4] B. Haskell and A. Sedrakian, Superfluidity and Superconductivity in Neutron Stars, Astrophys. Space Sci. Libr. 457 (2018) 401–454 [arXiv:1709.10340].
  * [5] A. Sedrakian and J. W. Clark, Superfluidity in nuclear systems and neutron stars, Eur. Phys. J. A 55 (2019) 167 [arXiv:1802.00017].
  * [6] B. Haskell and A. Melatos, Models of Pulsar Glitches, Int. J. Mod. Phys. D 24 (2015) 1530008 [arXiv:1502.07062].
  * [7] J. R. Fuentes, et al., The glitch activity of neutron stars, Astron. Astrophys. 608 (2017) A131 [arXiv:1710.00952].
  * [8] D. Antonopoulou, B. Haskell, and C. M. Espinoza, Pulsar glitches: observations and physical interpretation, Rept. Prog. Phys. 85 (2022) 126901.
  * [9] S. Zhou, E. Gügercinoğlu, J. Yuan, M. Ge, and C. Yu, Pulsar Glitches: A Review, Universe 8 (2022) 641 [arXiv:2211.13885].
  * [10] P. W. Anderson and N. Itoh, Pulsar glitches and restlessness as a hard superfluidity phenomenon, Nature 256 (1975) 25–27.
  * [11] M. A. Alpar, Pinning and Threading of Quantized Vortices in the Pulsar Crust Superfluid, Astrophys. J. 213 (1977) 527–530.
  * [12] M. A. Alpar, D. Pines, P. W. Anderson, and J. Shaham, Vortex creep and the internal temperature of neutron stars. I - General theory, Astrophys. J. 276 (1984) 325–334.
  * [13] M. A. Alpar, P. W. Anderson, D. Pines, and J. Shaham, Vortex creep and the internal temperature of neutron stars. II. VELA pulsar., Astrophys. J. 278 (1984) 791–805.
  * [14] M. A. Alpar, R. Nandkumar, and D. Pines, Vortex creep and the internal temperature of neutron stars : the Crabpulsar and PSR 0525+21., Astrophys. J. 288 (1985) 191–195.
  * [15] N. Shibazaki and F. K. Lamb, Neutron Star Evolution with Internal Heating, Astrophys. J. 346 (1989) 808.
  * [16] M. A. Alpar, K. S. Cheng, and D. Pines, Vortex Creep and the Internal Temperature of Neutron Stars: Linear and Nonlinear Response to a Glitch, Astrophys. J. 346 (1989) 823.
  * [17] K. A. van Riper, R. I. Epstein, and G. S. Miller, Soft X-Ray Pulses from Neutron Star Glitches, Astrophys. J. Lett. 381 (1991) L47.
  * [18] H. Umeda, N. Shibazaki, K. Nomoto, and S. Tsuruta, Thermal Evolution of Neutron Stars with Internal Frictional Heating, Astrophys. J. 408 (1993) 186.
  * [19] K. A. van Riper, B. Link, and R. I. Epstein, Frictional Heating and Neutron Star Thermal Evolution, Astrophys. J. 448 (1995) 294 [astro-ph/9404060].
  * [20] C. Schaab, A. Sedrakian, F. Weber, and M. K. Weigel, Impact of internal heating on the thermal evolution of neutron stars, Astron. Astrophys. 346 (1999) 465–480 [astro-ph/9904127].
  * [21] M. B. Larson and B. Link, Superfluid friction and late-time thermal evolution of neutron stars, Astrophys. J. 521 (1999) 271 [astro-ph/9810441].
  * [22] D. Gonzalez and A. Reisenegger, Internal Heating of Old Neutron Stars: Contrasting Different Mechanisms, Astron. Astrophys. 522 (2010) A16 [arXiv:1005.5699].
  * [23] O. Kargaltsev, G. G. Pavlov, and R. W. Romani, Ultraviolet emission from the millisecond pulsar j0437-4715, Astrophys. J. 602 (2004) 327–335 [astro-ph/0310854].
  * [24] R. P. Mignani, G. G. Pavlov, and O. Kargaltsev, A possible optical counterpart to the old nearby pulsar J0108-1431, Astron. Astrophys. 488 (2008) 1027 [arXiv:0805.2586].
  * [25] M. Durant, et al., The spectrum of the recycled PSR J0437-4715 and its white dwarf companion, Astrophys. J. 746 (2012) 6 [arXiv:1111.2346].
  * [26] B. Rangelov, et al., Hubble Space Telescope Detection of the Millisecond Pulsar J2124$-$3358 and its Far-ultraviolet Bow Shock Nebula, Astrophys. J. 835 (2017) 264 [arXiv:1701.00002].
  * [27] G. G. Pavlov, et al., Old but still warm: Far-UV detection of PSR B0950+08, Astrophys. J. 850 (2017) 79 [arXiv:1710.06448].
  * [28] V. Abramkin, Y. Shibanov, R. P. Mignani, and G. G. Pavlov, Hubble Space Telescope Observations of the Old Pulsar PSR J0108–1431, Astrophys. J. 911 (2021) 1 [arXiv:2103.00332].
  * [29] V. Abramkin, G. G. Pavlov, Y. Shibanov, and O. Kargaltsev, Thermal and Nonthermal Emission in the Optical-UV Spectrum of PSR B0950+08*, Astrophys. J. 924 (2022) 128 [arXiv:2111.08801].
  * [30] D. G. Yakovlev, K. P. Levenfish, and Y. A. Shibanov, Cooling neutron stars and superfluidity in their interiors, Phys. Usp. 42 (1999) 737–778 [astro-ph/9906456].
  * [31] D. G. Yakovlev, A. D. Kaminker, O. Y. Gnedin, and P. Haensel, Neutrino emission from neutron stars, Phys. Rept. 354 (2001) 1 [astro-ph/0012122].
  * [32] D. G. Yakovlev and C. J. Pethick, Neutron star cooling, Ann. Rev. Astron. Astrophys. 42 (2004) 169–210 [astro-ph/0402143].
  * [33] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, Minimal cooling of neutron stars: A New paradigm, Astrophys. J. Suppl. 155 (2004) 623–650 [astro-ph/0403657].
  * [34] M. E. Gusakov, A. D. Kaminker, D. G. Yakovlev, and O. Y. Gnedin, Enhanced cooling of neutron stars via Cooper-pairing neutrino emission, Astron. Astrophys. 423 (2004) 1063–1072 [astro-ph/0404002].
  * [35] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, Neutrino Emission from Cooper Pairs and Minimal Cooling of Neutron Stars, Astrophys. J. 707 (2009) 1131–1140 [arXiv:0906.1621].
  * [36] J. M. Lattimer, K. A. van Riper, M. Prakash, and M. Prakash, Rapid Cooling and the Structure of Neutron Stars, Astrophys. J. 425 (1994) 802.
  * [37] O. Y. Gnedin, D. G. Yakovlev, and A. Y. Potekhin, Thermal relaxation in young neutron stars, Mon. Not. Roy. Astron. Soc. 324 (2001) 725 [astro-ph/0012306].
  * [38] K. S. Thorne, The relativistic equations of stellar structure and evolution., Astrophys. J. 212 (1977) 825–831.
  * [39] G. Glen and P. Sutherland, On the cooling of neutron stars, Astrophys. J. 239 (1980) 671–684.
  * [40] A. Y. Potekhin, G. Chabrier, and D. G. Yakovlev, Internal temperatures and cooling of neutron stars with accreted envelopes, Astron. Astrophys. 323 (1997) 415 [astro-ph/9706148].
  * [41] F. Köpp, J. E. Horvath, D. Hadjimichef, C. A. Z. Vasconcellos, and P. O. Hess, Internal heating mechanisms in neutron stars, Int. J. Mod. Phys. D 32 (2023) 2350046 [arXiv:2208.07770].
  * [42] G. Baym, C. Pethick, D. Pines, and M. Ruderman, Spin Up in Neutron Stars : The Future of the Vela Pulsar, Nature, 224 (1969) 872–874.
  * [43] J. A. Sauls, D. L. Stein, and J. W. Serene, Magnetic vortices in a rotating ${}^{3}P_{2}$ neutron superfluid, Phys. Rev. D 25 (1982) 967–975.
  * [44] M. A. Alpar, S. A. Langer, and J. A. Sauls, Rapid postglitch spin-up of the superfluid core in pulsars., Astrophys. J. 282 (1984) 533–541.
  * [45] E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics: Theory of the Condensed State (Part 2), $\S~{}29$, (1980).
  * [46] G. A. Williams and R. E. Packard, Photographs of Quantized Vortex Lines in Rotating He II, Phys. Rev. Lett. 33 (1974) 280–283.
  * [47] D. V. Osborne, The Rotation of Liquid Helium II, Proceedings of the Physical Society. Section A 63 (1950) 909.
  * [48] A. J. Manninen, et al., First optical observations of superfluid ${}^{3}\mathrm{He}$, Phys. Rev. Lett. 69 (1992) 2392–2395.
  * [49] G. Baym, R. I. Epstein, and B. Link, Dynamics of vortices in neutron stars, Physica B Condensed Matter 178 (1992) 1–12.
  * [50] B. Link, R. I. Epstein, and G. Baym, Superfluid Vortex Creep and Rotational Dynamics of Neutron Stars, Astrophys. J. 403 (1993) 285.
  * [51] B. K. Link and R. I. Epstein, Mechanics and Energetics of Vortex Unpinning in Neutron Stars, Astrophys. J. 373 (1991) 592.
  * [52] T. Hamada and I. D. Johnston, A potential model representation of two-nucleon data below $315\,\mathrm{MeV}$, Nucl. Phys. 34 (1962) 382.
  * [53] R. Machleidt, The Meson Theory of Nuclear Forces and Nuclear Structure, in Advances in Nuclear Physics, Volume 19, J. W. Negele, et al., eds., ch. 2, p. 189.  Springer US, Boston, 1989.
  * [54] R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, An Accurate nucleon-nucleon potential with charge independence breaking, Phys. Rev. C 51 (1995) 38–51 [nucl-th/9408016].
  * [55] T. H. R. Skyrme, CVII. The nuclear surface, Phil. Mag. 1 (1956) 1043–1054.
  * [56] T. H. R. Skyrme, The effective nuclear potential, Nucl. Phys. 9 (1958) 615.
  * [57] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, A Skyrme parametrization from subnuclear to neutron star densities, Nuclear Physics A 627 (1997) 710–746.
  * [58] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Self-consistent mean-field models for nuclear structure, Rev. Mod. Phys. 75 (2003) 121–180.
  * [59] M. Dutra, et al., Skyrme interaction and nuclear matter constraints, Phys. Rev. C 85 (2012) 035201.
  * [60] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, A Skyrme parametrization from subnuclear to neutron star densities. 2. Nuclei far from stablities, Nucl. Phys. A 635 (1998) 231–256. [Erratum: Nucl.Phys.A 643, 441–441 (1998)].
  * [61] J. Bartel, P. Quentin, M. Brack, C. Guet, and H. B. Håkansson, Towards a better parametrisation of Skyrme-like effective forces: A critical study of the SkM force, Nuclear Physics A 386 (1982) 79–100.
  * [62] J. Dechargé and D. Gogny, Hartree-Fock-Bogolyubov calculations with the $D1$ effective interaction on spherical nuclei, Phys. Rev. C 21 (1980) 1568–1593.
  * [63] P. Avogadro, F. Barranco, R. A. Broglia, and E. Vigezzi, Quantum calculation of vortices in the inner crust of neutron stars, Phys. Rev. C 75 (2007) 012805 [nucl-th/0602028].
  * [64] P. Avogadro, F. Barranco, R. A. Broglia, and E. Vigezzi, Microscopic quantal calculation of the vortex-nucleus interaction in the inner crust of neutron stars, Nucl. Phys. A 788 (2007) 130–135.
  * [65] P. Avogadro, F. Barranco, R. A. Broglia, and E. Vigezzi, Vortex-nucleus interaction in the inner crust of neutron stars, Nucl. Phys. A 811 (2008) 378–412 [arXiv:0804.1765].
  * [66] P. Klausner, F. Barranco, P. M. Pizzochero, X. Roca-Maza, and E. Vigezzi, Microscopic calculation of the pinning energy of a vortex in the inner crust of a neutron star, arXiv:2303.18151 (2023).
  * [67] R. I. Epstein and G. Baym, Vortex Pinning in Neutron Stars, Astrophys. J. 328 (1988) 680.
  * [68] P. M. Pizzochero, L. Viverit, and R. A. Broglia, Vortex nucleus interaction and pinning forces in neutron stars, Phys. Rev. Lett. 79 (1997) 3347–3350 [astro-ph/9709060].
  * [69] P. Donati and P. M. Pizzochero, Is there Nuclear Pinning of Vortices in Superfluid Pulsars?, Phys. Rev. Lett. 90 (2003) 211101.
  * [70] P. Donati and P. M. Pizzochero, Realistic energies for vortex pinning in intermediate-density neutron star matter, Phys. Lett. B 640 (2006) 74–81.
  * [71] P. Donati and P. M. Pizzochero, Fully consistent semi-classical treatment of vortex–nucleus interaction in rotating neutron stars, Nucl. Phys. A 742 (2004) 363–379.
  * [72] S. Seveso, P. M. Pizzochero, F. Grill, and B. Haskell, Mesoscopic pinning forces in neutron star crusts, Monthly Notices of the Royal Astronomical Society 455 (2015) 3952–3967.
  * [73] E. Garrido, P. Sarriguren, E. Moya de Guerra, and P. Schuck, Effective density dependent pairing forces in the T = 1 and T = 0 channels, Phys. Rev. C 60 (1999) 064312 [nucl-th/9909026].
  * [74] B. Link, Dynamics of Quantum Vorticity in a Random Potential, Phys. Rev. Lett. 102 (2009) 131101 [arXiv:0807.1945].
  * [75] A. Bulgac, M. M. Forbes, and R. Sharma, Strength of the Vortex-Pinning Interaction from Real-Time Dynamics, Phys. Rev. Lett. 110 (2013) 241102 [arXiv:1302.2172].
  * [76] G. Wlazłowski, K. Sekizawa, P. Magierski, A. Bulgac, and M. M. Forbes, Vortex pinning and dynamics in the neutron star crust, Phys. Rev. Lett. 117 (2016) 232701 [arXiv:1606.04847].
  * [77] B. Link and Y. Levin, Vortex Pinning in Neutron Stars, Slipstick Dynamics, and the Origin of Spin Glitches, Astrophys. J. 941 (2022) 148 [arXiv:2208.11494].
  * [78] M. Antonelli and B. Haskell, Superfluid vortex-mediated mutual friction in non-homogeneous neutron star interiors, Mon. Not. Roy. Astron. Soc. 499 (2020) 3690–3705 [arXiv:2007.11720].
  * [79] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, The Equation of state of nucleon matter and neutron star structure, Phys. Rev. C 58 (1998) 1804–1828 [nucl-th/9804027].
  * [80] NSCool, http://www.astroscu.unam.mx/neutrones/NSCool/.
  * [81] J. W. Negele and D. Vautherin, Neutron star matter at subnuclear densities, Nucl. Phys. A 207 (1973) 298–320.
  * [82] P. Haensel, J. L. Zdunik, and J. Dobaczewski, Composition and equation of state of cold catalyzed matter below neutron drip, Astron. Astrophys. 222 (1989) 353–357.
  * [83] K. E. McGowan, et al., XMM - Newton observations of PSR B1706-44, Astrophys. J. 600 (2004) 343–350 [astro-ph/0309445].
  * [84] O. Kargaltsev, M. Durant, Z. Misanovic, and G. Pavlov, Absorption Features in the X-ray Spectrum of an Ordinary Radio Pulsar, Science 337 (2012) 946–949 [arXiv:1208.5400].
  * [85] K. E. McGowan, et al., Evidence for surface cooling emission in the xmm-newton spectrum of the x-ray pulsar psr b2334+61, Astrophys. J. 639 (2006) 377–381 [astro-ph/0508439].
  * [86] P. Arumugasamy, O. Kargaltsev, B. Posselt, G. G. Pavlov, and J. Hare, Possible phase-dependent absorption feature in the x-ray spectrum of the middle-aged PSR J0659+1414, Astrophys. J. 869 (2018) 97 [arXiv:1810.11814].
  * [87] K. Mori et al., A broadband X-ray study of the Geminga pulsar with NuSTAR and XMM-Newton, Astrophys. J. 793 (2014) 88 [arXiv:1408.1644].
  * [88] C.-Y. Ng, R. W. Romani, W. F. Brisken, S. Chatterjee, and M. Kramer, The Origin and Motion of PSR J0538+2817 in S147, Astrophys. J. 654 (2006) 487–493 [astro-ph/0611068].
  * [89] A. De Luca, P. A. Caraveo, S. Mereghetti, M. Negroni, and G. F. Bignami, On the polar caps of the 3 Musketeers, Astrophys. J. 623 (2005) 1051 [astro-ph/0412662].
  * [90] N. Tetzlaff, J. G. Schmidt, M. M. Hohle, and R. Neuhaeuser, Neutron stars from young nearby associations the origin of RXJ1605.3+3249, Publ. Astron. Soc. Austral. 29 (2012) 98 [arXiv:1202.1388].
  * [91] C. Malacaria et al., A joint NICER and XMM-Newton view of the ”Magnificent” thermally emitting X-ray Isolated Neutron Star RX J1605.3+3249, Astrophys. J. 880 (2019) 74 [arXiv:1906.02806].
  * [92] S. Beloin, S. Han, A. W. Steiner, and D. Page, Constraining Superfluidity in Dense Matter from the Cooling of Isolated Neutron Stars, Phys. Rev. C 97 (2018) 015804 [arXiv:1612.04289].
  * [93] K. Auchettl et al., X-ray analysis of the proper motion and pulsar wind nebula for PSR J1741-2054, Astrophys. J. 802 (2015) 68 [arXiv:1501.03225].
  * [94] A. Kirichenko, et al., Deep optical observations of the $\gamma$-ray pulsar J0357+3205, Astron. Astrophys. 564 (2014) A81 [arXiv:1402.2246].
  * [95] D. L. Kaplan and M. H. van Kerkwijk, A Coherent Timing Solution for the Nearby, Thermally Emitting Isolated Neutron Star RX J0420.0-5022, Astrophys. J. Lett. 740 (2011) L30 [arXiv:1109.2105].
  * [96] C. Motch, A. M. Pires, F. Haberl, A. Schwope, and V. E. Zavlin, Proper motions of thermally emitting isolated neutron stars measured with Chandra, Astron. Astrophys. 497 (2009) 423 [arXiv:0901.1006].
  * [97] A. D. Schwope, V. Hambaryan, F. Haberl, and C. Motch, The complex X-ray spectrum of the isolated neutron star RBS1223, Astrophys. Space Sci. 308 (2007) 619–623 [astro-ph/0609705].
  * [98] N. Tetzlaff, T. Eisenbeiss, R. Neuhaeuser, and M. M. Hohle, The origin of RXJ1856.5-3754 and RXJ0720.4-3125 – updated using new parallax measurements, Mon. Not. Roy. Astron. Soc. 417 (2011) 617 [arXiv:1107.1673].
  * [99] V. Hambaryan, et al., The compactness of the isolated neutron star RX J0720.4$-$3125, Astron. Astrophys. 601 (2017) A108 [arXiv:1702.07635].
  * [100] N. Sartore, et al., Spectral monitoring of RX J1856.5-3754 with XMM-Newton. Analysis of EPIC-pn data, Astron. Astrophys. 541 (2012) A66 [arXiv:1202.2121].
  * [101] D. L. Kaplan and M. H. van Kerkwijk, Constraining the Spin-down of the Nearby Isolated Neutron Star RX J2143.0+0654, Astrophys. J. Lett. 692 (2009) L62–L66 [arXiv:0901.4133].
  * [102] A. Schwope et al., The isolated neutron star RBS1774 revisited. Revised XMM-Newton X-ray parameters and an optical counterpart from deep LBT-observations, Astron. Astrophys. 499 (2009) 267 [arXiv:0902.4110].
  * [103] D. L. Kaplan and M. H. van Kerkwijk, Constraining the Spin-down of the Nearby Isolated Neutron Star RX J0806.4-4123, and Implications for the Population of Nearby Neutron Stars, Astrophys. J. 705 (2009) 798–808 [arXiv:0909.5218].
  * [104] S. Guillot, et al., Hubble Space Telescope Nondetection of PSR J2144–3933: The Coldest Known Neutron Star, Astrophys. J. 874 (2019) 175 [arXiv:1901.07998].
  * [105] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, The Australia Telescope National Facility pulsar catalogue, Astron. J. 129 (2005) 1993 [astro-ph/0412641].
  * [106] ATNF Pulsar Catalogue, https://www.atnf.csiro.au/research/pulsar/psrcat.
  * [107] J. A. Pons, J. A. Miralles, and U. Geppert, Magneto–thermal evolution of neutron stars, Astron. Astrophys. 496 (2009) 207–216 [arXiv:0812.3018].
  * [108] S. Dall’Osso, J. Granot, and T. Piran, Magnetic field decay in neutron stars: from soft gamma repeaters to ’weak-field magnetars’, Mon. Not. Roy. Astron. Soc. 422 (2012) 2878–2903 [arXiv:1110.2498].
  * [109] D. Viganò, et al., Unifying the observational diversity of isolated neutron stars via magneto-thermal evolution models, Mon. Not. Roy. Astron. Soc. 434 (2013) 123 [arXiv:1306.2156].
  * [110] A. Y. Potekhin, D. A. Zyuzin, D. G. Yakovlev, M. V. Beznogov, and Y. A. Shibanov, Thermal luminosities of cooling neutron stars, Mon. Not. Roy. Astron. Soc. 496 (2020) 5052–5071 [arXiv:2006.15004].
  * [111] A. Reisenegger, Deviations from chemical equilibrium due to spindown as an internal heat source in neutron stars, Astrophys. J. 442 (1995) 749 [astro-ph/9410035].
  * [112] P. Haensel, Non-equilibrium neutrino emissivities and opacities of neutron star matter, Astron. Astrophys. 262 (1992) 131–137.
  * [113] E. Gourgoulhon and P. Haensel, Upper bounds on the neutrino burst from collapse of a neutron star into a black hole, Astron. Astrophys. 271 (1993) 187.
  * [114] R. Fernandez and A. Reisenegger, Rotochemical heating in millisecond pulsars. Formalism and non-superfluid case, Astrophys. J. 625 (2005) 291–306 [astro-ph/0502116].
  * [115] L. Villain and P. Haensel, Non-equilibrium beta processes in superfluid neutron star cores, Astron. Astrophys. 444 (2005) 539 [astro-ph/0504572].
  * [116] C. Petrovich and A. Reisenegger, Rotochemical heating in millisecond pulsars: modified Urca reactions with uniform Cooper pairing gaps, Astron. Astrophys. 521 (2010) A77 [arXiv:0912.2564].
  * [117] C.-M. Pi, X.-P. Zheng, and S.-H. Yang, Neutrino Emissivity of Non-equilibrium beta processes With Nucleon Superfluidity, Phys. Rev. C 81 (2010) 045802 [arXiv:0912.2884].
  * [118] N. González-Jiménez, C. Petrovich, and A. Reisenegger, Rotochemical heating of millisecond and classical pulsars with anisotropic and density-dependent superfluid gap models, Mon. Not. Roy. Astron. Soc. 447 (2015) 2073 [arXiv:1411.6500].
  * [119] K. Yanagi, N. Nagata, and K. Hamaguchi, Cooling Theory Faced with Old Warm Neutron Stars: Role of Non-Equilibrium Processes with Proton and Neutron Gaps, Mon. Not. Roy. Astron. Soc. 492 (2020) 5508–5523 [arXiv:1904.04667].
  * [120] D. Toyouchi, K. Hotokezaka, and M. Takada, Hunting wide-area optical surveys for high proper motion isolated neutron stars, Mon. Not. Roy. Astron. Soc. 510 (2021) 611–619 [arXiv:2106.04846].
  * [121] A. Schwenk, B. Friman, and G. E. Brown, Renormalization group approach to neutron matter: Quasiparticle interactions, superfluid gaps and the equation of state, Nucl. Phys. A 713 (2003) 191–216 [nucl-th/0207004].
  * [122] J. M. C. Chen, J. W. Clark, R. D. Davé, and V. V. Khodel, Pairing gaps in nucleonic superfluids, Nucl. Phys. A 555 (1993) 59–89.
  * [123] C. Kouvaris, WIMP Annihilation and Cooling of Neutron Stars, Phys. Rev. D 77 (2008) 023006 [arXiv:0708.2362].
  * [124] G. Bertone and M. Fairbairn, Compact Stars as Dark Matter Probes, Phys. Rev. D 77 (2008) 043515 [arXiv:0709.1485].
  * [125] C. Kouvaris and P. Tinyakov, Can Neutron stars constrain Dark Matter?, Phys. Rev. D 82 (2010) 063531 [arXiv:1004.0586].
  * [126] A. de Lavallaz and M. Fairbairn, Neutron Stars as Dark Matter Probes, Phys. Rev. D 81 (2010) 123521 [arXiv:1004.0629].
  * [127] J. Bramante, A. Delgado, and A. Martin, Multiscatter stellar capture of dark matter, Phys. Rev. D 96 (2017) 063002 [arXiv:1703.04043].
  * [128] M. Baryakhtar, J. Bramante, S. W. Li, T. Linden, and N. Raj, Dark Kinetic Heating of Neutron Stars and An Infrared Window On WIMPs, SIMPs, and Pure Higgsinos, Phys. Rev. Lett. 119 (2017) 131801 [arXiv:1704.01577].
  * [129] N. Raj, P. Tanedo, and H.-B. Yu, Neutron stars at the dark matter direct detection frontier, Phys. Rev. D 97 (2018) 043006 [arXiv:1707.09442].
  * [130] C.-S. Chen and Y.-H. Lin, Reheating neutron stars with the annihilation of self-interacting dark matter, JHEP 08 (2018) 069 [arXiv:1804.03409].
  * [131] N. F. Bell, G. Busoni, and S. Robles, Heating up Neutron Stars with Inelastic Dark Matter, JCAP 09 (2018) 018 [arXiv:1807.02840].
  * [132] R. Garani, Y. Genolini, and T. Hambye, New Analysis of Neutron Star Constraints on Asymmetric Dark Matter, JCAP 05 (2019) 035 [arXiv:1812.08773].
  * [133] D. A. Camargo, F. S. Queiroz, and R. Sturani, Detecting Dark Matter with Neutron Star Spectroscopy, JCAP 09 (2019) 051 [arXiv:1901.05474].
  * [134] N. F. Bell, G. Busoni, and S. Robles, Capture of Leptophilic Dark Matter in Neutron Stars, JCAP 06 (2019) 054 [arXiv:1904.09803].
  * [135] K. Hamaguchi, N. Nagata, and K. Yanagi, Dark Matter Heating vs. Rotochemical Heating in Old Neutron Stars, Phys. Lett. B 795 (2019) 484–489 [arXiv:1905.02991].
  * [136] R. Garani and J. Heeck, Dark matter interactions with muons in neutron stars, Phys. Rev. D 100 (2019) 035039 [arXiv:1906.10145].
  * [137] J. F. Acevedo, J. Bramante, R. K. Leane, and N. Raj, Warming Nuclear Pasta with Dark Matter: Kinetic and Annihilation Heating of Neutron Star Crusts, JCAP 03 (2020) 038 [arXiv:1911.06334].
  * [138] A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Relativistic capture of dark matter by electrons in neutron stars, Phys. Lett. B (2020) 135767 [arXiv:1911.13293].
  * [139] W.-Y. Keung, D. Marfatia, and P.-Y. Tseng, Heating neutron stars with GeV dark matter, JHEP 07 (2020) 181 [arXiv:2001.09140].
  * [140] K. Yanagi, Thermal Evolution of Neutron Stars as a Probe of Physics beyond the Standard Model, arXiv:2003.08199 (2020).
  * [141] A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Dark kinetic heating of neutron stars from contact interactions with relativistic targets, Phys. Rev. D 102 (2020) 123002 [arXiv:2004.09539].
  * [142] N. F. Bell, G. Busoni, S. Robles, and M. Virgato, Improved Treatment of Dark Matter Capture in Neutron Stars, JCAP 09 (2020) 028 [arXiv:2004.14888].
  * [143] N. F. Bell, G. Busoni, S. Robles, and M. Virgato, Improved Treatment of Dark Matter Capture in Neutron Stars II: Leptonic Targets, JCAP 03 (2021) 086 [arXiv:2010.13257].
  * [144] F. Anzuini, et al., Improved treatment of dark matter capture in neutron stars III: nucleon and exotic targets, JCAP 11 (2021) 056 [arXiv:2108.02525].
  * [145] Y.-P. Zeng, X. Xiao, and W. Wang, Constraints on Pseudo-Nambu-Goldstone dark matter from direct detection experiment and neutron star reheating temperature, Phys. Lett. B 824 (2022) 136822 [arXiv:2108.11381].
  * [146] J. Bramante, B. J. Kavanagh, and N. Raj, Scattering Searches for Dark Matter in Subhalos: Neutron Stars, Cosmic Rays, and Old Rocks, Phys. Rev. Lett. 128 (2022) 231801 [arXiv:2109.04582].
  * [147] P. Tinyakov, M. Pshirkov, and S. Popov, Astroparticle Physics with Compact Objects, Universe 7 (2021) 401 [arXiv:2110.12298].
  * [148] T. N. Maity and F. S. Queiroz, Detecting bosonic dark matter with neutron stars, Phys. Rev. D 104 (2021) 083019 [arXiv:2104.02700].
  * [149] M. Fujiwara, K. Hamaguchi, N. Nagata, and J. Zheng, Capture of electroweak multiplet dark matter in neutron stars, Phys. Rev. D 106 (2022) 055031 [arXiv:2204.02238].
  * [150] K. Hamaguchi, N. Nagata, and M. E. Ramirez-Quezada, Neutron star heating in dark matter models for the muon g $-$ 2 discrepancy, JHEP 10 (2022) 088 [arXiv:2204.02413].
  * [151] S. Chatterjee, et al., Faint light of old neutron stars and detectability at the James Webb Space Telescope, Phys. Rev. D 108 (2023) L021301 [arXiv:2205.05048].
  * [152] J. Coffey, D. McKeen, D. E. Morrissey, and N. Raj, Neutron star observations of pseudoscalar-mediated dark matter, Phys. Rev. D 106 (2022) 115019 [arXiv:2207.02221].
  * [153] J. T. Acuña, P. Stengel, and P. Ullio, Neutron Star Heating in Dark Matter Models for Muon g-2 with Scalar Lepton Partners up to the TeV Scale, arXiv:2209.12552 (2022).
  * [154] G. Alvarez, A. Joglekar, M. Phoroutan-Mehr, and H.-B. Yu, Heating neutron stars with inelastic dark matter and relativistic targets, Phys. Rev. D 107 (2023) 103024 [arXiv:2301.08767].
  * [155] J. Bramante and N. Raj, Dark matter in compact stars, arXiv:2307.14435 (2023).
  * [156] M. Fujiwara, K. Hamaguchi, N. Nagata, and E. M. Ramirez-Quezada, to be published,.
  * [157] W. Becker, et al., Revealing the x-ray emission processes of old rotation-powered pulsars: XMM-Newton observations of PSR B0950+08, PSR B0823+26 and PSR J2043+2740, Astrophys. J. 615 (2004) 908–920 [astro-ph/0405180].
  * [158] V. E. Zavlin and G. G. Pavlov, X-ray emission from the old pulsar B0950+08, Astrophys. J. 616 (2004) 452–462 [astro-ph/0405212].
  * [159] B. Posselt, et al., XMM-Newton observation of the very old pulsar J0108-1431, Astrophys. J. 761 (2012) 117 [arXiv:1210.7179].
  * [160] J. P. W. Verbiest, D. R. Lorimer, and M. A. McLaughlin, Lutz-Kelker bias in pulsar parallax measurements, Mon. Not. Roy. Astron. Soc. 405 (2010) 564–572 [arXiv:1002.1213].
  * [161] A. Tiengo, et al., The first deep X-ray and optical observations of the closest isolated radio pulsar, AIP Conf. Proc. 1357 (2011) 269–272 [arXiv:1103.4538].
  * [162] Z. Misanovic, G. Pavlov, and G. Garmire, Chandra observations of the pulsar PSR B1929+10 and its environment, Astrophys. J. 685 (2008) 1129–1142 [arXiv:0711.4171].
  * [163] W. W. Tian, D. A. Leahy, M. Haverkorn, and B. Jiang, Discovery of the counterpart of TeV Gamma-ray source HESS J1731-347: a new SNR G353.6-0.7 with radio and X-ray images, Astrophys. J. Lett. 679 (2008) L85 [arXiv:0801.3254].
  * [164] D. Klochkov, et al., The neutron star in HESSJ1731-347: Central compact objects as laboratories to study the equation of state of superdense matter, Astron. Astrophys. 573 (2015) A53 [arXiv:1410.1055].