Rotochemical heating in millisecond pulsars: modified Urca reactions with uniform Cooper pairing gaps.

# Rotochemical heating in millisecond pulsars: modified Urca reactions with uniform Cooper pairing gaps.

Cristobal Petrovich Departamento de Astronomía y Astrofísica, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile. 11email: cpetrovi@astro.puc.cl    Andreas Reisenegger Departamento de Astronomía y Astrofísica, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile. 11email: cpetrovi@astro.puc.cl
###### Key Words.:
stars: neutron — dense matter — relativity — stars: rotation — pulsars: general — pulsars: individual (PSR J0437-4715)
###### Abstract

Context:When a neutron star’s rotation slows down, its internal density increases, causing deviations from beta equilibrium that induce reactions, heating the stellar interior. This mechanism, named rotochemical heating, has previously been studied for non-superfluid neutron stars. However, the likely presence of superfluid nucleons will affect the thermal evolution of the star by suppressing the specific heat and the usual neutrino-emitting reactions, while at the same time opening new Cooper pairing reactions.

Aims:We describe the thermal effects of Cooper pairing with spatially uniform and isotropic energy gaps of neutrons and protons , on the rotochemical heating in millisecond pulsars (MSPs) when only modified Urca reactions are allowed. In this way, we are able to determine the amplitude of the superfluid energy gaps for the neutron and protons needed to produce different thermal evolution of MSPs.

Methods:We integrate numerically, and analytically in some approximate cases, the neutrino reactions for the modified Urca processes with superfluid nucleons to include them in the numerical simulation of rotochemical heating.

Results:We find that the chemical imbalances in the star grow up to the threshold value , which is higher than the quasi-steady state achieved in the absence of superfluidity. Therefore, the superfluid MSPs will take longer to reach the quasi-steady state than their nonsuperfluid counterparts, and they will have a higher a luminosity in this state, given by , where is the period derivative in units of and is the period in milliseconds. We are able to explain the UV emission of the PSR J0437-4715 for . These results are valid if the energy gaps are uniform and isotropic.

Conclusions:

## 1 Introduction

The observation of thermal emission from the surface of a neutron star (NS) has the potential to provide constraints on its inner structure. In the existing literature, several detailed cooling calculations have been compared to the few estimates available for the surface temperatures of neutron stars (see Yakovlev & Pethick 2004 for a review and references). These calculations are based on passive cooling, at first neutrino-dominated, and later driven by photon emission at the age of yr.

Several mechanisms can keep NSs hot beyond the standard cooling timescale yr. One of them is rotochemical heating, which has previously been studied for neutron stars of non-superfluid matter. It was first proposed in Reisenegger (1995) and then developed in Fernández & Reisenegger (2005) by considering the internal structure via realistic EOSs, in the framework of general relativity. It works as follows. The reduction in the centrifugal force makes the NS contract. This perturbs each fluid element, raising the local pressure and causing deviations from beta equilibrium. The system eventually reaches a new quasi-equilibrium configuration where the rate at which spin-down modifies the equilibrium concentrations is the same as that at which neutrino reactions restore the equilibrium. This implies that rotational energy is converted into thermal energy and an enhanced neutrino emission is produced by a departure from the beta equilibrium. Thus, this mechanism can keep old millisecond pulsars (MSPs) warm, at temperatures K, which implies that the surface temperature of the MSP J0437-4715 should be below its measurement (Fernández & Reisenegger 2005).

On the other hand, cooling curves usually consider the effects of nucleon superfluidity on the thermal evolution of NSs. Superfluidity is produced by Cooper pairing of baryons due to the attractive component of their interaction, and it is present only when the temperature of the matter falls below a critical temperature . The physics of these interactions is rather uncertain and very model-dependent, and so is the critical temperature obtained from theory (see Lombardo & Schulze 2001). An important microscopic effect is that the onset of superfluidity leads to the appearance of a gap in the spectrum of excitations around the Fermi surface. This considerably reduces the neutrino reactions and the specific heat involving superfluid species (Yakovlev et al. 2001), and additionally opens new neutrino emission processes, namely pair breaking and formation reactions (Flowers et al. 1976). Taking these effects into account, NS cooling has been used to constrain the amplitude of the energy gaps by comparing predictions with the surface temperatures measured from young neutron stars (see Yakovlev & Pethick 2004 for a complete discussion). In particular, Page et al. (2009) considers the minimal cooling paradigm, where no direct Urca processes are allowed and the cooling is enhanced by Cooper-pair emission processes. They find this mechanism to be consistent with observations as long as the critical temperature of the neutrons covers a range of values between K and K in the core of the star.

In this paper, we include the effects of superfluidity in rotochemical heating, building on the framework of Fernández & Reisenegger (2005) and Reisenegger et al. (2006). The order of magnitude of these effects was estimated in Reisenegger (1997), by considering the reactions to be totally forbidden, until the chemical imbalances ( ) exceed a combination of the energy gaps of either when direct Urca is allowed or when only modified Urca operates, where . This effect lengthens the equilibrium timescale and raises the surface temperature.

We consider how exactly the neutrino reactions are suppressed by superfluidity in the core, by avoiding the step-like approximation of Reisenegger (1997). We calculate instead the reduction in the net modified Urca reaction rate and the emissivity in the regime where the energy gaps, the chemical imbalances, and the temperature are all relevant quantities. In this scenario, Villain & Haensel (2005) numerically computed the phase-space integrals of the net reaction rate for direct Urca and modified Urca processes, finding that Cooper pairing does not strongly inhibit these reactions when the energy gaps are not too large compared to both the temperature and the chemical imbalances.

