EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 521 (October 2010) A&A, 521 (2010) A77 Full HTML
Free Access
Issue
A&A
Volume 521, October 2010
Article Number A77
Number of page(s) 12
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/200913861
Published online 22 October 2010
A&A 521, A77 (2010)

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 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 $\Delta _n$ and protons $\Delta_p$, 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 $\Delta_{\rm thr}=\mbox{min}\left( \Delta_n+3\Delta_p,3\Delta_n+\Delta_p \right)$, 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 $L_\gamma^{\infty,qs} = \left(1{-}4\right) \times 10^{32}\left(\Delta_{\rm thr}/{\rm MeV}\right) \left(\dot{P}_{-20}/P_{{\rm ms}}^3\right)~{\rm erg}~{\rm s}^{-1}$, where $\dot{P}_{-20}$ is the period derivative in units of 10-20 and $P_{{\rm ms}}$ is the period in milliseconds. We are able to explain the UV emission of the PSR J0437-4715 for $0.05~[ \mbox{MeV}]\lesssim\Delta_{\rm thr}\lesssim 0.45~[ \mbox{MeV}]$. 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 J0437-4715

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 ${\sim}10^5$ yr.

Several mechanisms can keep NSs hot beyond the standard cooling timescale ${\sim}10^7$ 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 ${\sim}10^5$ K, which implies that the surface temperature of the MSP J0437-4715 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 $T_{\rm c}$. 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 $\Delta$ 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 3P2 of the neutrons covers a range of values between $T_{\rm c}^{\min}\lesssim0.2\times10^9$ K and $T_{\rm c}^{\max}\gtrsim0.5\times10^9$ 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 $\eta_{npl}= \mu_n- \mu_p- \mu_l$ ($l=e,\mu$) exceed a combination of the energy gaps of either $\Delta_n+\Delta_p$ when direct Urca is allowed or $\Delta_n+\Delta_p+2\Delta_<$ when only modified Urca operates, where $\Delta_<=\min (\Delta_n,\Delta_p )$. 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 ($npe\mu$ matter), ignoring the potential presence of exotic particles.

The internal temperature, redshifted to a distant observer, $T_\infty $, 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)

$\displaystyle \dot{T}_\infty = \frac{1}{C} \left( L_{\rm H}^{\infty}-L_\nu^\infty -L_\nu^{{\rm BF},\infty}- L_\gamma^\infty\right),$     (1)

where C is the total heat capacity of the star, $L_{\rm H}^{\infty}$ the total power released by the heating mechanism, $L_\nu^\infty$ the total power emitted as neutrinos due to Urca reactions, $L_\nu^{{\rm BF},\infty}$ the total power emitted as neutrinos due to Cooper pair-breaking and pair-formation (PB-PF) processes, and $L_\gamma^\infty$ the power released as thermal photons. We briefly discuss the quantities C and [0pt] $L_\nu^{{\rm BF},\infty}$ in Sects. 2.4 and 2.5, respectively.

The amount of energy released by each Urca-type reaction is $\eta_{npl}= \mu_n- \mu_p- \mu_l$ ($l=e,\mu$). Thus, we write the total energy dissipation rate as

$\displaystyle L_{\rm H}^{\infty}= \eta_{npe}^{\infty}\Delta\tilde{\Gamma}_{npe}+\eta_{np\mu}^{\infty}\Delta\tilde{\Gamma}_{np\mu},$     (2)

where $\Delta\tilde{\Gamma}_{npl}=\tilde{\Gamma}_{n\rightarrow pl}-\tilde{\Gamma}_{pl\rightarrow n}$ 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 black-body radiation $L_\gamma^\infty=4\pi\sigma (R^\infty)^2 T_{s,\infty}^4$, where $R^\infty$ and $T_{s,\infty }$ 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

                           $\displaystyle \dot{\eta}^\infty_{npe}$ = $\displaystyle -Z_{npe}\Delta\tilde{\Gamma}_{npe} - Z_{np}\Delta\tilde{\Gamma}_{np\mu} + 2W_{npe}\Omega \dot{\Omega}$ (3)
$\displaystyle \dot{\eta}^\infty_{np\mu}$ = $\displaystyle -Z_{np}\Delta\tilde{\Gamma}_{npe} - Z_{np\mu}\Delta\tilde{\Gamma}_{np\mu} + 2W_{np\mu}\Omega \dot{\Omega},$ (4)

