Rotochemical heating in millisecond pulsars: modified Urca reactions with uniform Cooper pairing gaps
C. Petrovich^{}  A. Reisenegger
Departamento de Astronomía y Astrofísica, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile
Received 13 December 2009 / Accepted 3 March 2010
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
nonsuperfluid 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 neutrinoemitting
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 quasisteady state achieved in the absence of
superfluidity. Therefore, the superfluid MSPs will take longer to reach
the quasisteady 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 10^{20} and
is the period in milliseconds. We are able to explain the UV emission of the PSR J04374715 for
.
These results are valid if the energy gaps are uniform and isotropic.
Key words: stars: neutron  dense matter  stars: rotation  pulsars: general  pulsars: individual: PSR J04374715
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 neutrinodominated, 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 nonsuperfluid 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 quasiequilibrium configuration where the rate at which spindown 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 J04374715 should be 20% 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 Tof the matter falls below a critical temperature . The physics of these interactions is rather uncertain and very modeldependent, 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 Cooperpair emission processes. They find this mechanism to be consistent with observations as long as the critical temperature ^{3}P_{2} 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 steplike 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 phasespace 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 J04374715 (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)
where C 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 pairbreaking and pairformation (PBPF) processes, and the power released as thermal photons. We briefly discuss the quantities C and [0pt] in Sects. 2.4 and 2.5, respectively.
The amount of energy released by each Urcatype reaction is
().
Thus, we write the total energy dissipation rate as
where is the net reaction rate of the Urca reaction integrated over the core (indicated by the tilde) involving the lepton l.
The photon luminosity is calculated by assuming blackbody 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
where the terms Z_{np}, Z_{npe}, , W_{npe}, 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 spindown 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 ^{3}P_{2} state, while protons form singlet ^{1}S_{0} pairs. In addition, neutrons in the outermost core and inner crust are believed to form singletstate ^{1}S_{0} pairs (type A) (Yakovlev et al. 2001). The ^{3}P_{2} (type B and C) state description is rather uncertain in the sense that the energetically most probable state of nnpairs ( ) is not known, being extremely sensitive to the still unknown nninteraction (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 (^{1}S_{0} 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 ^{1}S_{0},
.
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 befitted by the practical formula of
Levenfish & Yakovlev (1994) for the isotropic gap given by
where is the variable used in the phasespace 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 5%, which is accurate enough for the purpose of this work.
Having defined the energy gap of the nucleon
with i=n,p, 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
where p_{Fi}, v_{Fi}, and are the Fermi momentum, the Fermi velocity, and the chemical potential of species i=n,p, 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. In additional these reactions can be kinematically forbidden for a wide range of
EOSs and central densities.
If this were the case, the socalled modified Urca reactions (or Murca reactions)
would prevail (Yakovlev et al. 2001) because they involve
an additional spectator nucleon N that allows energy and momentum conservation.
In general, these reactions are
(7)  
(8) 
where the subindices i and f represent the 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 l and
integrated over the core, respectively, as
where constants and are defined in terms of the neutrino luminosities for a nonsuperfluid NS in beta equilibrium, as
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 n (e.g., Yakovlev et al. 2001), and and are the usual Schwarzschild metric terms. The quantities and are dimensionless phasespace 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 kT, of
which represent the energy of the nonsuperfluid degenerated particle j, the neutrino, and the chemical imbalance involving the lepton l, respectively, while for the superfluid nucleon i we write
In terms of these variables,
and
where is the Fermi function f(x)=1/(1+e^{x}) 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)
which are the same for the neutron branch and the proton branch.
Since the zvariables 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 phasespace 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 socalled reduction factors as the ratio of these superfluid integrals
to their nonsuperfluid limits
In principle, to calculate these reduction factors, a fivedimensional 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 twodimensional 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).
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 timestep 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 GaussLaguerre 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)
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. (16) and (17) tend to the highest power of the polynomial, i.e., and . Thus, the integrals can be expressed as
=  (22)  
=  (23) 
where, doing some algebra, one can verify that these approximations satisfy the identity of FloresTulián & Reisenegger (2006), i.e., . In Appendix B, we compare these formulae with the numerical nonzero temperature calculations, finding that they agree to a very good approximation when .
Figure 1: Contour plot of the reduction factor in the net reaction rate in the neutron branch of Murca in the limit of zero temperature. The numbers are . 