The structure of this paper is as follows. In Sect. 2 we present the basic equations of rotochemical heating and describe the superfluid input to compute these numerically. We explain our approach to calculating the reduction factors and how they behave in the regime of interest. In Sect. 3 we describe our results and compare our prediction with the likely thermal emission of PSR J0437-4715 (Kargaltsev et al. 2004) to constrain the values of the energy gaps. We summarize our main conclusions in Sect. 4. Finally, in Appendices A and B we explain in detail our numerical and analytical approaches to computing the net reaction rate and emissivity.

## 2 Theoretical framework

### 2.1 Rotochemical heating: basic equations

The framework of rotochemical heating used in this work is described in detail in Fernández & Reisenegger (2005). Here we just point out the fundamental equations and the modifications required to introduce the effects of superfluidity.

We consider the simplest model of a neutron star core, composed of neutrons, protons, electrons, and muons ( matter), ignoring the potential presence of exotic particles.

The internal temperature, redshifted to a distant observer, , is taken to be uniform inside the star because we are modeling the thermal evolution of a MSP over timescales much longer than the diffusion time (Reisenegger 1995). Thus, the evolution of the internal temperature for an isothermal interior is given by the thermal balance equation (Thorne 1977)

 ˙T∞ = 1C(L∞H−L∞ν−LBF,∞ν−L∞γ), (1)

where is the total heat capacity of the star, the total power released by the heating mechanism, the total power emitted as neutrinos due to Urca reactions, the total power emitted as neutrinos due to Cooper pair-breaking and pair-formation (PB-PF) processes, and the power released as thermal photons. We briefly discuss the quantities and in Sect. 2.4 and Sect. 2.5, respectively.

The amount of energy released by each Urca-type reaction is ( ). Thus, we write the total energy dissipation rate as

 L∞H=η∞npeΔ~Γnpe+η∞npμΔ~Γnpμ, (2)

where is the net reaction rate of the Urca reaction integrated over the core (indicated by the tilde) involving the lepton .

The photon luminosity is calculated by assuming black-body radiation , where and are the radius and the surface temperature of the star measured from an observer at infinity, respectively. To relate the internal and the surface temperatures, the fully accreted envelope model of Potekhin et al. (1997) is used.

The evolution of the redshifted chemical imbalances, which are also uniform throughout the core, are given by

 ˙η∞npe = −ZnpeΔ~Γnpe−ZnpΔ~Γnpμ+2WnpeΩ˙Ω (3) ˙η∞npμ = −ZnpΔ~Γnpe−ZnpμΔ~Γnpμ+2WnpμΩ˙Ω, (4)

where the terms , , , , and are constants that depend on the stellar structure and are unchanged with respect to their latest definition in Reisenegger et al. (2006), and is the product of the angular velocity and its time derivative (proportional to the spin-down power).

The key new contribution in this work is the recalculation of and , which differ substantially from those of Fernández & Reisenegger (2005), beacuse superfluidity strongly inhibits these reactions.

### 2.2 Cooper pairing

In the core, neutrons are believed to form Cooper pairs due to their interaction in the triplet state, while protons form singlet pairs. In addition, neutrons in the outermost core and inner crust are believed to form singlet-state pairs (type A) (Yakovlev et al. 2001). The (type B and C) state description is rather uncertain in the sense that the energetically most probable state of -pairs () is not known, being extremely sensitive to the still unknown -interaction (see, e.g., Amundsen & Østgaard 1985). Taking this classification into account, Villain & Haensel (2005) solve numerically the suppression caused by each type of superfluidity of the net reaction rate for direct Urca and modified Urca reactions out of beta equilibrium, finding that the suppression caused by type A superfluidity is in between the suppression by anisotropic types (type B) and (type C) superfluidity, respectively. For simplicity, we consider the energy gaps for the neutrons and the protons at zero temperature, redshifted to a distant observer, as parameters that are isotropic ( pairs) and uniform throughout the core of the NS.

The phase transition for a nucleon species into a superfluid state takes place when its temperature falls below a critical value . This temperature is related to the energy gap at zero temperature ; for the isotropic pairing channel , . In addition, when the transition occurs the amplitude of the energy gap depends on the temperature by means of the BCS equation (Yakovlev et al. 2001), which can be fitted by the practical formula of Levenfish & Yakovlev (1994) for the isotropic gap given by

 δ≡Δ(T)kT=√1−T/Tc(1.456−0.157√T/Tc+1.764T/Tc), (5)

where is the variable used in the phase-space integrals in Sect. 2.3 that depends on the input parameters and used throughout this paper. It is straightforward to check that the limiting cases are reproduced by the Eq. (5), i.e., when and when . In addition, Levenfish & Yakovlev (1994) claim that intermediate values of are also reproduced by this formula with a maximum error smaller than , which is accurate enough for the purpose of this work.

Having defined the energy gap of the nucleon with , it is possible to express the momentum dependence of the nucleon energy near the Fermi level, i.e. , as in Yakovlev et al. (2001) as

 ϵi(pi)=μi−√v2Fi(pi−pFi)2+Δ2iifpipFi, (6)

where , , and are the Fermi momentum, the Fermi velocity, and the chemical potential of species , respectively.

### 2.3 Neutrino emissivity

