Vortex creep heating in neutron stars

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 proportionality 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.


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, 1 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.

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 redshifted temperature of NSs, T ∞ (r, t) = T (r, t)e ϕ (r) , where ϕ(r) specifies the gravitational redshift, reaches a constant value T ∞ (r, t) ≃ T ∞ (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] 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 ∞ ν , photon cooling L ∞ γ , and heating source L ∞ H .At temperatures below a few × 10 9 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 ≳ 10 5 yrs,2 photon emission dominates neutrino emission.The thermal photon emission follows a blackbody spectrum.The photon luminosity L γ can be described by the Stefan-Boltzmann law and related to surface temperature T s as in the local reference frame of the NS, where σ SB is the Stephan-Boltzmann constant and R NS the NS radius.To relate the internal temperature T of a NS to its surface temperature T 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, 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 yrs.This can be seen by estimating the typical timescale τ eq for the relaxation into the equilibrium state: where we have used T ∞ ∼ 10 6 K corresponding to the surface temperature, at the equilibrium phase, of T s ∼ 10 5 K. 3 It is found that this timescale is shorter than the age of old NSs, assuring the equilibrium condition (2.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.

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. 4This 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.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 ∼ 1 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 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, φ, z).Since we treat the crust component as a rigid body, its angular velocity Ω c (t) is independent of r.On the other hand, the superfluid angular velocity varies with r, and we denote it by Ω s (t, r).The equation of motion for the crust component is where I 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 ext (t) represents the external torque acting on the system.The second term N int (t) corresponds to the effects of internal torques in the system 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 ω characterizes the microscopic features of rotation and is a locally determined value defined as the curl of the fluid velocity v, ω ≡ ∇ × v. (3.8) In contrast, the circulation Γ measures the macroscopic rotation and is defined over a finite region.It is given by the line integral of the fluid velocity v around a closed path C, which can be expressed as a surface integral over any surface S with the boundary C, where dℓ and dS denote the line and surface elements, respectively.We used Stokes' theorem to obtain the last expression.Superfluid motion obeys the potential flow condition, where v 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 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π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, 5⟨v s ⟩ = Ω s × r , (3.12) where 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 Ω s .By substituting Eq. (3.12) in Eq. (3.9), we obtain 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, where v r is the vortex velocity in the radial direction 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: 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 ∂Ω s /∂t = 0.In this case, N int turns out to be zero, and the two equations of motion are decoupled.

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: where v 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.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, where E NP (E 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

Interstitial pinning Nuclear pinning
where ∆r is the distance between the nuclear and interstitial pinning positions and ∆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 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, where we use Eqs.(3.12) and (3.16) to obtain the last expression and introduce the relative angular velocity, This velocity difference induces the Magnus force per unit length of a vortex line, where ρ is the superfluid density and κ 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 |κ| ≡ κ.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 where Λ characterizes the vortex tension, T v = ρκ 2 Λ/(4π), and is evaluated as 2 ≲ Λ ≲ 10 in the inner crust [51].The transition rate is then given by [49,50] where T eff is defined by with T being the temperature of the inner crust and As can be seen from these expressions, the unpinning occurs predominantly through thermal activation (quantum tunneling) for T ≫ T q (T ≪ T q ).In particular, for NSs as old as those considered in the following analysis, their internal temperature is ≲ 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 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: where we set |E pin | = 1 MeV as a representative value.We find that this distance is considerably larger than the crust thickness (≃ 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., Note that Ω c and | Ω∞ | 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 pin = f Mag [50]: with

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] where we use Eqs.(3.6) and (3.7), and the inertial momenta I 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), where we define, and δΩ ∞ 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 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 s can be estimated using Eq.(2.3) as (3.33) (3.34) In the steady vortex creep scenario, therefore, the surface temperature is predicted as a function of J, | Ω∞ |, and R NS , free from the uncertainty of the initial condition and the subsequent temperature evolution.We will examine this prediction against observation in Sec. 5.

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 pin available in the literature in Sec.4.1; the values of f 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.

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 momentumdependent 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 manybody 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. 7This 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 O(R 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 (≃ 1/R 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 ∼ (10 2 -10 3 ) × R 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.
All in all, we have 2 × 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 pin /R WS for this configuration in Ref. [71].To convert this into the pinning force per unit length f pin , we multiply this by a factor of 1/(2R WS ), as cubic cells with the side length 2R WS are used in the calculation of Ref. [71].As a result, we obtain the values of f pin = (1-7) × 10 −3 MeV • fm −2 , depending on the position in the inner crust.We list quantities relevant to this calculation in Table 3 in Appendix A. 8 Notice that the values of the pinning force shown here should be regarded as a ballpark estimate.The lower limit (f pin = 1 × 10 −3 MeV • 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 ρ ∼ 10 13 g • cm −3 around which E pin ≃ 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 densitydependent 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

Mesoscopic
Mean field: Pairing: Hartree-Fock: Pairing: Woods-Saxon Argonne SLy4 contact force 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 MeV • fm −2 .
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). 9 The values of the pinning force estimated from the pinning energy and the Wigner-Seitz radius as in the semi-classical calculation are f pin = (3-4) × 10 −4 MeV • 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 pin = 8 × 10 −7 -4 × 10 −4 MeV • 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 10 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]. 11Such 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.

Theoretical evaluation of J
For the evaluation of J in Eq. (3.32), it is convenient to change the coordinate from cylindrical coordinates (r, φ, z) to spherical coordinates (R, θ, ϕ): where we approximate δΩ ∞ in Eq. (3.32) by δΩ 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 in , R 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 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 × 10 40 − 1.9 × 10 43 erg • s for the semi-classic mesoscopic calculation and J = 1.7 × 10 40 − 2.7 × 10 42 erg • s for the quantum mesoscopic calculation.

Vortex creep heating vs. observation
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:12

.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 ∼ 10 5 years.Since NSs are comparable in size and mass, we expect that J obs is also roughly equal (up to a factor of O(1)) for every NS.We test this expectation by using the data of J obs for old NSs.We also compare the values of J obs with the theoretical computations given in Sec. 4 [105,106].The value of J obs is evaluated as in Eq. (5.36) with R NS = 11.43 km.
In Table 2, we list the values of J obs for the NSs we consider in this paper.We select isolated NSs older than 10 4 yrs.In evaluating J obs , we have just assumed R NS = 11.43 km for all NSs, as the radius is poorly known for most NSs; we keep in mind that this may introduce an O(1) error in the determination of J obs .We also show the age, surface temperature, and Ω = 2π Ṗ /P 2 of the NSs, (P and Ṗ are the period and its time derivative, respectively).Regarding the NS age, we use the kinetic age t kin if available.Otherwise, we use the spin-down age t sd = P/(2 Ṗ ).We calculate t sd and Ω from P and Ṗ 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 vortexcreep heating.As we have discussed in Sec. 2, photon emission becomes the dominant cooling source for NSs older than ∼ 10 5 years.For younger NSs, this may not be the case, so 3), for which the values of J 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 [95,103] that experienced the decay of magnetic fields before.This may make these NSs hotter than ordinary NSs of the same age [95,[107][108][109], resulting in an overestimate of J obs .We denote these NSs by the type The uncertainty in the determination of J obs stems mainly from that in the surface temperature, which is significant due to its quartic dependence on T 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 O(1) uncertainty, as can be seen in, e.g., Ref. [110].Motivated by this, we include a factor of two uncertainty in T 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].
In Fig. 5, we show the range of J 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 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 obs for these NSs are in the same ballpark, J ∼ 10 43 erg • s, even though their | Ω|'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 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 O( 1)    2, with the same colors as in Fig. 5.
coming from the integration in Eq. (4.35), in addition to that from the estimation of f pin .The quantum mesoscopic calculation can explain some of the points with a small J 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. 13Nevertheless, obtaining a current compilation of the value of J suggested by the observation is intriguing.Considering intrinsic 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 J ≃ 10 42.9−43.8erg • s , ( which we show as the red band 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). 14The 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. 6a, we take P Ṗ = 10 −15 s and the initial period P 0 = 10 ms to calculate | Ω(t)|, where we assume that the external torque is dominated by magnetic dipole radiation. 15These values are typical for ordinary pulsars.Nevertheless, we note that | Ω(t)| obtained with these parameters do not exactly agree with the observed values of | Ω| in Table 2, so the data points shown in this figure should be regarded as just an eye guide.In Fig. 6b, we set P Ṗ = 3.3 × 10 −22 s, which is the observed value for PSR J2124-3358, and P 0 = 1 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 ∼ 10 5 years and remains high enough at later times to be compatible with the observed data.

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 | Ω|, 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 ∼ 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 .In this case, the surface temperature at late times is predicted to be a few × 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 15 In this case, we have Ω ∝ −Ω 3 , i.e., P Ṗ = constant, and by solving this we obtain with t sd,0 ≡ P 2 0 /(2P Ṗ ).For the choice of parameters in Fig. 6a and Fig. 6b, t sd,0 ≃ 2 × 10 3 and 5 × 10 7 years, respectively.The value of P Ṗ is related to the surface magnetic flux density B s .In the ATNF pulsar catalogue [105,106], B s = 3.2 × 10 19 (P Ṗ ) 1/2 G is used for this relation, with which we have B s ≃ 1.0 × 10 12 G and 5.8 × 10 8 G in Fig. 6a and Fig. 6b  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.
only for the nuclear pinning, just for easy comparison with the semi-classical calculations.
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 β = 1 introduced in Ref. [72].Because of the averaging procedure, we find that a larger L results in a smaller value of f pin .
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., β = 1 as in the previous case.
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 pin are distributed in the range 10 −8 -10 −2 MeV • fm −2 , depending on the evaluation scheme and the selected nuclear potential.