where the terms Znp, Znpe, $Z_{np\mu}$, Wnpe, and $W_{np\mu}$are constants that depend on the stellar structure and are unchanged with respect to their latest definition in Reisenegger et al. (2006), and $\Omega \dot{\Omega}$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 $L_\nu^\infty$ and $\Delta\tilde{\Gamma}_{npe}$, 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 3P2 state, while protons form singlet 1S0 pairs. In addition, neutrons in the outermost core and inner crust are believed to form singlet-state 1S0 pairs (type A) (Yakovlev et al. 2001). The 3P2 (type B and C) state description is rather uncertain in the sense that the energetically most probable state of nn-pairs ( $\vert m_{\rm J}\vert=0,1,2$) is not known, being extremely sensitive to the still unknown nn-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 $\vert m_{\rm J}\vert=0$ (type B) and $\vert m_{\rm J}\vert=2$ (type C) superfluidity, respectively. For simplicity, we consider the energy gaps for the neutrons $\Delta _n$ and the protons $\Delta_p$ at zero temperature, redshifted to a distant observer, as parameters that are isotropic (1S0 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 $T_{\rm c}$. This temperature is related to the energy gap at zero temperature $\Delta(T=0)$; for the isotropic pairing channel 1S0, $\Delta(T=0)=1.764 k T_{\rm c}$. 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

$\displaystyle \delta \equiv \frac{\Delta(T)}{kT}=\sqrt{1-T/T_{\rm c}}\left(1.456-\frac{0.157}{\sqrt{T/T_{\rm c}}}+\frac{1.764}{T/T_{\rm c}}\right),$   (5)

where $\delta$ is the variable used in the phase-space integrals in Sect. 2.3 that depends on the input parameters $\Delta _n$ and $\Delta_p$ used throughout this paper. It is straightforward to check that the limiting cases are reproduced by the Eq. (5), i.e., $\delta=0$ when $T=T_{\rm c}$ and $\delta=\Delta(T=0)/kT$ when $T\ll T_{\rm c}$. In addition, Levenfish & Yakovlev (1994) claim that intermediate values of $T/T_{\rm c}$ 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 $\Delta_i$ with i=n,p, it is possible to express the momentum dependence of the nucleon energy $\epsilon_i(p_i)$ near the Fermi level, i.e. $\vert p_i-p_{F_i}\vert\ll p_{F_i}$, as in Yakovlev et al. (2001) as

$\displaystyle \epsilon_i(p_i)=\mu_i- \sqrt{v_{F_i}^2(p_i-p_{F_i})^2+\Delta_i^2}\quad\mbox{if}\quad p_i<p_{F_i} ,$    
$\displaystyle \epsilon_i(p_i)=\mu_i+ \sqrt{v_{F_i}^2(p_i-p_{F_i})^2+\Delta_i^2}\quad\mbox{if}\quad p_i>p_{F_i},$   (6)

where pFi, vFi, and $\mu_i$ 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 so-called 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

$\displaystyle n+N_i\rightarrow p+N_f+{\rm e}^-+\bar{\nu}_e$   (7)
$\displaystyle p+N_i+{\rm e}^-\rightarrow n+N_f+\nu_e,$   (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

                     $\displaystyle L_{\nu,l}^\infty$ = $\displaystyle \tilde{L}_{nl}I_{M,\epsilon}^n T_\infty^8 +\tilde{L}_{pl}I_{M,\epsilon}^p T_\infty^8,$ (9)
$\displaystyle \Delta\tilde{\Gamma}_{npl}$ = $\displaystyle \frac{\tilde{L}_{nl}}{k}I_{M,\Gamma}^n T_\infty^7 +\frac{\tilde{L}_{pl}}{k}I_{M,\Gamma}^p T_\infty^7,$ (10)

where constants $\tilde{L}_{nl}$ and $\tilde{L}_{pl}$ are defined in terms of the neutrino luminosities for a nonsuperfluid NS in beta equilibrium, as
$\displaystyle \tilde{L}_\alpha \equiv \frac{L_\alpha^{\rm eq}}{T_\infty^8} = \int_{{\rm core}} 4\pi r^2 {\rm e}^{\Lambda} S_\alpha(n) {\rm e}^{-6\Phi} {\rm d}r,$   (11)

where $\alpha$ indicates both the branch of the Murca process and the lepton involved in the reaction (Fernández & Reisenegger 2005). The term $S_\alpha$ is a slowly varying function of the baryon number density n (e.g., Yakovlev et al. 2001), and $\Lambda$ and $\phi$ are the usual Schwarzschild metric terms. The quantities $I_{M,\epsilon}^N$ and $I_{M,\Gamma}^N$are dimensionless phase-space integrals that contain the dependence of the emissivity and the net reaction rate, respectively, on the chemical imbalances $\eta_{npl}^\infty$ and on the energy gaps $\Delta _n$ and $\Delta_p$.

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

\begin{displaymath} x_j\equiv \frac{\epsilon_j-\mu_j}{kT} , \mbox{ }\mbox{ } x_... ...ox{ }\mbox{ and } \mbox{ } \xi_l\equiv \frac{\eta_{npl}}{kT} , \end{displaymath} (12)

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

\begin{displaymath} x_i\equiv\frac{v_{F_i}(p_i-p_{F_i})}{kT} \mbox{ and } z_i \equiv \mbox{sgn}(x_i)\sqrt{x_i^2+\delta_i^2}. \end{displaymath} (13)

In terms of these variables,
                                    $\displaystyle I_{M,\epsilon}^N$ = $\displaystyle \frac{60~480}{11~513\pi^8} \int_0^{\infty}dx_{\nu} x_{\nu}^3 \! \...  
    $\displaystyle \times ~ {\rm d}x_{p}{\rm d}x_{e}f(z_n)f(z_{N_i})f(z_{N_f}) f(z_p)f(x_e)$  
    $\displaystyle \times~ \left[\delta(x_{\nu}+\xi_l-z_n-z_{N_i}-z_{N_f}-z_p-x_e) \right.$  
    $\displaystyle +~ \left.\delta(x_{\nu}-\xi_l-z_n-z_{N_i}-z_{N_f}-z_p-x_e) \right]$ (14)

and
                                  $\displaystyle I_{M,\Gamma}^N$ = $\displaystyle \frac{60~480}{11~513\pi^8} \int_0^{\infty}{\rm d}x_{\nu} x_{\nu}^... ...}^{\infty}\! \int_{-\infty}^{\infty}\! {\rm d}x_{n}{\rm d}x_{N_i}{\rm d}x_{N_f}$  
    $\displaystyle \times ~ {\rm d}x_{p}{\rm d}x_{e}f(z_n)f(z_{N_i})f(z_{N_f}) f(z_p)f(x_e)$  
    $\displaystyle \times~ \left[\delta(x_{\nu}+\xi_l-z_n-z_{N_i}-z_{N_f}-z_p-x_e) \right.$  
    $\displaystyle -~ \left.\delta(x_{\nu}-\xi_l-z_n-z_{N_i}-z_{N_f}-z_p-x_e) \right],$ (15)

where $f(\cdot)$ is the Fermi function f(x)=1/(1+ex) 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., $\delta_n=\delta_p=0$), these integrals reduce to the polynomials calculated by Reisenegger (1995)

                                          $\displaystyle I_{M,\epsilon}^N (\delta_n=\delta_p=0)=F_M(\xi_l)=$  
    $\displaystyle \quad \quad 1 + \frac{22~020\xi_l^2}{11~513\pi^2} + \frac{5670\xi... (16)
    $\displaystyle I_{M,\Gamma}^N (\delta_n=\delta_p=0)=H_M(\xi_l)=$  
    $\displaystyle \quad \quad\frac{14~680\xi_l}{11~513\pi^2}+\frac{7560\xi_l^3}{11~513\pi^4}+\frac{840\xi_l^5}{11~513\pi^6} +\frac{24\xi_l^7}{11~513\pi^8},$ (17)

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

Since the z-variables defined in Eq. (13) depend on the energy gaps, the equality between $I_{M,\epsilon}^n$ and $I_{M,\epsilon}^p$, or between $I_{M,\Gamma}^n$ and $I_{M,\Gamma}^p$, 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

          $\displaystyle R_{M,\epsilon}^N(\xi_l,\delta_n,\delta_p)$ = $\displaystyle I_{M,\epsilon}^N(\xi_l,\delta_n,\delta_p)/F_M(\xi_l),$ (18)
$\displaystyle R_{M,\Gamma}^N(\xi_l,\delta_n,\delta_p)$ = $\displaystyle I_{M,\Gamma}^N(\xi_l,\delta_n,\delta_p)/H_M(\xi_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 $\eta_l$, $\delta_n$, and $\delta_p$, 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 $\eta_l$, $\delta_n$, $\delta_p$. In this way they cover a wide range of the infinite integration domain, for the parameters in Eq. (19) in the range $\delta_N \in [0,10^3]$ and $\xi_l\in[0,10^4]$.

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 $\eta\gtrsim\Delta$ and $\Delta>10kT$, the reduction factors become functions of $\Delta/\eta$ only. By using our numerical results, we show in Appendix B that this is a good approximation if $\Delta\gtrsim30kT$ and $\eta>\Delta$. 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)

                            $\displaystyle R_{M,\epsilon}(\Delta/\eta)$ = $\displaystyle \left( 28 (\Delta/\eta)^2\! +\! 105 (\Delta/\eta)^4 \! +\! \frac{105}{2}(\Delta/\eta)^6\! +\! \frac{35}{16}(\Delta/\eta)^8\right)$  
    $\displaystyle \times \ln\left(\frac{\Delta/\eta}{1+\sqrt{1-(\Delta/\eta)^2}}\right) + \sqrt{1-(\Delta/\eta)^2}$  
    $\displaystyle \times \left(1\! +\! \frac{551}{10}(\Delta/\eta)^2 \! \!+\!\! \frac{4327}{40}(\Delta/\eta)^4\! +\! \frac{1873}{80}(\Delta/\eta)^6\right)$ (20)
$\displaystyle R_{M,\Gamma}(\Delta/\eta)$ = $\displaystyle \left(21(\Delta/\eta)^2+ \frac{105}{2}(\Delta/\eta)^4 + \frac{105}{8}(\Delta/\eta)^6\right)$  
    $\displaystyle \times \ln\left(\frac{\Delta/\eta}{1+\sqrt{1-(\Delta/\eta)^2}}\right) + \sqrt{1-(\Delta/\eta)^2}$  
    $\displaystyle \times\left(1+\frac{759}{20}(\Delta/\eta)^2 \!\! +\!\! \frac{1779}{40}(\Delta/\eta)^4\! \! +\!\! \frac{16}{5}(\Delta/\eta)^6\right),$ (21)

where $\Delta$ 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 $\lim_{\Delta/\eta\rightarrow0}R=1$ and $\lim_{\Delta/\eta\rightarrow1}R=0$as expected. In addition, in the limit $kT\ll1$ the functions in Eqs. (16) and (17) tend to the highest power of the polynomial, i.e., $F(\xi)\cong\frac{9\xi^8}{11~513\pi^8}$ and $H(\xi)\cong\frac{24\xi^8}{11~513\pi^8}$. Thus, the integrals can be expressed as
                $\displaystyle I_{\epsilon}(\delta,\xi)$ = $\displaystyle \frac{9\xi^8}{11~513\pi^8}R_{\epsilon}(\delta/\xi), \mbox{and}$ (22)
$\displaystyle I_{\Gamma}(\delta,\xi)$ = $\displaystyle \frac{24\xi^8}{11~513\pi^8}R_{\Gamma}(\delta/\xi),$ (23)

where, doing some algebra, one can verify that these approximations satisfy the identity of Flores-Tulián & Reisenegger (2006), i.e., $\frac{{\rm d}I_{\epsilon}(\delta,\xi)}{{\rm d}\xi}=3I_{\Gamma}(\delta,\xi)$. In Appendix B, we compare these formulae with the numerical non-zero temperature calculations, finding that they agree to a very good approximation when $\Delta/T\gtrsim 30$.

\begin{figure} \par\includegraphics[width=9cm,height=5.5cm,clip]{13861fg01.eps} \end{figure} 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 $\log(R_{M,\Gamma}^n)$.

Open with DEXTER

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 $\Delta _{\rm thr}$ imposed by the energy gaps $\Delta _n$ and $\Delta_p$. If we consider the reaction $n+n_i\rightarrow n_f+p+e^-+\bar{\nu}_e$ at zero temperature, the neutron is energetically allowed to decay only when $2\mu_n-2\Delta_n>\mu_n+\Delta_n+\mu_p+\Delta_p+\mu_e$, which implies the simple condition $\eta_{npe}>3\Delta_n+\Delta_p$ (see Reisenegger 1997, for a schematicjustification of this). Finally, the threshold imposed by the gaps $\Delta _{\rm thr}$ for both branches of Murca reactions are

                             $\displaystyle \Delta_{\rm thr}$ = $\displaystyle 3\Delta_n+\Delta_p \quad \mbox{for the neutron branch and},$ (24)
$\displaystyle \Delta_{\rm thr}$ = $\displaystyle \Delta_n+3\Delta_p \quad \mbox{for 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 $\eta<\Delta_{\rm thr}=3\Delta_n+\Delta_p$. Remarkably, the contour lines of the reduction factor are almost straight lines that coincide with the slope of the contour levels of equation $3\Delta_n+\Delta_p$. This means that the reduction factor in this regime is close to a certain function of $\Delta _{\rm thr}$, which is valid for both branches and also for the reduction factor of the emissivity.

\begin{figure} \par\includegraphics[width=8.9cm,height=5.6cm,clip]{13861fg02.eps} \end{figure} Figure 2:

Ratio of the reduction factor of the emissivity $R_\epsilon $ to the reduction factor of the net reaction rate $R_\Gamma $ for the proton branch of Murca reactions with superfluid neutrons as a function of $\eta /\Delta _n$ for different values of $\Delta _n/T$. Ratio of the zero temperature approximations given by Eqs. (20) and (21) is given by the black dash-dotted 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 $R_\epsilon/R_\Gamma$ 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 $R_\epsilon/R_\Gamma$ approaches a value of one when $\eta/\Delta\rightarrow\infty$. 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 ${\sim}(\eta-\Delta_{\rm thr})$ that can be arbitrarily small as $\eta$ approaches  $\Delta _{\rm thr}$.

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 $\eta<\Delta_{\rm thr}$, we use the Gauss-Laguerre quadrature method explained in Appendix A, which happens for the early evolution of rotochemical heating (see Fig. 3). When $\eta>\Delta_{\rm thr}$, we use the analytical approximation of the integral presented in Appendix B.

2.4 Specific heat suppression

When T reaches $T_{\rm c}$, 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 $\sim$2.4 with respect to the normal-matter specific heat. Subsequently, when $T\ll T_{\rm c}$, 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 $C_{\rm V}$. These have been calculated by several authors for various temperature regimes, e.g. Pizzochero et al. (2002) for $T\ll T_{\rm c}$ and Maxwell (1979) for $0.2T_{\rm c}<T<T_{\rm c}$. The former set of authors obtain an exponential decaying control function of the form $\frac{\Delta}{kT}e^{-\Delta/kT}$ and the latter, which is the result we adopt, obtain the function

                           $\displaystyle C_{{\rm V},N}^{\rm sup}$ = $\displaystyle C_{{\rm V},N}\times3.15\frac{T_{{\rm c},N}}{T}e^{-1.76T_{{\rm c},N}/T}\times$  
    $\displaystyle \times \left[2.5-1.66 \frac{T}{T_{{\rm c},N}} +3.64 \left(\frac{T}{T_{{\rm c},N}}\right)^2 \right],$ (26)

where $T_{{\rm c},N}$ is the critical temperature of the nucleon N and $C_{{\rm V},N}$ 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 $C_{{\rm V},N}^{\rm sup}(T=0.2T_{\rm c})/C_{{\rm V},N}=0.005$, which is sufficiently small to ensure a the lepton contribution to the specific heat dominant. This formula also represents the discontinuous increase at $T=T_{\rm c}$, which implies that $C_{{\rm V},N}^{\rm sup}(T=T_{\rm c})/C_{{\rm V},N}=2.42$, 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 $C_{\rm V}=C_{{\rm V},e}+C_{{\rm V},\mu}$.

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 $T_{\rm c}$ and successively suppressed at lower temperatures because of an exponential decrease in the number of unpaired particles. Schematically, these neutrino reactions are

                          NN $\textstyle \rightarrow$ $\displaystyle N+N+\nu+\bar{\nu} \quad \mbox{ pair breaking (PB)},$ (27)
N+N $\textstyle \rightarrow$ $\displaystyle {NN}+\nu+\bar{\nu} \quad \mbox{ pair formation (PF)},$ (28)

where NN denotes the Cooper pair and N the excitation.
These authors find an emissivity $ Q_{\rm PB,PF}\sim 10^{28}(kT/\mbox{[MeV]})^7$ $\left[\mbox{erg}\mbox{ cm}^{-3}\mbox{ s}^{-1}\right]$ when $T \lesssim T_{\rm c}$. 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 $Q_{\rm MU}\sim 10^{29{-}30} (kT/\mbox{[MeV]})^8 \left[\mbox{erg}\mbox{ cm}^{-3}\mbox{ s}^{-1}\right]$. 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 $kT\sim10^{-3}$ MeV (see Fernández & Reisenegger 2005). In this limit, in which $T\ll T_{\rm c}$ (Flowers et al. 1976),
$\displaystyle Q_{\rm PB,PF}\sim 10^{28}\left(\frac{\Delta}{\mbox{MeV}}\right)^7...   (29)

which clearly becomes negligible when the values of the energy gaps are relatively large compared to the thermal scale, say $\Delta/kT\sim10^{2}{-}10^3$, because the exponential behavior rapidly suppresses the effects of these reactions. For $\Delta>0.01$ 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 + $\delta \upsilon $ + 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 $\rho_D=1.59\times10^{15}$ g cm-3, which corresponds to the central density of a 2 $M_\odot $ 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 $\Delta _n=0.1$ MeV by integrating numerically the emissivity and net reaction rate without any approximation.

\begin{figure} \includegraphics[width=9cm,height=5.5cm,clip]{13861fg03.eps}\par\includegraphics[width=9cm,height=5.5cm,clip]{13861fg04.eps} \end{figure} Figure 3:

Evolution of the internal temperature $T_\infty $, the surface temperature $T_{s,\infty }$ and the chemical imbalances $\eta _{npe}^\infty $, $\eta _{npu}^\infty $ for a star with the mass fixed to the PSR J0437-4715, i.e. $1.76~M_{\odot}$ (Verbiest et al. 2008), built with the A18 + $\delta \upsilon $ + UIX* EOS, with the initial condition $T_{\infty }=10^9$ K and null chemical imbalances. The spin-down is assumed to be due to the magnetic dipole radiation, with dipolar field strength $B=2.8 \times 10^8$ G and initial period of P0=1 ms. The error bar is the $90\%$ confidence level for the surface temperature measured for the PSR J0437-4715 (Kargaltsev et al. 2004) at the current spin-down parameters. Upper panel: nonsuperfluid case (null energy gaps). Lower panel: superfluidity of neutrons with $\Delta _n=0.1$ 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 105 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; 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 $\Delta _n$. 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 $\Delta_n=$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 $\Delta _{\rm thr}$, that for Fig. 3 is $\Delta_{\rm thr}=\Delta_n$ (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 $\Delta _n$ and $\Delta_p$ in the range of $\Delta_{\rm thr}=0.05{-}1$ MeV, we verify that for the Murca processes, the zero temperature approach of Reisenegger (1997) is valid, because in the quasi-steady state $\eta_{npe},\eta_{np\mu}\gtrsim\Delta_{\rm thr}$. 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 Quasi-steady state

Sufficiently old stars will have reached the quasi-steady state, where $\dot{T}_\infty=0$, $\dot{\eta}_{npe}^{\infty}=0$, and $\dot{\eta}_{np\mu}^{\infty}=0$. In this state, the Eqs. (1), (3), and (4) can be written as

                         $\displaystyle L_\gamma^{\infty,{\rm qs}}$ = $\displaystyle \eta_{npe}^{\infty,qs}\Delta\tilde{\Gamma}_{npe}+\eta_{np\mu}^{\infty,qs}\Delta\tilde{\Gamma}_{np\mu}- L_\nu^\infty,$ (30)
$\displaystyle 2 W_{npe}\Omega \dot{\Omega}$ = $\displaystyle Z_{npe}\Delta\tilde{\Gamma}_{npe}+ Z_{np}\Delta\tilde{\Gamma}_{np\mu},$ (31)
$\displaystyle 2 W_{np\mu}\Omega \dot{\Omega}$ = $\displaystyle Z_{np}\Delta\tilde{\Gamma}_{npe}+ Z_{np\mu}\Delta\tilde{\Gamma}_{np\mu},$ (32)

where the superscript qs 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 $\Omega \dot{\Omega}$. Figure 4 shows the solution to this equations when only the neutrons are superfluid, such that $\Delta_{\rm thr}=\Delta_n$. This illustrates that the chemical imbalance always stabilizes at a value slightly higher than $\Delta _{\rm thr}$,
$\displaystyle \eta_{npe}^{\infty,{\rm qs}}\approx\eta_{np\mu}^{\infty,{\rm qs}}... ...elta_{\rm thr}= \mbox{min}\left( \Delta_n+3\Delta_p,3\Delta_n+\Delta_p \right).$     (33)

We can also observe from Fig. 4 that the solution becomes closer to $\Delta _{\rm thr}$ 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 $\sim\eta^7$. Hence, the reduction factors need to block more reactions, which happens for values of the chemical imbalance closer to $\Delta _{\rm thr}$. 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.

\begin{figure} \par\includegraphics[width=9cm,height=5.5cm,clip]{13861fg05.eps} \vspace*{2mm}\end{figure} Figure 4:

Chemical imbalance of the electrons in the quasi-steady state with only neutrons as superfluid particles as a function of $\Omega \dot{\Omega}$for three values of the energy gap $\Delta _n$ (solid lines) and the nonsuperfluid case (dashed line) for a $1.76 M_{\odot }$ star calculated with the A18 + $\delta \upsilon $ + 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 $\Delta\tilde{\Gamma}_{npe}$ and $\Delta\tilde{\Gamma}_{np\mu}$are determined entirely by Eqs. (31) and (32). Solving these equations, we obtain

                         $\displaystyle \Delta\tilde{\Gamma}_{npe}$ = $\displaystyle 2 \Omega \dot{\Omega} \frac{Z_{np\mu} W_{npe}-Z_{np} W_{np\mu}} {Z_{npe}Z_{np\mu}-Z_{np}^2}\equiv -K_{e} \Omega \dot{\Omega},$ (34)
$\displaystyle \Delta\tilde{\Gamma}_{np\mu}$ = $\displaystyle 2 \Omega \dot{\Omega} \frac{Z_{npe} W_{np\mu}-Z_{np} W_{npe}} {Z_{npe}Z_{np\mu}-Z_{np}^2}\equiv -K_{\mu} \Omega \dot{\Omega},$ (35)

where the constants Ke and $K_\mu $ 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 $\xi_l\gtrsim10^2$, the polynomials in Eqs. (16) and (17) are in a very good approximation $F(\xi_l)=\frac{9\xi_l^8}{11~513\pi^8}$ and $H(\xi_l)=\frac{24\xi_l^7}{11~513\pi^8}$, 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

                                   $\displaystyle L_\gamma^{\infty,{\rm qs}}$ = $\displaystyle \left[ \eta_{npe}^{\infty,{\rm qs}} K_e\left(1-\frac{9}{24} \frac... ...L}_{pe} R_{M,\Gamma}^p\left(\eta_{npe}^{\infty,{\rm qs}}\right)}\right) \right.$  
    $\displaystyle + \left. \eta_{np\mu}^{\infty,{\rm qs}}K_\mu\left(1-\frac{9}{24} ... ...{p\mu} R_{M,\Gamma}^p\left(\eta_{np\mu}^{\infty,{\rm qs}}\right)}\right)\right]$  
    $\displaystyle \times \left\vert\Omega \dot{\Omega}\right\vert,$ (36)