The most rapid reactions in NS cores are the direct Urca processes. However, as already mentioned, we assume in this work that these reactions are forbidden because direct Urca with superfluidity might drastically change the behavior of rotochemical heating, leading to new conclusions that will be discussed in a forthcoming work. Additionally, these reactions can be kinematically forbidden for a wide range of EOSs and central densities. If this were the case, the so-called modified Urca reactions (or Murca reactions) would prevail (Yakovlev et al. 2001) because they involve an additional spectator nucleon that allows energy and momentum conservation. In general, these reactions are

 n+Ni→p+Nf+e−+¯νe (7) p+Ni+e−→n+Nf+νe, (8)

where the subindices and represent the for its initial and final states of the spectator nucleon. If this nucleon is a proton, these reactions are called the proton branch of Murca, and if it is a neutron, they are called the neutron branch of Murca.

We write the neutrino emissivity and the net reaction rate due to Murca reactions involving the lepton and integrated over the core, respectively, as

 L∞ν,l = ~LnlInM,ϵT8∞+~LplIpM,ϵT8∞, (9) Δ~Γnpl = ~LnlkInM,ΓT7∞+~LplkIpM,ΓT7∞, (10)

where constants and are defined in terms of the neutrino luminosities for a nonsuperfluid NS in beta equilibrium, as

 ~Lα ≡ LeqαT8∞=∫core4πr2eΛSα(n)e−6Φdr, (11)

where indicates both the branch of the Murca process and the lepton involved in the reaction (Fernández & Reisenegger 2005). The term is a slowly varying function of the baryon number density (e.g., Yakovlev et al. 2001), and and are the usual Schwarzschild metric terms. The quantities and are dimensionless phase-space integrals that contain the dependence of the emissivity and the net reaction rate, respectively, on the chemical imbalances and on the energy gaps and .

To introduce these integrals, it is useful to define the usual dimensionless variables normalized by the thermal energy , of

 xj≡ϵj−μjkT,  xν≡ϵνkT,  and  ξl≡ηnplkT, (12)

which represent the energy of the nonsuperfluid degenerated particle , the neutrino, and the chemical imbalance involving the lepton , respectively, while for the superfluid nucleon we write

 xi≡vFi(pi−pFi)kT and zi≡% sgn(xi)√x2i+δ2i. (13)

In terms of these variables,

 INM,ϵ = 6048011513π8⋅∫∞0dxνx3ν∫∞−∞∫∞−∞∫∞−∞∫∞−∞∫∞−∞dxndxNidxNf (14) dxpdxef(zn)f(zNi)f(zNf)f(zp)f(xe) ×[δ(xν+ξl−zn−zNi−zNf−zp−xe) +δ(xν−ξl−zn−zNi−zNf−zp−xe)]

and

 INM,Γ = 6048011513π8⋅∫∞0dxνx2ν∫∞−∞∫∞−∞∫∞−∞∫∞−∞∫∞−∞dxndxNidxNf (15) dxpdxef(zn)f(zNi)f(zNf)f(zp)f(xe) ×[δ(xν+ξl−zn−zNi−zNf−zp−xe) −δ(xν−ξl−zn−zNi−zNf−zp−xe)],

where is the Fermi function and the numerical factor in front of the integral is to normalize it to 1 when the energy gaps and the chemical imbalances are zero.

In the nonsuperfluid case (i.e., ), these integrals reduce to the polynomials calculated by Reisenegger (1995)

 INM,ϵ(δn=δp=0)=FM(ξl)= 1+22020ξ2l11513π2+5670ξ4l11513π4+420ξ6l11513π6+9ξ8l11513π8, (16) INM,Γ(δn=δp=0)=HM(ξl)= 14680ξl11513π2+7560ξ3l11513π4+840ξ5l11513π6+24ξ7l11513π8, (17)

which are the same for the neutron branch and the proton branch.

Since the -variables defined in Eq. (13) depend on the energy gaps, the equality between and , or between and , is no longer satisfied if the gaps are different.

#### 2.3.1 Reduction factors

The phase-space integrals in Eqs. (14) and (15) do not have an analytical expression when one or both of the reacting nucleons are superfluid. Thus, their calculation must be done numerically, as in Villain & Haensel (2005). A natural way to account for the suppression produced by Cooper pairing is to define the so-called reduction factors as the ratio of these superfluid integrals to their non-superfluid limits

 RNM,ϵ(ξl,δn,δp) = INM,ϵ(ξl,δn,δp)/FM(ξl), (18) RNM,Γ(ξl,δn,δp) = INM,Γ(ξl,δn,δp)/HM(ξl). (19)

In principle, to calculate these reduction factors, a five-dimensional integral needs to be computed, because only one dimension can be eliminated by integrating out the electron variable in Eqs. (14) and (15). On the other hand, if one of the nucleons is not superfluid, one can eliminate more dimensions by integrating out the nonsuperfluid variables. For instance, if we consider both the neutron branch and the protons as the only superfluid species, we can integrate out the electron plus the three nonsuperfluid neutrons, obtain a two-dimensional integral that has to be calculated numerically (see Villain & Haensel 2005 for the formulae and details). However, it may be a problem to calculate these integrals efficiently including the new superfluid luminosities and net reaction rates in the rotochemical evolution equations.

The main difficulties in performing the integration of Eqs. (14) and (15) are

• infinite sizes of the integration domains;

• external free parameters , , and , which can be very large;