B Selection criteria of NS data
We explain how we choose the range of uncertainty of T 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 atmo-sphere+PL models, and only BB+PL and atmosphere+PL models result in acceptable χ 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 ∞ = (0.48 − 2.2) × 10 6 K, to include the uncertainty coming from the choice of fitting models.
model zone  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.The filled and opened markers correspond to the nuclear and interstitial pinnings.
• 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 ∞ s < 6.27 × 10 5 K, where R NS = 10 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 ∞ s ≃ 9×10 5 K with the radiation radius R ∞ ≃ 2 km [158].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 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 s = (6 − 12) × 10 4 K, and we decided to use the proposed value.Note that the conservative upper bound is also derived as T s < 1.7 × 10 5 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 s ∈ [0.5, 2.1] × 10 5 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 B T ∞ = 104.0± 0.4 eV with the BB radius as R ∞ = (3.10 ± 0.04) km, where distance is fixed to be d = 500 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 B T ∞ = 105.1 ± 0.9 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 B T = 104 eV (k B T = 40 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 B T s = 40-106 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 B T = 110 +30 −10 eV and R NS = 43 +16 −9 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 s < 5.9 × 10 4 K.To derive this conservative bound, they included uncertainty from the parallax distance [160].Furthermore, they obtain the value of T s by assuming the FUV flux is dominated by a thermal component, T s = 27000 − 55000 K.We selected this range to represent the uncertainty.where a range of NS radius R NS = [11,13] km and parallax distance d = 172 +20 −15 pc is considered to estimate the uncertainty [160].In this analysis, NS mass is fixed as 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 | Ω| 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].