where we have just rearranged the quasi-steady 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 $\eta_{npe}^{\infty,{\rm qs}}\approx\eta_{np\mu}^{\infty,{\rm qs}}\approx\Delta_{\rm thr}^\infty$. 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 $\Delta_{\rm thr}>0.05$ MeV.

Our second approximation is to neglect the ratios of the reduction factors $R_{M,\epsilon}^N/R_{M,\Gamma}^N$. As we discussed in Sect. 2.3.1, the emissivity is more suppressed by the gaps than the net reaction rate when $\eta>\Delta$. It becomes an acceptable limit when the values of the gaps are relatively large compared to the thermal energy, say $\Delta/kT\gtrsim10^2$, which is the limit of interest for rotochemical heating. We checked that when we consider $\Delta_{\rm thr}\sim0.05{-}1$ MeV, the ratio is $R_{M,\epsilon}^N/R_{M,\Gamma}^N\lesssim0.1$.

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

$\displaystyle L_\gamma^{\infty,{\rm qs}} \approx\left( K_e+K_\mu \right)\Delta_{\rm thr}\left\vert\Omega \dot{\Omega}\right\vert,$     (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 $\Delta _{\rm thr}$ to reheat the star, as in the description of Reisenegger (1997). The approximation is even more accurate for the surface temperature, since $T_{s,\infty}^{{\rm qs}}\propto \left(L_\gamma^{\infty,{\rm qs}}\right)^{1/4}$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 $L_\gamma^{\infty,{\rm qs}}$ and $T_{s,\infty}^{{\rm qs}}$. 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 $T_{s,\infty}^{{\rm qs}}\propto \vert\Omega \dot{\Omega}\vert^{1/4}$, while for the nonsuperfluid case, Fernández & Reisenegger (2005) obtain $T_{s,\infty}^{{\rm qs}}\propto \vert\Omega \dot{\Omega}\vert^{2/7}$.

\begin{figure} \par\includegraphics[width=9cm,height=6cm,clip]{13861fg06.eps}\vspace*{2mm} \end{figure} Figure 5:

Surface temperature predicted for the quasi-steady state for several values of $\Delta_{\rm thr}=\min \left(\Delta_n+3\Delta_p,3\Delta_n+\Delta_p \right)$(solid lines) and the nonsuperfluid case (dash-dotted line) as a function of the spin-down rate $\Omega \dot{\Omega}$ with the same stellar parameters as in Fig. 4.

Open with DEXTER

\begin{figure} \par\includegraphics[width=9cm,height=6cm,clip]{13861fg07.eps} \end{figure} Figure 6:

Values of Ke (solid lines) and $K_\mu $(dash-dotted lines) for different EOSs for a range of masses from 1 $M_\odot $ to the mass at which direct Urca with electrons is opened for each EOS.

Open with DEXTER

In Fig. 6, we plot Ke and $K_\mu $ 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 + $\delta \upsilon $ and A18 + $\delta \upsilon $ + 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_\odot $. For this set of EOSs, we show in Fig. 6 that $K_e+K_\mu\simeq(2{-}6)\times10^{47}\left[\mbox{s}^2\right]$, which, using Eq. (37), infers a bolometric luminosity of

$\displaystyle L_\gamma^{\infty,{\rm qs}} \simeq \left(1{-}4\right) \times 10^{3... ...ht) \left(\frac{\dot{P}_{-20}}{P_{{\rm ms}}^3}\right)\mbox{ erg}\mbox{ s}^{-1},$     (38)

where $\dot{P}_{-20}$ is the period derivative in units of 10-20 and $P_{{\rm ms}}$ is the period in milliseconds. Finally, we can express the surface temperature using the expression for the luminosity given by Eq. (37) as
$\displaystyle T_{s,\infty}^{{\rm qs}} \simeq \left(\frac{ K_e+K_\mu}{4\pi\sigma... ...ght)^{1/4}\Delta_{\rm thr}^{1/4}\left\vert\Omega \dot{\Omega}\right\vert^{1/4},$     (39)

where $\sigma$ is the Stefan-Boltzmann constant. For the set of EOSs that we use, we obtain
$\displaystyle T_{s,\infty}^{{\rm qs}} \simeq \left(5.7{-}6.6\right) \times 10^{...   (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 $\Delta _{\rm thr}$, using one particular EOS. In this case, we obtain the allowed range

$\displaystyle 0.05~[\mbox{MeV}]\lesssim\mbox{min}\left( \Delta_n+3\Delta_p,3\Delta_n+\Delta_p \right)\lesssim0.2~[\mbox{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
$\displaystyle 0.05~[\mbox{MeV}]\lesssim\mbox{min}\left( \Delta_n+3\Delta_p,3\Delta_n+\Delta_p \right)\lesssim0.45~[\mbox{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).

\begin{figure} \par\includegraphics[width=9cm,clip]{13861fg08.eps}\vspace*{2mm} \end{figure} Figure 7:

Surface temperature predicted in the quasi-steady state as a function of the the stellar radius for a star built with the A18 + $\delta \upsilon $ + UIX* EOS for several values of $\Delta _{\rm thr}$ in units of MeV (solid lines) and the current spin-down rate of the MSP J0437. Vertical dotted lines indicates the mass range $1.76\pm0.2~M_\odot$ allowed for MSP J0437 (Verbiest et al. 2008). Dashed-dotted lines are the 90$\%$ confidence contours of the black-body 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 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 $\Omega \dot{\Omega}$ 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

$\displaystyle \dot{\eta}^\infty_{npl} = 2W_{npl}\Omega \dot{\Omega},$     (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 $\eta^\infty_{npl}=\Delta_{\rm thr}$, 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 Pi 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
$\displaystyle P_i < \left(\frac{1}{P^2} + \frac{\Delta_{\rm thr}}{4 \pi^2\left\... ...}{4} \left(\frac{\Delta_{\rm thr}}{\mbox{MeV}}\right)\right)^{-1/2} \mbox{ ms},$     (44)

where $P_{{\rm ms}}$ 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., $\Delta_{\rm thr}=0.45$ MeV, should be $P_i\lesssim2.7$ 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 $\Delta_{\rm thr}=\mbox{min}\left( \Delta_n+3\Delta_p,3\Delta_n+\Delta_p \right)$, 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

$\displaystyle L_\gamma^{\infty,{\rm qs}} \simeq \left(1{-}4\right) \times 10^{3...     (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
$\displaystyle 0.05 [\mbox{MeV}]\lesssim\Delta_{\rm thr}\lesssim0.45~[\mbox{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 (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 PFB-06/2007.

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 $[0, \infty]$, 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 nth-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 $\int_0^{\infty}f(x){\rm d}x$ is that $f(x)\cdot W(x)$ has to be a smooth function, such as a polynomial (Davis & Rabinowitz 1975).

\begin{figure} \par\includegraphics[width=12cm,clip]{13861fg09.eps}\vspace*{-1.5mm} \end{figure} 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 $\xi =\delta _n=$ $\delta _p=10$.

Open with DEXTER

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

                                   $\displaystyle I_{D,\Gamma}$ = $\displaystyle \! \! \int_0^{\infty}\! {\rm d}x_{\nu} x_{\nu}^2 \int_{0}^{\infty... ...{j_n=\pm 1} \sum_{j_p=\pm 1} f\left(j_nz_n\right) f\left(j_p z_p\right) \right.$  
    $\displaystyle \times \left. \left[ f\left(x_{\nu}\! +\! \xi-j_n z_n\! -\! j_p z... ...! f\left(x_{\nu}\! -\! \xi\! -\! j_n z_n\! -\! j_p z_p\right) \right] \right\},$ (A.1)

where $z_i=\sqrt{x_i^2+\delta_i^2}$ (i=n,p) and $f(\cdot)$ 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 $f(\pm\sqrt{x_n^2+\delta_n^2})\cdot f(\alpha \mp \sqrt{x_n^2+\delta_n^2})$, where $\alpha$ 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 $x_n \gg \delta_n$, which is the relevant scale for this variable. Thus, in this limit $f(-\sqrt{x_n^2+\delta_n^2})\sim e^{-x_n}$ and $f(\sqrt{x_n^2+\delta_n^2})\sim 1$, while the other Fermi functions have the opposite behaviors $f(\alpha - \sqrt{x_n^2+\delta_n^2})\sim 1$ and $f(\alpha + \sqrt{x_n^2+\delta_n^2})\sim e^{-x_n}$, 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 $\alpha$ take. In this sense, higher values of $\alpha$, which depends on the inputs $\delta_n$, $\delta_n$, and $\xi$, 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 $f\left(x_{\nu}\pm \xi-j_n z_n-j_p z_p\right)$ in the integrand of Eq. (A.1), as in Villain & Haensel (2005). For the neutrino variable, this scale is $S_{x_\nu}\equiv \xi+\delta_n+\delta_p$, and for the superfluid nucleons they are $S_{x_n}\equiv \xi+\sqrt{2\delta_n}+\delta_p$ and $S_{x_p}\equiv \xi+\sqrt{2\delta_p}+\delta_n$. The numerical formula is, finally

                               $\displaystyle I_{D,\Gamma}$ = $\displaystyle \sum_{i_{\nu}=1}^{n_{\nu}} \sum_{i_n=1}^{n_n} \sum_{i_p=1}^{n_p} ... ..._n}S_{x_p} x_{i_\nu}^2\cdot f\left(j_n z_{i_n}\right) f\left(j_p z_{i_p}\right)$  
    $\displaystyle \times \left[f\left(x_{i_\nu}\! +\! \xi\! -\! j_n z_n\! -\! j_p z... ...! f\left(x_{i_\nu}\! -\! \xi\! -\! j_n z_{i_n}\! -\! j_p z_{i_p}\right)\right],$ (A.2)

where $z_{i_{n}}=\sqrt{x_{i_{n}}^2+\delta_{n}^2}$ and $z_{i_{p}}=\sqrt{x_{i_{p}}^2+\delta_{p}^2}$. The values $x_{i_{\nu}}$, xin, and xip are the roots of the Laguerre polynomials of order $n_\nu$, nn, and np, multiplied by $S_{x_\nu}$, Sxn, and Sxp, respectively . The values $W_{i_{\nu}}$, Win, and Wip 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 Legendre-like 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 ( $\xi=\frac{\eta}{kT}\gtrsim100$). 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

\begin{displaymath}\xi \sim \delta_{n,p}=\frac{\Delta_{n,p}}{kT}\gg1. \end{displaymath} (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 phase-space integral for the net reaction rate is
                               $\displaystyle I_{D,\Gamma}$ = $\displaystyle \int_0^{\infty}{\rm d}x_{\nu} x_{\nu}^2 \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} {\rm d}x_{n}{\rm d}x_{p}f(z_n)f(z_p)$  
    $\displaystyle \times \{f(x_{\nu}+\xi-z_n-z_p) - f(x_{\nu}-\xi-z_n-z_p) \},$ (B.2)

where $z_{i}=\mbox{sgn}(x_{i})\sqrt{x_{i}^2+\delta_{i}^2}$, 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 $\log(\Delta/T)\geq 1$, and as soon as $\eta>\Delta$, 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., $f(x)\longrightarrow \Theta(-x)$, the integral gives

                               $\displaystyle I_{D,\Gamma}$ = $\displaystyle \int_0^{\infty}{\rm d}x_{\nu} x_{\nu}^2 \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} {\rm d}x_{n}{\rm d}x_{p}\Theta(-z_n)\Theta(-z_p)$  
    $\displaystyle \times~ \{\Theta(-x_{\nu}-\xi+z_n+z_p) - \Theta(-x_{\nu}+\xi+z_n+z_p) \}.$ (B.3)

Rearranging the limits of integration and eliminating the $\mbox{sgn}(\cdot)$ function from $z_{i}=\mbox{sgn}(x_{i})\sqrt{x_{i}^2+\delta_{i}^2}$, the integral can be reduced to
                           $\displaystyle I_{D,\Gamma}$ = $\displaystyle \int_0^{\infty}\int_0^{\infty}\int_0^{\infty}{\rm d}x_{\nu} {\rm d}x_{n}{\rm d}x_{p} x_{\nu}^2$  
    $\displaystyle \times~ \Theta\left(\xi-x_{\nu}-\sqrt{x_n^2+\delta_n^2}-\sqrt{x_p^2+\delta_p^2}\right).$ (B.4)

By applying these approximations, the reactions can be seen to be open only when $\delta_{n}+\delta_{p}<\xi$. (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 $\xi, \delta_n, \delta_p$, 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

$\displaystyle u_{\nu}\equiv \frac{x_{\nu}}{\xi} \mbox{ } \mbox{ and } \mbox{ } u_{i}\equiv \sqrt{ \frac{x_{i} ^2+\delta_{i}^2}{\xi^2}} \mbox{, } (i=n,p) \mbox{,}$     (B.5)

where the Jacobian of this transformation is $J(u_{\nu},u_n,u_p)=\xi^3\frac{u_nu_p}{\sqrt{u_n^2-(\Delta_n/\eta)^2}\sqrt{u_p^2-(\Delta_p/\eta)^2} }$, considering that $\Delta_p/\eta=\delta_p/\xi$, 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
                           $\displaystyle I_{D,\Gamma}$ = $\displaystyle \xi^5\int_{\Delta_n/\eta}^{1}\int_{\Delta_p/\eta}^{1-u_n}\int_0^{1-u_n-u_p} {\rm d}u_{n} {\rm d}u_{p} {\rm d}u_{\nu} u_{\nu}^2$  
    $\displaystyle \times \frac{u_n u_p}{\sqrt{u_n^2-(\Delta_n/\eta)^2}\sqrt{u_p^2-(\Delta_p/\eta)^2} } \cdot$ (B.6)

Integrating over the neutrino variable, we obtain
                          $\displaystyle I_{D,\Gamma}$ = $\displaystyle \frac{\xi^5}{3}\int_{\Delta_n/\eta}^{1}\int_{\Delta_p/\eta}^{1-u_n} {\rm d}u_{n} {\rm d}u_{p} (1-u_n-u_p)^3$  
    $\displaystyle \times \frac{u_n u_p}{\sqrt{u_n^2-(\Delta_n/\eta)^2}\sqrt{u_p^2-(\Delta_p/\eta)^2} } \cdot$ (B.7)

The integration over the proton variable yields
$\displaystyle I_{D,\Gamma} = \frac{\xi^5}{3}\int_{\Delta_n/\eta}^{1-\Delta_p/\e... ...}u_{n}\frac{u_n}{\sqrt{u_n^2-(\Delta_n/\eta)^2}} K_{\Gamma}(u_n,\Delta_p/\eta),$     (B.8)

where
                  $\displaystyle K_{\Gamma}(u_n,\Delta_p/\eta)$ = $\displaystyle -\frac{3(\Delta_p/\eta)^2}{2}\left[\frac{(\Delta_p/\eta)^2}{4} +(1-u_n)^2 \right]$  
    $\displaystyle \times \ln\left[\frac{1-u_n+\sqrt{(1-u_n)^2-(\Delta_p/\eta)^2}}{\Delta_p/\eta}\right]$  
    $\displaystyle +\sqrt{(1-u_n)^2-(\Delta_p/\eta)^2}$  
    $\displaystyle \times \left[\frac{13(\Delta_p/\eta)^2}{8}\left(1-u_n \right) + \frac{\left(1-u_n \right)^3}{4} \right]\cdot$ (B.9)

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 $x_\nu^2$ with $x_\nu^3$ in Eq. (B.4). The formula obtained in this case is
                                $\displaystyle I_{D,\epsilon}$ = $\displaystyle \frac{\xi^6}{4}\int_{\Delta_n/\eta}^{1-\Delta_p/\eta} {\rm d}u_{n}\frac{u_n} {\sqrt{u_n^2-(\Delta_n/\eta)^2}} K_{\epsilon}(u_n,\Delta_p/\eta),$ (B.10)

where
                         $\displaystyle K_{\epsilon}(u_n,\Delta_p/\eta)$ = $\displaystyle 2(\Delta_p/\eta)^2\left[\frac{3}{4}(\Delta_p/\eta)^2(1-u_n)+(1-u_n)^3\right]$  
    $\displaystyle \times \ln\left[\frac{1-u_n+\sqrt{(1-u_n)^2-(\Delta_p/\eta)^2}}{\Delta_p/\eta}\right]$  
    $\displaystyle + \frac{\sqrt{(1-u_n)^2-(\Delta_p/\eta)^2}}{5} \Bigg\{ (1-u_n)^4+\frac{(\Delta_p/\eta)^2}{6}$  
    $\displaystyle \times \left[ 83(1-u_n)^2+16 (\Delta_p/\eta)^2\right] \Bigg\}.$ (B.11)

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 $I_{M,\Gamma}^{p}(\delta_n,\delta_p)=I_{M,\Gamma}^{n}(\delta_p,\delta_n)$ only one expression is needed. Thus, without loss of generality, we only write the neutron branch expression of the net reaction rate of Murca
                                    $\displaystyle I_{M,\Gamma}^{n}$ = $\displaystyle \frac{\xi^7}{3}\int_{\Delta_n/\eta}^{1}\int_{\Delta_n/\eta}^{1-u_n} \int_{\Delta_n/\eta}^{1-u_n-u_{n_i}}{\rm d}u_{n} {\rm d}u_{n_i} {\rm d}u_{n_f}$  
    $\displaystyle \times ~ \frac{u_n u_{n_i} u_{n_f}} {\sqrt{u_n^2-(\Delta_n/\eta)^2}\sqrt{u_{n_i}^2-(\Delta_n/\eta)^2}\sqrt{u_{n_f}^2-(\Delta_n/\eta)^2}}$  
    $\displaystyle \times~ K_{\Gamma}(u_n+u_{n_i}+u_{n_f},\Delta_p/\eta),$ (B.12)

where $K_{\Gamma}(\cdot,\cdot)$ is the function given by Eq. (B.9) and ni and nf represent the initial and final neutron spectator, respectively. In the same way, the emissivity has an expression in terms of the function $K_{\epsilon}(\cdot,\cdot)$defined by Eq. (B.11):
                                      $\displaystyle I_{M,\epsilon}^n$ = $\displaystyle \frac{\xi^8}{4}\int_{\Delta_n/\eta}^{1}\int_{\Delta_n/\eta}^{1-u_n}\int_{\Delta_n/\eta}^{1-u_n-u_{n_i}}{\rm d}u_{n} {\rm d}u_{n_i} {\rm d}u_{n_f}$  
    $\displaystyle \times ~ \frac{u_n u_{n_i} u_{n_f}} {\sqrt{u_n^2-(\Delta_n/\eta)^2}\sqrt{u_{n_i}^2-(\Delta_n/\eta)^2}\sqrt{u_{n_f}^2-(\Delta_n/\eta)^2}}$  
    $\displaystyle \times~ K_{\epsilon}(u_n+u_{n_i}+u_{n_f},\Delta_p/\eta).$ (B.13)

In both cases, the net reaction rate and the emissivity, the condition to recover the non-superfluid 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 tanh-sinh method (Takahasi & Mori 1973), which integrates the singularities in $\Delta_n/\eta$ and/or in $\Delta_p/\eta$ quite rapidly.

\begin{figure} \par\includegraphics[width=10cm, height=5cm,clip]{13861fg10.eps} \end{figure} Figure B.1:

Reduction factor R for the net reaction rate of direct Urca with one superfluid particle as a function of the ratio $\eta /\Delta $ for several values of the ratio $\Delta /T$ (dashed lines) using the method detailed in Appendix A, and the limit case $T/\Delta \rightarrow 0$ 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 $\Delta$. The reduction factor R (for all the previous cases) obtained from this simplification is given by the integral


 		R = $\displaystyle (\alpha+1) \int_{\Delta/\eta}^{1} {\rm d}u\frac{u (1-u)^\alpha}{\sqrt{u^2-(\Delta/\eta)^2}} \mbox{ }\mbox{ ,}$
(B.14)

where $\alpha=4$ and $\alpha=5$ for Durca processes (net-reaction rate and emissivity, respectively), and $\alpha=6$ and $\alpha=7$ for Murca processes (net-reaction rate and emissivity, respectively). The solution to this integral is
$\displaystyle R(\Delta_n/\eta)$ = $\displaystyle P(\Delta/\eta)\ln\left(\frac{1}{1+\sqrt{1-(\Delta/\eta)^2}}\right)$  
    $\displaystyle +~ Q(\Delta/\eta)\sqrt{1-(\Delta/\eta)^2},$ (B.15)

where $P(\cdot)$ and $Q(\cdot)$ are the polynomials
                                       P(x) = $\displaystyle 10x^2+ \frac{15}{2}x^4 \mbox{ } \mbox{ } \mbox{ and }\mbox{ } \mbox{ } Q(x)=1+\frac{83}{6}x^2 +\frac{8}{3}x^4$  
    $\displaystyle \mbox{for the Durca net reaction rate,}$ (B.16)
P(x) = $\displaystyle 15x^2+ \frac{45}{2}x^4+ \frac{15}{8}x^6\mbox{ } \mbox{ } \mbox{ and }\mbox{ } \mbox{ } Q(x)=1+\frac{97}{4}x^2 +\frac{113}{8}x^4$  
    $\displaystyle \mbox{for the Durca emissivity,}$ (B.17)
P(x) = $\displaystyle 21x^2+ \frac{105}{2}x^4 + \frac{105}{8}x^6$  
    $\displaystyle \mbox{ and }\mbox{ } \mbox{ } Q(x)=1+\frac{759}{20}x^2 +\frac{1779}{40}x^4 +\frac{16}{5}x^6$  
    $\displaystyle \mbox{for the Murca net reaction rate,}$ (B.18)
P(x) = $\displaystyle 28 x^2+ 105 x^4 + \frac{105}{2}x^6 + \frac{35}{16}x^8$  
    $\displaystyle \mbox{ and }\mbox{ } \mbox{ } Q(x)=1+\frac{551}{10}x^2 +\frac{4327}{40}x^4+ \frac{1873}{80}x^6$  
    $\displaystyle \mbox{for the Murca emissivity.}$ (B.19)

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 $\Delta/T\gtrsim 30$.


References

  1. Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (New York: Dover), 890 [Google Scholar]
  2. Amundsen, L., & Østgaard, E. 1985, Nucl. Phys. A, 442, 163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804 [NASA ADS] [CrossRef] [Google Scholar]
  4. Davis, P., & Rabinowitz P. 1975, Methods of numerical integration, (New York: Academic Press) [Google Scholar]
  5. Deller, A., Verbiest, J., Tingay, S., & Bailes, M., 2008, A&A, 685, 67 [Google Scholar]
  6. Fernández, R., & Reisenegger, A. 2005, A&A, 625, 291 [Google Scholar]
  7. Flores-Tulián, S., & Reisenegger, A. 2006, MNRAS, 372, 276 [NASA ADS] [CrossRef] [Google Scholar]
  8. Flowers, E., Ruderman, M., & Sutherland, P. 1976, A&A, 205, 541 [Google Scholar]
  9. Haensel, P., & Pichon, B. 1994, A&A, 283, 313 [NASA ADS] [Google Scholar]
  10. Hansen, B. M. S., & Phinney, E. S. 1998, MNRAS, 294, 569 [NASA ADS] [CrossRef] [Google Scholar]
  11. Kargaltsev, O., Pavlov, G. G., & Romani, R. 2004, A&A, 602, 327 [Google Scholar]
  12. Levenfish, K. P., & Yakovlev, D. G. 1994, Astron. Rep., 38, 247 [NASA ADS] [Google Scholar]
  13. Lombardo, U., & Schulze, H. 2001, Lect. Not. Phys., 578, 30 [Google Scholar]
  14. Maxwell, O. V. 1979, ApJ, 231, 201 [NASA ADS] [CrossRef] [Google Scholar]
  15. Page, D., Lattimer, J. M., Prakash, M., & Steiner, A. W. 2009, ApJ, 707, 1131 [NASA ADS] [CrossRef] [Google Scholar]
  16. Pethick, C. J., Ravenhall, D. G., & Lorentz, C. P. 1995, Nucl. Phys. A, 584, 675 [NASA ADS] [CrossRef] [Google Scholar]
  17. Pizzochero, P. M. 2002, A&A, 569, 381 [Google Scholar]
  18. Potekhin, A. Y., Chabrier, G., & Yakovlev, D. G. 1997, A&A, 323, 415 [NASA ADS] [Google Scholar]
  19. Prakash, M., Ainsworth, T. L., & Lattimer, J. M. 1988, , 61, 2518 [Google Scholar]
  20. Reisenegger, A. 1995, A&A, 442, 749 [Google Scholar]
  21. Reisenegger, A. 1997, A&A, 485, 313 [Google Scholar]
  22. Reisenegger, A., Jofré, P., Fernández, R., & Kantor, E. 2006, A&A, 653, 568 [Google Scholar]
  23. Takahasi, H., & Mori, M. 1973, PRIMS, 9, 721 [Google Scholar]
  24. Thorne, K. S. 1977, ApJ, 212, 825 [NASA ADS] [CrossRef] [Google Scholar]
  25. Villain, L., & Haensel, P. 2005, A&A, 444, 539 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Verbiest, J. P. W., Bailes, M., van Straten, W., et al. 2008, A&A, 679, 675 [Google Scholar]
  27. Yakovlev, D. G. 2001, Phys. Rep., 354, 1 [NASA ADS] [CrossRef] [Google Scholar]
  28. Yakovlev, D. G., & Pethick, C. J. 2004, ARA&A, 42, 169 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]

Footnotes

... Petrovich[*]
Present address: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

All Figures

  \begin{figure} \par\includegraphics[width=9cm,height=5.5cm,clip]{13861fg01.eps} \end{figure} 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 $\log(R_{M,\Gamma}^n)$.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.9cm,height=5.6cm,clip]{13861fg02.eps} \end{figure} Figure 2:

Ratio of the reduction factor of the emissivity $R_\epsilon $ to the reduction factor of the net reaction rate $R_\Gamma $ for the proton branch of Murca reactions with superfluid neutrons as a function of $\eta /\Delta _n$ for different values of $\Delta _n/T$. Ratio of the zero temperature approximations given by Eqs. (20) and (21) is given by the black dash-dotted line.

Open with DEXTER
In the text

  \begin{figure} \includegraphics[width=9cm,height=5.5cm,clip]{13861fg03.eps}\par\includegraphics[width=9cm,height=5.5cm,clip]{13861fg04.eps} \end{figure} Figure 3:

Evolution of the internal temperature $T_\infty $, the surface temperature $T_{s,\infty }$ and the chemical imbalances $\eta _{npe}^\infty $, $\eta _{npu}^\infty $ for a star with the mass fixed to the PSR J0437-4715, i.e. $1.76~M_{\odot}$ (Verbiest et al. 2008), built with the A18 + $\delta \upsilon $ + UIX* EOS, with the initial condition $T_{\infty }=10^9$ K and null chemical imbalances. The spin-down is assumed to be due to the magnetic dipole radiation, with dipolar field strength $B=2.8 \times 10^8$ G and initial period of P0=1 ms. The error bar is the $90\%$ confidence level for the surface temperature measured for the PSR J0437-4715 (Kargaltsev et al. 2004) at the current spin-down parameters. Upper panel: nonsuperfluid case (null energy gaps). Lower panel: superfluidity of neutrons with $\Delta _n=0.1$ MeV (dashed line).

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,height=5.5cm,clip]{13861fg05.eps} \vspace*{2mm}\end{figure} Figure 4:

Chemical imbalance of the electrons in the quasi-steady state with only neutrons as superfluid particles as a function of $\Omega \dot{\Omega}$for three values of the energy gap $\Delta _n$ (solid lines) and the nonsuperfluid case (dashed line) for a $1.76 M_{\odot }$ star calculated with the A18 + $\delta \upsilon $ + 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

  \begin{figure} \par\includegraphics[width=9cm,height=6cm,clip]{13861fg06.eps}\vspace*{2mm} \end{figure} Figure 5:

Surface temperature predicted for the quasi-steady state for several values of $\Delta_{\rm thr}=\min \left(\Delta_n+3\Delta_p,3\Delta_n+\Delta_p \right)$(solid lines) and the nonsuperfluid case (dash-dotted line) as a function of the spin-down rate $\Omega \dot{\Omega}$ with the same stellar parameters as in Fig. 4.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,height=6cm,clip]{13861fg07.eps} \end{figure} Figure 6:

Values of Ke (solid lines) and $K_\mu $(dash-dotted lines) for different EOSs for a range of masses from 1 $M_\odot $ to the mass at which direct Urca with electrons is opened for each EOS.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm,clip]{13861fg08.eps}\vspace*{2mm} \end{figure} Figure 7:

Surface temperature predicted in the quasi-steady state as a function of the the stellar radius for a star built with the A18 + $\delta \upsilon $ + UIX* EOS for several values of $\Delta _{\rm thr}$ in units of MeV (solid lines) and the current spin-down rate of the MSP J0437. Vertical dotted lines indicates the mass range $1.76\pm0.2~M_\odot$ allowed for MSP J0437 (Verbiest et al. 2008). Dashed-dotted lines are the 90$\%$ confidence contours of the black-body 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

  \begin{figure} \par\includegraphics[width=12cm,clip]{13861fg09.eps}\vspace*{-1.5mm} \end{figure} 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 $\xi =\delta _n=$ $\delta _p=10$.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=10cm, height=5cm,clip]{13861fg10.eps} \end{figure} Figure B.1:

Reduction factor R for the net reaction rate of direct Urca with one superfluid particle as a function of the ratio $\eta /\Delta $ for several values of the ratio $\Delta /T$ (dashed lines) using the method detailed in Appendix A, and the limit case $T/\Delta \rightarrow 0$ calculated from Eqs. (B.15) and (B.16) (solid line).

Open with DEXTER
In the text

Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.