• many integration dimensions (up to five).

Villain & Haensel (2005) calculate these integrals numerically via the so-called Gauss-Legendre quadrature, using logarithmic variables scaled to external parameters , , . In this way they cover a wide range of the infinite integration domain, for the parameters in Eq. (19) in the range and .

We require a fast integration method because we need to evaluate eight integrals (emissivity and net reaction rate for two nucleon branches and two leptons branches) for each time-step of the thermal evolution. For this reason, we used a slightly different method than that of Villain & Haensel (2005), although we also implemented their code to calibrate ours. We chose the Gauss-Laguerre method, which is accurate enough for a few evaluations when the integrands are asymptotically exponentially decaying functions, as is the case for the Fermi functions (see Appendix A for a detailed explanation).

One striking feature of the reduction factors is that, for relatively high values of the energy gaps and chemical imbalances compared to the thermal scale, they tend to be independent of the temperature. Villain & Haensel (2005), considering only one superfluid particle species, claim by graphical inspection that, when and , the reduction factors become functions of only. By using our numerical results, we show in Appendix B that this is a good approximation if and . In the limit of zero temperature and in the presence of one superfluid variable we find the analytical expressions (see Appendix B for more details)

 RM,ϵ(Δ/η) = (28(Δ/η)2+105(Δ/η)4+1052(Δ/η)6+3516(Δ/η)8) (20) × ln(Δ/η1+√1−(Δ/η)2)+√1−(Δ/η)2 × (1+55110(Δ/η)2+432740(Δ/η)4+187380(Δ/η)6) RM,Γ(Δ/η) = (21(Δ/η)2+1052(Δ/η)4+1058(Δ/η)6) (21) × ln(Δ/η1+√1−(Δ/η)2)+√1−(Δ/η)2 × (1+75920(Δ/η)2+177940(Δ/η)4+165(Δ/η)6),

where is the energy gap of one of the nucleons, i.e., the proton if we consider the neutron branch and the neutron for the proton branch. It is straightforward to verify that these expressions satisfy the limiting cases and as expected. In addition, in the limit the functions in Eqs. (2.3) and (2.3) tend to the highest power of the polynomial, i.e., and . Thus, the integrals can be expressed as

 Iϵ(δ,ξ) = 9ξ811513π8Rϵ(δ/ξ),and (22) IΓ(δ,ξ) = 24ξ811513π8RΓ(δ/ξ), (23)

where, doing some algebra, one can verify that these approximations satisfy the identity of Flores-Tulián & Reisenegger (2006), i.e., . In Appendix B, we compare these formulae with the numerical non-zero temperature calculations, finding that they agree to a very good approximation when .

When two nucleons are superfluid, we cannot find an analytical expression, but we can simplify the problem to a three-dimensional, bounded integral. This method is obviously much faster than the finite temperature one with the quadrature integration in Appendix A. However, it is an approximation to the problem, and only works when the temperature is quite low compared to the chemical imbalance and energy gaps. Furthermore, if we analyze the case with two superfluid nucleons, the limit works when the chemical imbalance overcomes a certain energy threshold imposed by the energy gaps and . If we consider the reaction at zero temperature, the neutron is energetically allowed to decay only when , which implies the simple condition (see Reisenegger 1997 for a schematic justification of this). Finally, the threshold imposed by the gaps for both branches of Murca reactions are

 Δthr = 3Δn+Δpfor the neutron branch and, (24) Δthr = Δn+3Δpfor the proton branch. (25)

In Fig. 1, we plot the reduction factor in the net reaction rate of the neutron branch of Murca under the approximation of zero temperature. As shown in this figure, there are no reactions when . Remarkably, the contour lines of the reduction factor are almost straight lines that coincide with the slope of the contour levels of equation . This means that the reduction factor in this regime is close to a certain function of , which is valid for both branches and also for the reduction factor of the emissivity.

Since we obtain the reduction factors in the net reaction rate that are similar to those found by Villain & Haensel (2005), these can be seen in their work. The behavior of the reduction factor in the emissivity is also similar to that in the net reaction rate. Thus, to illustrate the difference between these two quantities, we plot their ratio in Fig. 2 for the range where the chemical imbalance slightly exceeds the energy gaps for different values of the temperature compared to the gaps. We show in Sect. 3 that this is the regime of interest for rotochemical heating. We can verify using the formulae in Eqs. (20) and (21) that the ratio approaches a value of one when . The ratio also decreases as the temperature diminishes compared to the energy gap until the limiting case of zero temperature. This is even more dramatic if we have two superfluid particle species. In conclusion, for very low temperatures relative to the other energy scales and chemical imbalances slightly above the energy gap, the suppression of the emissivity is much higher than that of the net reaction rate. Physically what happens is that the neutrino released in the reaction escapes with a very small amount of energy, which at zero temperature is proportional to the excess of energy that can be arbitrarily small as approaches .

We can finally incorporate the recalculation of the emissivity and net reaction rate into the rotochemical equations when one of the nucleon species is superfluid, by applying the Gauss-Laguerre quadrature method to the entire evolution. In contrast, when two superfluid nucleons are present, we separate the evolution into two regimes and follow their respective approaches. While , we use the Gauss-Laguerre quadrature method explained in Appendix A, which happens for the early evolution of rotochemical heating (see Fig. 3). When , we use the analytical approximation of the integral presented in Appendix B.

### 2.4 Specific heat suppression