Figure 1 :
Fig. 1 in paper1Figure1: 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.

Figure 2 :
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.

Figure 3 :
Figure3: The Magnus force acting on a vortex line.We define the direction of κ as the direction of the right-hand screw.

Figure 4 :
Figure 4: The landscape of a vortex-line configuration at different length-scales.
[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].

Figure 5 :PP = 3 . 3 ×
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 ≃ 10 42.9-43.8erg • s.For comparison, we also show the values of J estimated with the mesoscopic calculations by black bars.

Figure 6 :
Figure6: 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 Table2, with the same colors as in Fig.5.

1 ρ 2 ]n [fm - 3 ]Figure 7 :
Figure 7: The values of f pin given in the tables in this appendix against the density ρ.The filled and opened markers correspond to the nuclear and interstitial pinnings.

•
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 s < 2.3 × 10 5 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, T s < 4.2 × 10 4 K,(B.38) .2.

Table 2 :
The data of the NSs considered in this paper.We classify them into four types-ordinary pulsars younger than 10 5 yrs [Y], ordinary pulsars older than 10 5 yrs [O], XDINSs [X], and millisecond pulsars [M].The values of t sd and Ω without references are computed from the data given in the ATNF pulsar catalogue

Table 3 :
[156]pectively.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.ρ,RWS , and ξ are the mass density, Wigner-Seitz radius, and coherence length, respectively.oldNSs.A detailed study of this issue will be given in the forthcoming paper[156].

Table 4 :
[71]tities 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.