Open with DEXTER 
When two nucleons are superfluid, we cannot find an analytical
expression, but we can simplify the problem to a threedimensional,
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 schematicjustification of this).
Finally, the threshold imposed by the gaps
for both branches
of Murca reactions are
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.
Figure 2: Ratio of the reduction factor of the emissivity to the reduction factor of the net reaction rate for the proton branch of Murca reactions with superfluid neutrons as a function of for different values of . Ratio of the zero temperature approximations given by Eqs. (20) and (21) is given by the black dashdotted line. 

Open with DEXTER 
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 GaussLaguerre 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 GaussLaguerre 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 T 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 2.4 with respect
to the normalmatter specific heat.
Subsequently, when
,
an exponentiallike 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
where is the critical temperature of the nucleon N 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 lowtemperature 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 twocomponent system,
which consists of paired quasiparticles and
elementary excitations above the condensate.
Their associated quasiequilibrium 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  (27)  
N+N  (28) 
where NN denotes the Cooper pair and N 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 PBPF emissivity for a a typical scenario of rotochemical heating, where the temperature that it reaches in the late quasisteady state is MeV (see Fernández & Reisenegger 2005). In this limit, in which (Flowers et al. 1976),
(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 PBPF 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^{3}, 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 noncausal 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.
Figure 3: Evolution of the internal temperature , the surface temperature and the chemical imbalances , for a star with the mass fixed to the PSR J04374715, i.e. (Verbiest et al. 2008), built with the A18 + + UIX* EOS, with the initial condition K and null chemical imbalances. The spindown is assumed to be due to the magnetic dipole radiation, with dipolar field strength G and initial period of P_{0}=1 ms. The error bar is the confidence level for the surface temperature measured for the PSR J04374715 (Kargaltsev et al. 2004) at the current spindown parameters. Upper panel: nonsuperfluid case (null energy gaps). Lower panel: superfluidity of neutrons with MeV (dashed line). 

Open with DEXTER 
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 10^{5} 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 quasisteady state, in which the rate at which the spindown 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; Fernández & Reisenegger 2005, for more details). The timescale on which the superfluid case reaches this quasisteady 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 quasisteady state. The subsequent evolution can be approximated by the simultaneous solution of Eqs. (1), (3), and (4) with their lefthand 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 quasisteady 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 quasisteady 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 J04374715 (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 quasisteady 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.
3.2 Quasisteady state
Sufficiently old stars will have reached the quasisteady state,
where
,
,
and
.
In this state, the Eqs. (1), (3),
and (4) can be written as
where the superscript qs stands for quasisteady 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 ,
(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 quasisteady state are higher. For a given spindown 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. B.1 (Appendix B) that the reduction factors indeed become smaller as the chemical imbalance approaches the energy gap from above.
Figure 4: Chemical imbalance of the electrons in the quasisteady state with only neutrons as superfluid particles as a function of for three values of the energy gap (solid lines) and the nonsuperfluid case (dashed line) for a star calculated with the A18 + + UIX* EOS. The horizontal lines indicate the values of the energy gaps: 0.01,0.05, and 0.1 MeV (dotted lines). 

Open with DEXTER 
3.2.1 Analytical approximation
The net reaction rates
and
are determined entirely by Eqs. (31) and (32).
Solving these equations, we obtain
where the constants K_{e} 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 quasisteady state, in which
,
the polynomials in Eqs. (16) and
(17) 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
where we have just rearranged the quasisteady Eqs. (30)(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
whose physical interpretation is that the spindown 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 quasisteady state as a function of the spindown rate by solving the equilibrium equations without approximations and using the blackbody 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 spindown 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 .
Figure 5: Surface temperature predicted for the quasisteady state for several values of (solid lines) and the nonsuperfluid case (dashdotted line) as a function of the spindown rate with the same stellar parameters as in Fig. 4. 

Open with DEXTER 
Figure 6: Values of K_{e} (solid lines) and (dashdotted lines) for different EOSs for a range of masses from 1 to the mass at which direct Urca with electrons is opened for each EOS. 

Open with DEXTER 
In Fig. 6, we plot K_{e} 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 .
For this set of EOSs, we show in Fig. 6 that
,
which,
using Eq. (37), infers a bolometric luminosity of
where is the period derivative in units of 10^{20} and is the period in milliseconds. Finally, we can express the surface temperature using the expression for the luminosity given by Eq. (37) as
where is the StefanBoltzmann constant. For the set of EOSs that we use, we obtain
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
By adding two more EOSs (BPAL 21, BPAL 9) that forbid direct Urca reactions in the allowed mass range, we expand the constraint to
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).
Figure 7: Surface temperature predicted in the quasisteady state as a function of the the stellar radius for a star built with the A18 + + UIX* EOS for several values of in units of MeV (solid lines) and the current spindown rate of the MSP J0437. Vertical dotted lines indicates the mass range allowed for MSP J0437 (Verbiest et al. 2008). Dasheddotted lines are the 90 confidence contours of the blackbody fit of Kargaltsev et al. (2004) to probable thermal emission from this pulsar corrected for the latest distance of 157 pc (Deller et al. 2008). 

Open with DEXTER 
3.2.3 Condition for arrival at the quasisteady state
To estimate the time taken for the star to arrive at the quasisteady
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 l evolves as
until it exceeds the threshold imposed by gaps, at the moment the quasisteady state is reached. Thus, we set the condition for arrival at the quasisteady 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 P_{i} to the current value of the spin period P and using the previous condition, we obtain the upper limit to the initial spin period of
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 quasisteady 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 quasisteady
state achieved in the absence of superfluidity. Therefore, the superfluid MSPs
will take longer to reach the quasisteady state than their nonsuperfluid
counterparts, and they will have a higher luminosity in this state, given
by
The constraint that we found for the energy gaps using our predicted effective temperatures and the blackbody fit of Kargaltsev et al. (2004) to the UV emission of PSR J04374715 is
(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 (González et al., in prep.). 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, 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 PFB06/2007.
Appendix A: GaussLaguerre integration method
We describe the numerical method used in this work to calculate the phasespace 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 phasespace integrals numerically using the socalled GaussLegendre quadrature, logarithmic variables, and cutoff 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 socalled Laguerre polynomials because these are defined on the basis of the continuous functions in , whose inner product is defined with weight function W(x)=e^{x} (Abramowitz & Stegun 1972). If the function to be integrated is the product of an nthorder polynomial and an exponential function, this method is exact using an nthorder 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).
Figure A.1: Relative error in the integral for the Legendre quadrature and Laguerre quadrature as a function of the number of evaluations considered in each dimension for the integral of the direct Urca net reaction rate with . 

Open with DEXTER 
We present how this method has been applied to the phasespace 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
where (i=n,p) and is the Fermi function. In what follows, we argue that this integral satisfies the Laguerre quadrature condition explained above, and we then express the numerical formula to compute it.
The integrand is an exponentially decaying function of the neutrino variable, and thus, it satisfies the quadrature condition. The integrand for the nucleons has the shape of a multiplication of two Fermi functions. For instance, if we consider the neutron variable we have , where depends on the neutrino variable, the proton variable, and the chemical imbalance. It is easy to see that the asymptotic behavior of this function is exponential by separating both cases and imposing the limiting case , which is the relevant scale for this variable. Thus, in this limit and , while the other Fermi functions have the opposite behaviors and , respectively. Therefore, the nucleons also satisfy the quadrature condition, which may, however, also be applicable to higher values of these variables depending on the values that take. In this sense, higher values of , which depends on the inputs , , and , will make the quadrature less accurate for a given amount of point evaluations.
An additional step in applying this quadrature method is to define a
relevant scale for each integration variable, as the scale of
the argument of the Fermi function
in the integrand of Eq. (A.1), as in Villain & Haensel (2005).
For the neutrino variable, this scale is
,
and for the
superfluid nucleons they are
and
.
The numerical formula is, finally
where and . The values , x_{in}, and x_{ip} are the roots of the Laguerre polynomials of order , n_{n}, and n_{p}, multiplied by , S_{xn}, and S_{xp}, respectively . The values , W_{in}, and W_{ip} are the weight factors of the Laguerre polynomials, which can easily be obtained from tables or using their usual recursive formulae.
In practice, this method is used for Modified Urca processes where there are four superfluid nucleon variables and the integration dimensions increase to five.
In Fig. A.1, we provide an example where the Laguerre quadrature is more accurate than the Legendre quadrature over a few evaluations. However, to achieve high precision the Legendrelike integration method is preferable because it can be used to perform many evaluations raising the number of point evaluations, while it is more difficult to compute the roots and weights of the Laguerre method using recursive methods to high orders of these polynomials.
Appendix B: Analytical approximation
In the context of rotochemical heating, most of the MSP's lifetime
is in a regime where the chemical imbalance is much higher
than the temperature (
).
On invoking superfluidity, the chemical imbalances increase
until they become of the order of the superfluidity energy gaps.
Thus, a relevant case to analyze and calculate in this work is the limit
(B.1) 
The first analysis is given for the simplest case: the direct Urca process. The idea is to extend this analysis to the most general and complicated case: modified Urca with two superfluid nucleons. After integrating over the electron variable, the phasespace integral for the net reaction rate is
where , with i=n,p.
In the regime of interest, the Fermi functions behave almost like
a step function, since the thermal scale is negligible compared
to the other relevant scales.
Moreover, one of the conclusions
of Villain & Haensel (2005), obtained by means of their numerical results,
is that for large enough energy gaps, say
,
and as soon
as
,
the reduction factors of the net reaction rates for
direct Urca and modified Urca reactions become very nearly independent of the temperature.
This motivates us to take the limit of zero temperature.
Thus, replacing the Fermi functions by step functions, i.e.,
,
the integral gives
Rearranging the limits of integration and eliminating the function from , the integral can be reduced to
By applying these approximations, the reactions can be seen to be open only when . (see Reisenegger 1997, for a schematic justification).
The integral in Eq. (B.4) appears not to have a general, closed expression in terms of the three parameters , but it is possible to perform two out of three variable integrations. By considering only one superfluid nucleon, an analytical expression can be found. In what follows, we explain these results.
The first step is to define an appropriate change of variables
where the Jacobian of this transformation is , considering that , and the new limits of integration must satisfy the upper limit imposed by the step function and the lower limit defined by the change of variables. Thus, the new integral is
Integrating over the neutrino variable, we obtain
The integration over the proton variable yields
where
The latter integral in Eq. (B.9) appears not to have an analytical closed formula. For this reason, the solution hereafter is taken numerically. The method for the emissivity integral is completely analogous, but with the only difference that we replace with in Eq. (B.4). The formula obtained in this case is
where
The final step is to extend the same reasoning to the modified Urca processes, and is natural on the basis of the previous analysis. Here one has to distinguish between the neutron branch (superscript n) and the proton branch (superscript p), but for the similarity relation only one expression is needed. Thus, without loss of generality, we only write the neutron branch expression of the net reaction rate of Murca
where is the function given by Eq. (B.9) and n_{i} and n_{f} represent the initial and final neutron spectator, respectively. In the same way, the emissivity has an expression in terms of the function defined by Eq. (B.11):
In both cases, the net reaction rate and the emissivity, the condition to recover the nonsuperfluid case asymptotically is satisfied.
The advantage of this treatment is that the numerical integration of the net reaction rate and the emissivity for the Urca processes is much faster than those without these approximations. First, the region of integration is bounded, while those of the initial integrals are not. This reduces considerably the number of points used in each integral. Secondly, the dimensions to integrate are reduced to one in the case of direct Urca and to three in the case of modified Urca. In practice, to compute these integrals we use the tanhsinh method (Takahasi & Mori 1973), which integrates the singularities in and/or in quite rapidly.
Figure B.1: Reduction factor R for the net reaction rate of direct Urca with one superfluid particle as a function of the ratio for several values of the ratio (dashed lines) using the method detailed in Appendix A, and the limit case calculated from Eqs. (B.15) and (B.16) (solid line). 

Open with DEXTER 
The integrals in Eqs. (B.8), (B.10),
(B.12), and (B.13)
can be solved analytically when only one superfluid reactant is present, which is
characterized by the energy gap .
The reduction factor R (for all the previous cases) obtained from this simplification
is given by the integral
where and for Durca processes (netreaction rate and emissivity, respectively), and and for Murca processes (netreaction rate and emissivity, respectively). The solution to this integral is
where and are the polynomials
Finally, we show in Fig. B.1 a comparison between our zero temperature approximation and the exact calculations for finite values of temperature, which demonstrates the accuracy of our approximation even for values of the energy gap compared that are not very large relative to the temperature, say .
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (New York: Dover), 890 (In the text)
 Amundsen, L., & Østgaard, E. 1985, Nucl. Phys. A, 442, 163 [NASA ADS] [CrossRef] (In the text)
 Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804 [NASA ADS] [CrossRef] (In the text)
 Davis, P., & Rabinowitz P. 1975, Methods of numerical integration, (New York: Academic Press) (In the text)
 Deller, A., Verbiest, J., Tingay, S., & Bailes, M., 2008, A&A, 685, 67 (In the text)
 Fernández, R., & Reisenegger, A. 2005, A&A, 625, 291 [NASA ADS] [CrossRef] (In the text)
 FloresTulián, S., & Reisenegger, A. 2006, MNRAS, 372, 276 [NASA ADS] [CrossRef] (In the text)
 Flowers, E., Ruderman, M., & Sutherland, P. 1976, A&A, 205, 541 [NASA ADS] [CrossRef] (In the text)
 Haensel, P., & Pichon, B. 1994, A&A, 283, 313 [NASA ADS] (In the text)
 Hansen, B. M. S., & Phinney, E. S. 1998, MNRAS, 294, 569 [NASA ADS] [CrossRef] (In the text)
 Kargaltsev, O., Pavlov, G. G., & Romani, R. 2004, A&A, 602, 327 [NASA ADS] [CrossRef] (In the text)
 Levenfish, K. P., & Yakovlev, D. G. 1994, Astron. Rep., 38, 247 [NASA ADS] (In the text)
 Lombardo, U., & Schulze, H. 2001, Lect. Not. Phys., 578, 30 [NASA ADS] [CrossRef] (In the text)
 Maxwell, O. V. 1979, ApJ, 231, 201 [NASA ADS] [CrossRef] (In the text)
 Page, D., Lattimer, J. M., Prakash, M., & Steiner, A. W. 2009, ApJ, 707, 1131 [NASA ADS] [CrossRef] (In the text)
 Pethick, C. J., Ravenhall, D. G., & Lorentz, C. P. 1995, Nucl. Phys. A, 584, 675 [NASA ADS] [CrossRef] (In the text)
 Pizzochero, P. M. 2002, A&A, 569, 381 (In the text)
 Potekhin, A. Y., Chabrier, G., & Yakovlev, D. G. 1997, A&A, 323, 415 [NASA ADS] (In the text)
 Prakash, M., Ainsworth, T. L., & Lattimer, J. M. 1988, , 61, 2518 (In the text)
 Reisenegger, A. 1995, A&A, 442, 749 [NASA ADS] [CrossRef] (In the text)
 Reisenegger, A. 1997, A&A, 485, 313 [NASA ADS] [CrossRef] (In the text)
 Reisenegger, A., Jofré, P., Fernández, R., & Kantor, E. 2006, A&A, 653, 568 [NASA ADS] [CrossRef] (In the text)
 Takahasi, H., & Mori, M. 1973, PRIMS, 9, 721 (In the text)
 Thorne, K. S. 1977, ApJ, 212, 825 [NASA ADS] [CrossRef] (In the text)
 Villain, L., & Haensel, P. 2005, A&A, 444, 539 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Verbiest, J. P. W., Bailes, M., van Straten, W., et al. 2008, A&A, 679, 675 (In the text)
 Yakovlev, D. G. 2001, Phys. Rep., 354, 1 [NASA ADS] [CrossRef] (In the text)
 Yakovlev, D. G., & Pethick, C. J. 2004, ARA&A, 42, 169 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
Footnotes
 ... Petrovich^{}
 Present address: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA
All Figures
Figure 1: Contour plot of the reduction factor in the net reaction rate in the neutron branch of Murca in the limit of zero temperature. The numbers are . 

Open with DEXTER  
In the text 
Figure 2: Ratio of the reduction factor of the emissivity to the reduction factor of the net reaction rate for the proton branch of Murca reactions with superfluid neutrons as a function of for different values of . Ratio of the zero temperature approximations given by Eqs. (20) and (21) is given by the black dashdotted line. 

Open with DEXTER  
In the text 
Figure 3: Evolution of the internal temperature , the surface temperature and the chemical imbalances , for a star with the mass fixed to the PSR J04374715, i.e. (Verbiest et al. 2008), built with the A18 + + UIX* EOS, with the initial condition K and null chemical imbalances. The spindown is assumed to be due to the magnetic dipole radiation, with dipolar field strength G and initial period of P_{0}=1 ms. The error bar is the confidence level for the surface temperature measured for the PSR J04374715 (Kargaltsev et al. 2004) at the current spindown parameters. Upper panel: nonsuperfluid case (null energy gaps). Lower panel: superfluidity of neutrons with MeV (dashed line). 

Open with DEXTER  
In the text 
Figure 4: Chemical imbalance of the electrons in the quasisteady state with only neutrons as superfluid particles as a function of for three values of the energy gap (solid lines) and the nonsuperfluid case (dashed line) for a star calculated with the A18 + + UIX* EOS. The horizontal lines indicate the values of the energy gaps: 0.01,0.05, and 0.1 MeV (dotted lines). 

Open with DEXTER  
In the text 
Figure 5: Surface temperature predicted for the quasisteady state for several values of (solid lines) and the nonsuperfluid case (dashdotted line) as a function of the spindown rate with the same stellar parameters as in Fig. 4. 

Open with DEXTER  
In the text 
Figure 6: Values of K_{e} (solid lines) and (dashdotted lines) for different EOSs for a range of masses from 1 to the mass at which direct Urca with electrons is opened for each EOS. 

Open with DEXTER  
In the text 
Figure 7: Surface temperature predicted in the quasisteady state as a function of the the stellar radius for a star built with the A18 + + UIX* EOS for several values of in units of MeV (solid lines) and the current spindown rate of the MSP J0437. Vertical dotted lines indicates the mass range allowed for MSP J0437 (Verbiest et al. 2008). Dasheddotted lines are the 90 confidence contours of the blackbody fit of Kargaltsev et al. (2004) to probable thermal emission from this pulsar corrected for the latest distance of 157 pc (Deller et al. 2008). 

Open with DEXTER  
In the text 
Figure A.1: Relative error in the integral for the Legendre quadrature and Laguerre quadrature as a function of the number of evaluations considered in each dimension for the integral of the direct Urca net reaction rate with . 

Open with DEXTER  
In the text 
Figure B.1: Reduction factor R for the net reaction rate of direct Urca with one superfluid particle as a function of the ratio for several values of the ratio (dashed lines) using the method detailed in Appendix A, and the limit case calculated from Eqs. (B.15) and (B.16) (solid line). 

Open with DEXTER  
In the text 
Copyright ESO 2010