When reaches , there is a discontinuous increase in the specific heat that is characteristic of a second order phase transition. Yakovlev et al. (2001) claim that this increase is by a factor of with respect to the normal-matter specific heat. Subsequently, when , an exponential-like suppression occurs because of the gap in the energy spectrum. In practice, these effects are taken into account by using control functions that multiply the unpaired values of the specific heat at constant volume . These have been calculated by several authors for various temperature regimes, e.g. Pizzochero et al. (2002) for and Maxwell (1979) for . The former set of authors obtain an exponential decaying control function of the form and the latter, which is the result we adopt, obtain the function

 CsupV,N = CV,N×3.15Tc,NTe−1.76Tc,N/T× (26) ×[2.5−1.66TTc,N+3.64(TTc,N)2],

where is the critical temperature of the nucleon and is the specific heat of the normal matter case, as defined in Fernández & Reisenegger (2005). As might be noticed, this expression also repruduces the result of Pizzochero et al. (2002) in the low-temperature regime. The value at the lower limit of the expression found by Maxwell (1979) is indeed , which is sufficiently small to ensure a the lepton contribution to the specific heat dominant. This formula also represents the discontinuous increase at , which implies that , in agreement with Yakovlev et al. (2001). Finally, an important remark is that the minimum value that the specific heat can reach is the leptonic contribution, which is .

### 2.5 Cooper pairing emission

Another feature of the superfluid state is the appearance of new neutrino reaction mechanisms due to the Cooper pairs. These are the Cooper pair breaking and pair formation processes proposed by Flowers et al. (1976). These authors claim that a superfluid neutron star can be considered as a two-component system, which consists of paired quasiparticles and elementary excitations above the condensate. Their associated quasi-equilibrium densities are controlled by the processes, which are prevalent at temperatures close to and successively suppressed at lower temperatures because of an exponential decrease in the number of unpaired particles. Schematically, these neutrino reactions are

 NN → N+N+ν+¯ν pair breaking (PB), (27) N+N → NN+ν+¯ν pair formation (PF), (28)

where denotes the Cooper pair and the excitation.
These authors find an emissivity when . Thus, these reactions could certainly affect the thermal evolution of the early stage, since the order of magnitude of the emissivity for the Murca processes is . However, the relevant question here is whether or not it affects the late stage, i.e., when the photon luminosity is the dominant cooling mechanism. In this sense, we estimate the PB-PF emissivity for a a typical scenario of rotochemical heating, where the temperature that it reaches in the late quasi-steady state is MeV (see Fernández & Reisenegger 2005). In this limit, in which (Flowers et al. 1976),

 QPB,PF∼1028(ΔMeV)7√Δ/kTe−2Δ/kT  erg cm−3 s−1, (29)

which clearly becomes negligible when the values of the energy gaps are relatively large compared to the thermal scale, say , because the exponential behavior rapidly suppresses the effects of these reactions. For MeV, we check that the PB-PF processes become irrelevant compared to the photon luminosity, and therefore do not contribute to the total luminosity at this stage.

## 3 Results and discussion

### 3.1 Evolution

We show the results of numerically solving the rotochemical evolution Eqs. (1), (3), and (4), considering the inputs in Sect. 2.

Hereafter, we model the neutron star structure by the A18 + + UIX* equation of state (EOS) (Akmal et al. 1998). The most relevant feature of this EOS in this work is that it allows direct Urca reactions for electrons above a density g cm, which corresponds to the central density of a 2 star. On the other hand, the threshold for direct Urca reactions of muons lies at a higher density in a non-causal regime. We consider stars below this mass limit, in which no direct Urca processes occur.

The results are shown in Fig. 3, where in the upper panel we plot the solution to this set of equations for the nonsuperfluid case, and in the lower we compare these to our results considering superfluidity of neutrons for a constant energy gap of MeV by integrating numerically the emissivity and net reaction rate without any approximation.

Except for the superfluidity energy gap of the neutrons in the lower panel, both panels in Fig. 3 have the same input parameters. This figure indicates that during the neutrino cooling stage, the differences between the nonsuperfluid and superfluid cases are not noticeable. After yr, the temperature drops more rapidly with time in the superfluid case because the specific heat of the neutrons is suppressed. The photon cooling, which starts to dominate in both cases, is unaltered by the presence of the Cooper pairs.

Later on, the nonsuperfluid case reaches a quasi-steady state, in which the rate at which the spin-down modifies equilibrium concentrations is the same as the rate at which the reactions drive the system toward the equilibrium concentration, with heating and cooling also balancing each other (see Reisenegger 1995 and Fernández & Reisenegger 2005 for more details). The timescale on which the superfluid case reaches this quasi-steady state is longer because the Murca reactions are highly suppressed when the chemical imbalances are smaller than . But, immediately after the chemical imbalances overcome the value of the energy gap of the neutrons, several neutrino reactions are produced, which drastically reheat the star, causing a rapid increase in the temperature to finally reach the quasi-steady state. The subsequent evolution can be approximated by the simultaneous solution of Eqs. (1), (3), and (4) with their left-hand sides set equal zero, as discussed in Sect. 3.2. Finally, the chemical imbalances at this final evolution stage reach a higher value than in the nonsuperfluid case, lengthening its timescale to reach quasi-steady state. This implies that the star can store more chemical energy and it is released it later in the time evolution compared to the nonsuperfluid case. Therefore, the heating when Cooper pairing is present is more efficient to ensure that the MSP at higher temperatures during the quasi-steady state, compared to the normal matter case. Moreover, the choice of superfluid of neutrons with 0.1 MeV and nonsuperfluid protons can explain the 90 confidence level of the surface temperature of the MSP J0437-4715 (Kargaltsev et al. 2004), which the nonsuperfluid case cannot.

This aforementioned effect of superfluidity on rotochemical heating was already predicted by Reisenegger (1997) via a rough estimation; this author claimed that the neutrino reactions opened when they overcome the threshold , that for Fig. 3 is (see Sect. 2.3.1), ignoring the temperature dependence of these reactions. By looking at in Fig. 3 and solving the evolution equations for several combinations of and in the range of MeV, we verify that for the Murca processes, the zero temperature approach of Reisenegger (1997) is valid, because in the quasi-steady state . To a good approximation, this means that the net reactions rates and the emissivity indeed do not depend on the internal temperature, as we showed in Sect. 2.3.1.

Sufficiently old stars will have reached the quasi-steady state, where , , and . In this state, the Eqs. (1), (3), and (4) can be written as

 L∞,qsγ = η∞,qsnpeΔ~Γnpe+η∞,qsnpμΔ~Γnpμ−L∞ν, (30) 2WnpeΩ˙Ω = ZnpeΔ~Γnpe+ZnpΔ~Γnpμ, (31) 2WnpμΩ˙Ω = ZnpΔ~Γnpe+ZnpμΔ~Γnpμ, (32)

where the superscript stands for quasi-steady and we have neglected the neutrino contribution of the Cooper pairing emission as explained in Sect. 2.5. This system of equations totally specifies the temperature and chemical imbalances for a star with a certain value of . Figure 4 shows the solution to this equations when only the neutrons are superfluid, such that . This illustrates that the chemical imbalance always stabilizes at a value slightly higher than ,

 η∞,qsnpe≈η∞,qsnpμ≳Δthr=min(Δn+3Δp,3Δn+Δp). (33)

We can also observe from Fig. 4 that the solution becomes closer to as we increase its value. For higher values of the gap, this is because the chemical imbalances at the quasi-steady state are higher. For a given spin-down rate, the net reaction rate is then fixed and the nonsuperfluid reactions become more efficient, increasing as . Hence, the reduction factors need to block more reactions, which happens for values of the chemical imbalance closer to . We show in Fig. 9 (Appendix B) that the reduction factors indeed become smaller as the chemical imbalance approaches the energy gap from above.

#### 3.2.1 Analytical approximation

The net reaction rates and are determined entirely by Eqs. (31) and (32). Solving these equations, we obtain

 Δ~Γnpe = 2Ω˙ΩZnpμWnpe−ZnpWnpμZnpeZnpμ−Z2np≡−KeΩ˙Ω, (34) Δ~Γnpμ = 2Ω˙ΩZnpeWnpμ−ZnpWnpeZnpeZnpμ−Z2np≡−KμΩ˙Ω, (35)

where the constants and depend exclusively on the stellar structure, i.e., the EOS and the stellar mass (see Fernández & Reisenegger 2005 for the definition of these constants).

In the quasi-steady state, in which , the polynomials in Eqs. (2.3) and (2.3) are in a very good approximation and , respectively. Thus, from the definitions (9), (10), (18), (19), and doing some algebra we substitute the expressions of Eqs. (34) and (35) to Eq. (30) to obtain

 L∞,qsγ=⎡⎣η∞,qsnpeKe⎛⎝1−924~LneRnM,ϵ(η∞,qsnpe)+~LpeRpM,ϵ(η∞,qsnpe)~LneRnM,Γ(η∞,qsnpe)+~LpeRpM,Γ(η∞,qsnpe)⎞⎠+ η∞,qsnpμKμ⎛⎝1−924~LnμRnM,ϵ(η∞,qsnpμ)+~LpμRpM,ϵ(η∞,qsnpμ)~LnμRnM,Γ(η∞,qsnpμ)+~LpμRpM,Γ(η∞,qsnpμ)⎞⎠⎤⎦∣∣Ω˙Ω∣∣, (36)

where we have just rearranged the quasi-steady Eqs. (30), (31), and (32), and no approximation has yet been made.

The first approximation that we make, as previously discussed, is to consider that the chemical imbalances are . This limit is valid when the energy gaps are relatively large, as we illustrated in Fig. 4. We checked this assumption for several solutions of the rotochemical heating evolution equations and found that it is a reasonable approximation for MeV.

Our second approximation is to neglect the ratios of the reduction factors . As we discussed in Sect. 2.3.1, the emissivity is more suppressed by the gaps than the net reaction rate when . It becomes an acceptable limit when the values of the gaps are relatively large compared to the thermal energy, say , which is the limit of interest for rotochemical heating. We checked that when we consider MeV, the ratio is .

The first simplification tends to increase the predicted luminosity, while the second approximation tends to decrease it. Thus, both effects together tend to balance each other. Finally, placing both approximations together, we obtain the simple expression for the bolometric luminosity of

 L∞,qsγ ≈ (Ke+Kμ)Δthr∣∣Ω˙Ω∣∣, (37)

whose physical interpretation is that the spin-down compression sets the number of reactions per unit time and each one of these reactions releases an amount of energy to reheat the star, as in the description of Reisenegger (1997). The approximation is even more accurate for the surface temperature, since and the relative error in the lumimosity will lead to a smaller relative error in the surface temperature.

In Fig. 5, we show the results of the predicted surface temperature in the quasi-steady state as a function of the spin-down rate by solving the equilibrium equations without approximations and using the black-body law to relate and . From this plot, we can infer that the curves are parallel and the surface temperature depends on some power of the spin-down that differs from the power of the nonsuperfluid case, which is higher. We can explain this behavior using the approximate expression , while for the nonsuperfluid case, Fernández & Reisenegger (2005) obtain .

In Fig. 6, we plot and as function of the radius for two sets of EOSs. First, we show the two more realistic EOSs from Akmal et al. (1998) for the core (A18 + and A18 + + UIX*), supplemented with that of Pethick et al. (1995) and Haensel & Pichon (1994) for the inner and outer crust, respectively. Second, we plot four representative EOSs from Prakash et al. (1988), which open direct Urca reactions for stellar masses higher than 1 M. For this set of EOSs, we show in Fig. 6 that , which, using Eq. (37), infers a bolometric luminosity of

 L∞,qsγ ≃ (1−4)×1032(ΔthrMeV% )⎛⎜⎝˙P−20P3ms⎞⎟⎠ erg% s−1, (38)

where is the period derivative in units of and is the period in milliseconds. Finally, we can express the surface temperature using the expression for the luminosity given by Eq. (37) as

 Tqss,∞ ≃ (Ke+Kμ4πσ(R∞)2)1/4Δ1/4thr∣∣Ω˙Ω∣∣1/4, (39)

where is the Stefan-Boltzmann constant. For the set of EOSs that we use, we obtain

 Tqss,∞ ≃ (5.7−6.6)×105(Δthr% MeV)1/4⎛⎜⎝˙P−20P3ms⎞⎟⎠1/4 K. (40)

#### 3.2.2 Constraints on the energy gaps

To explain the thermal emission of PSR J0437, as we showed in Fig. 3, it is necessary to invoke superfluidity, and the required values of the energy gaps then fall in a theoretically interesting range. In Fig. 7, we compare the observational allowed range of surface temperatures (Kargaltsev et al. 2004) with the theoretical predictions for different values of , using one particular EOS. In this case, we obtain the allowed range

 0.05[MeV]≲min(Δn+3Δp,3Δn+Δp)≲0.2[MeV]. (41)

By adding two more EOSs (BPAL 21, BPAL 9) that forbid direct Urca reactions in the allowed mass range, we expand the constraint to

 0.05[MeV]≲min(Δn+3Δp,3Δn+Δp)≲0.45[MeV]. (42)

These results depend of course on the assumption of spatially uniform and isotropic energy gaps. In principle, these are unrealistic assumptions because the energy gaps depend on the local density and therefore on the radius of the star, and the neutrons in the core are expected to form anisotropic Cooper pairs (see Sect. 2.2).

#### 3.2.3 Condition for arrival at the quasi-steady state

To estimate the time taken for the star to arrive at the quasi-steady state, we consider that the chemical imbalances grow by the effect of in a unimpeded way, because the reactions are highly suppressed, until they overcome the threshold imposed by the gaps. Thus, the chemical imbalance associated with the lepton evolves as

 ˙η∞npl = 2WnplΩ˙Ω, (43)

until it exceeds the threshold imposed by gaps, at the moment the quasi-steady state is reached. Thus, we set the condition for arrival at the quasi-steady state as , which can be seen to be a good approximation from the lower panel of Fig. 3. Then, integrating the expression (43) over time from the initial spin period to the current value of the spin period and using the previous condition, we obtain the upper limit to the initial spin period of

 Pi<(1P2+Δthr4π2|Wnpl|)−1/2≈(1P2ms+14(ΔthrMeV))−1/2 ms, (44)

where is the current value of the spin period in milliseconds.

For the case of PSR J0437, in particular, we conclude that the initial spin period to reach the quasi-steady state with the upper limit of the constraint in Eq. (42), i.e., MeV, should be ms. This value could in principle be compared with the initial period constraints from the cooling age of its white dwarf companion. However, this constraint is highly uncertain (Hansen & Phinney 1998), and we are therefore unable to draw definite conclusions.

## 4 Conclusions

We have studied the rotochemical heating of millisecond pulsars with only modified Urca reactions in the presence of uniform and isotropic Cooper pairing gaps for one or two nucleon species.

Based on these assumptions, we have found that the chemical imbalances in the star grow to the threshold value , which is higher than in the quasi-steady state achieved in the absence of superfluidity. Therefore, the superfluid MSPs will take longer to reach the quasi-steady state than their nonsuperfluid counterparts, and they will have a higher luminosity in this state, given by

 L∞,qsγ ≃ (1−4)×1032(ΔthrMeV% )⎛⎜⎝˙P−20P3ms⎞⎟⎠ erg% s−1. (45)

The constraint that we found for the energy gaps using our predicted effective temperatures and the black-body fit of Kargaltsev et al. (2004) to the UV emission of PSR J0437-4715 is

 0.05[MeV]≲Δthr≲0.45[MeV]. (46)

In this sense, rotochemical heating presents an interesting tool for constraining the superfluidity parameters. In a future paper, we will include the density dependence of the energy gaps via the predictions of different theoretical models.

###### Acknowledgements.
We thank Claudio Dib, Guillermo Cabrera, Márcio Catelan, Olivier Espinosa, Rafael Benguria, and Sebastián Reyes, for discussions and comments that benefited the present paper, an anonymous referee for comments that improved the final version of this paper, and Rodrigo Fernández for letting us use his rotochemical heating code. This work was supported by Proyecto Regular FONDECYT 1060644, the FONDAP Center of Astrophysics (15010003), and Proyecto Basal PFB-06/2007.

## References

• Abramowitz & Stegun (1972) Abramowitz, M. and Stegun, I. A., 1972, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (New York: Dover), 890
• Amundsen & Østgaard (1985) Amundsen, L.; Østgaard, E., 1985, Nucl. Phys. A, 442, 163.
• Akmal et al. (1998) Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804
• Davis & Rabinowitz (1975) Davis, P. & Rabinowitz P., 1975, Methods of numerical integration, (New York: Academic Press)
• Deller et al. (2008) Deller, A., Verbiest, J., Tingay, S., & Bailes, M., 2008, A&A, 685, 67
• Fernández & Reisenegger (2005) Fernández, R., Reisenegger, A. 2005, A&A, 625, 291
• Flores-Tulián & Reisenegger (2006) Flores-Tulián, S., & Reisenegger, A., 2006, MNRAS, 372, 276
• Flowers et al. (1976) Flowers, E., Ruderman, M.,&, Sutherland, P., 1976, A&A., 205, 541
• Haensel & Pichon (1994) Haensel, P., & Pichon, B. 1994, A&A, 283, 313
• Hansen & Phinney (1998) Hansen, B. M. S., & Phinney, E. S. 1998, MNRAS, 294, 569
• Kargaltsev et al. (2004) Kargaltsev, O., Pavlov, G. G., & Romani, R. 2004, A&A, 602, 327
• Levenfish & Yakovlev (1994) Levenfish, K. P., & Yakovlev, D. G. 1994, Astronomy Reports, 38, 247
• Lombardo & Schulze (2001) Lombardo, U., & Schulze, H., 2001, Lecture Notes in Physics, 578, 30
• Maxwell (1979) Maxwell, O.V., 1979, ApJ, 231, 201
• Page et al. (2009) Page, D., Lattimer, J., Prakash, M., & Steiner, 2009, ApJ, 707, 1131
• Pethick et al. (1995) Pethick, C. J., Ravenhall, D. G., & Lorentz, C. P. 1995, Nucl. Phys. A, 584, 675
• Pizzochero et al. (2002) Pizzochero, P. M.,2002, A&A, 569, 381
• Potekhin et al. (1997) Potekhin, A. Y., Chabrier, G., & Yakovlev, D. G. 1997, A&A, 323, 415
• Prakash et al. (1988) Prakash, M., Ainsworth, T. L., & Lattimer, J. M. 1988, Phys. Rev. Lett., 61, 2518
• Reisenegger (1995) Reisenegger, A. 1995, A&A, 442, 749
• Reisenegger (1997) Reisenegger, A. 1997, A&A, 485, 313
• Reisenegger et al. (2006) Reisenegger, A., Jofré, P., Fernández, R. & Kantor, E. 2006, A&A, 653, 568
• Takahasi & Mori (1973) Takahasi, H. & Mori, M. , 1973, PRIMS, 9, 721
• Thorne (1977) Thorne, K. S. 1977, ApJ, 212, 825
• Villain & Haensel (2005) Villain, L., Haensel, P., 2005, A&A, 444, 539
• Verbiest et al. (2008) Verbiest, J., et al. 2008, A&A, 679, 675
• Yakovlev et al. (2001) Yakovlev, D. G. 2001, Phys. Rep, 354, 1
• Yakovlev & Pethick (2004) Yakovlev, D. G., & Pethick, C. J. 2004, ARA&A, 42, 169

## Appendix A Gauss-Laguerre integration method

We describe the numerical method used in this work to calculate the phase-space integrals involved in the net reaction rate and emissivity of the beta processes considering superfluid nucleon species. Villain & Haensel (2005) showed that Cooper pairing reduces these reactions and that there is no analytical expression to compute this suppression. These authors, therefore, solved the phase-space integrals numerically using the so-called Gauss-Legendre quadrature, logarithmic variables, and cut-off values to cover a wide range of the unbounded integration domain. In this work, it was more convenient for us to use the method presented here because it involves fewer evaluations than the method used in Villain & Haensel (2005) for a reasonable precision.

The method we show is motivated by the behavior of the integrand for the beta processes, which decays exponentially with each variable because of the Fermi functions. This happens asymptotically even when some particles are superfluid. For this reason, among the Gaussian quadratures, a natural candidate for an interpolation polynomial are the so-called Laguerre polynomials because these are defined on the basis of the continuous functions in , whose inner product is defined with weight function (Abramowitz & Stegun 1972). If the function to be integrated is the product of an th-order polynomial and an exponential function, this method is exact using an nth-order Laguerre polynomial in the Gaussian quadrature. Therefore, a necessary condition for this method to work properly for an integral is that has to be a smooth function, such as a polynomial (Davis & Rabinowitz 1975).

We present how this method has been applied to the phase-space integrals involved in the net reaction rate for the direct Urca process to save notation. An extension to the Modified Urca processes is straightforward by adding two more superfluid variables to the subsequent expressions. Considering the limits in the positive domain of each integration variable the integral is given by

 ID,Γ = ∫∞0dxνx2ν∫∞0∫∞0dxndxp×⎧⎨⎩∑jn=±1∑jp=±1f