Longperiod thermal oscillations in superfluid millisecond pulsars
Departamento de Astronomía y AstrofísicaPontificia Universidad Católica de Chile, Casilla 306 Santiago 22, Chile
email: cpetrovi@astro.princeton.edu
Received: 18 August 2010
Accepted: 29 January 2011
Context. In previous papers, we have shown that, as the rotation of a neutron star slows down, it will be internally heated as a consequence of the progressively changing mix of particles (rotochemical heating). In previously studied cases (nonsuperfluid neutron stars or superfluid stars with only modified Urca reactions), this leads to a quasisteady state in which the star radiates thermal photons for a long time, possibly accounting for the ultraviolet radiation observed from the millisecond pulsar J04374715.
Aims. For the first time, we explore the phenomenology of rotochemical heating with direct Urca reactions and uniform and isotropic superfluid energy gaps of different sizes.
Methods. We first do exploratory work by integrating the thermal and chemical evolution equations numerically for different energy gaps, which uncovers a rich phenomenology of stable and unstable solutions. To understand these, we perform a stability analysis around the quasisteady state, identifying the characteristic times of growing, decaying, and oscillating solutions.
Results. For small gaps, the phenomenology is similar to the previously studied cases, in the sense that the solutions quickly converge to a quasisteady state. For large gaps ( ≳ 0.05 MeV), these solutions become unstable, leading to a limitcycle behavior of periodicity ~10^{6−7} yr, in which the star is hot (T_{s} ≳ 10^{5} K) for a small fraction of the cycle (~5–20%), and cold for a longer time.
Key words: stars: neutron / dense matter / stars: rotation / pulsars: general / pulsars: individual: PSR J04374715
© ESO, 2011
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 ages ≳ 10^{5} yr.
Several proposed mechanisms could keep NSs hot beyond the standard cooling timescale ~10^{7} yr (see, e.g., Schaab et al. 1999, but note that this paper gives an incorrect parameterization of rotochemical heating, as explained in Gonzalez & Reisenegger 2010). Gonzalez & Reisenegger (2010) compared the abilities of these mechanisms to reheat NSs, concluding that two of them might be important in old NSs, namely, reheating by frictional motion of superfluid neutron vortices (Alpar et al. 1984; Shibazaki & Lamb 1989) and rotochemical heating. The latter was first proposed by Reisenegger (1995) and then improved by Fernández & Reisenegger (2005), and considered the internal structure of nonsuperfluid NSs with realistic equations of state (EOSs) in the framework of general relativity. It works as follows. As the NS’s rotation rate decreases, the reduction in the centrifugal force makes it contract. This perturbs each fluid element, raising the local pressure and causing deviations from beta equilibrium. Eventually, the system reaches a quasisteady configuration, where the rate at which spindown modifies the equilibrium concentrations is the same as that at which neutrino reactions restore the equilibrium. These reactions heat the stellar interior, making the star emit thermal radiation. The potential presence of superfluid nucleons in the NS’s interior has been widely considered to model their thermal evolution (see, e.g. Yakovlev & Pethick 2004, for cooling models) since it considerably reduces the neutrino reactions and the specific heat involving superfluid species (Yakovlev et al. 2001), and opens new neutrino emission processes, namely pair breaking and formation reactions (Flowers et al. 1976).
In a previous paper (Petrovich & Reisenegger 2010), we modeled rotochemical heating of millisecond pulsars with only modified Urca reactions in the presence of uniform and isotropic Cooper pairing gaps of neutrons Δ_{n} and protons Δ_{p}. We verified the orderofmagnitude predictions of Reisenegger (1997), finding that the chemical imbalances in the star grow up to the threshold value Δ_{thr} = min(Δ_{n} + 3Δ_{p},3Δ_{n} + Δ_{p}), which is higher than in the quasisteady state achieved in the absence of superfluidity. Therefore, the old superfluid NSs will take longer to reach the quasisteady state than their nonsuperfluid counterparts, and they have a higher a luminosity in this state, given by , where Ṗ_{20} is the period derivative in units of 10^{20} and P_{ms} is the period in milliseconds. With the previous relation, we found that energy gaps in the range 0.05 [MeV] ≲ Δ_{thr} ≲ 0.45 [MeV] are consistent with the ultraviolet emission of PSR J04374715 (Kargaltsev et al. 2004).
We extend our previous analysis to the case where the much faster direct Urca reactions are allowed. We find a qualitatively new behavior of rotochemical heating, where the temperature and chemical imbalances oscillate around the quasisteady state.
The structure of this paper is the following. In Sect. 2, we present the basic equations of rotochemical heating and the phasespace integrals for direct Urca reactions with Cooper pairing gaps and chemical imbalances. In Sect. 3, we describe our results for the thermal evolution and the linear stability analysis of old NSs. We summarize our main conclusions in Sect. 4.
2. Theoretical framework
2.1. Rotochemical heating: basic equations
The basic framework for rotochemical heating is explained in detail in Fernández & Reisenegger (2005), and the modifications made to consider the Cooper pairing effects are described in Petrovich & Reisenegger (2010). The latter is the framework we use throughout this paper. Therefore, we just point out the fundamental equations for completeness and to clarify the notation of the present paper.
We consider the simplest model of a neutron star core, composed of neutrons, protons, electrons, and muons (npeμ matter), ignoring the potential presence of exotic particles.
The internal temperature, redshifted to a distant observer, T_{∞}, is taken to be uniform inside the star because we are modeling the thermal evolution 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) (1)where C is the total heat capacity of the star, is the total power released by the heating mechanism, the total power emitted as neutrinos due to Urca reactions, and the power released as thermal photons.
The amount of energy released by each Urcatype reaction is η_{npl} = μ_{n} − μ_{p} − μ_{l} (l = e,μ), where μ_{i} is the chemical potential of the particle species i. Thus, we write the total energy dissipation rate as (2)where is the net reaction rate of the Urca reaction integrated over the core involving the lepton l.
The photon luminosity is calculated by assuming blackbody radiation , where R_{∞} and T_{s,∞} are the radius and the surface temperature of the star measured by 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, also uniform throughout the core, is given by where the terms Z_{np}, Z_{npe}, Z_{npμ}, W_{npe}, and W_{npμ} are constants that depend on the stellar structure and are kept 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).
As Petrovich & Reisenegger (2010) showed, the results of the evolution with rotochemical heating when computing and in the presence of superfluid nucleons are substantially different from their superfluid counterparts calculated by Fernández & Reisenegger (2005), since superfluidity strongly inhibits these reactions. Thus, the chemical imbalances become larger during the quasisteady state, lengthening the timescale to arrive at this state compared with the nonsuperfluid case, and predicting higher temperatures in old NSs. In this paper, we include the powerful direct Urca reactions to our previous study.
2.2. Cooper pairing
In the core, neutrons are believed to form Cooper pairs because of their interaction in the triplet states via the anisotropic channels m_{J} = 0 (type B) or m_{J} = 2 (type C), while protons form isotropic, singlet pairs (type A) (Yakovlev et al. 2001). Additionally, in the outermost core and inner crust, neutrons are believed to form singletstate pairs. The (types B and C) state description is rather uncertain in the sense that the energetically most probable state of nnpairs (m_{J} = 0,1,2) 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 due to each type of superfluidity of the net reaction rate for direct Urca and modified Urca reactions out of beta equilibrium, finding that the suppression due to type A superfluidity is of strength between the suppression due to anisotropic channels type B and type C superfluidity, respectively. For simplicity, we consider the energy gaps for the neutrons Δ_{n} and the protons Δ_{p} at zero temperature, redshifted to a distant observer, as parameters that are isotropic ( pairs) and uniform throughout the core of the NS.
The phase transition for a nucleon species into a superfluid state takes place when its temperature falls below a critical value T_{c}. This temperature is related to the energy gap at zero temperature Δ(T = 0); for the isotropic pairing channel , Δ(T = 0) = 1.764kT_{c}. Additionally, when the transition occurs, the amplitude of the energy gap depends on the temperature by means of the BCS equation (Yakovlev et al. 2001), which can be fitted by the practical formula of Levenfish & Yakovlev (1994) for the isotropic gap (5)where δ is the variable used in the phasespace integrals in Sect. 2.3. It is straightforward to check that the limiting cases are reproduced by Eq. (5), i.e. δ = 0 when T = T_{c} and δ = Δ(T = 0)/kT when T ≪ T_{c}. Levenfish & Yakovlev (1994) claim that intermediate values of T/T_{c} are also reproduced by this formula with a maximum error less than 5%, which is accurate enough for the purposes of this work.
Having defined the energy gap of the nucleon Δ_{i} with i = n, p, it is possible to express the momentum dependence of the nucleon energy ϵ_{i}(p_{i}) near the Fermi level, i.e. p_{i} − p_{Fi} ≪ p_{Fi}, as follows (Yakovlev et al. 2001) (6)where p_{i}, p_{Fi}, v_{Fi}, and μ_{i} are the momentum, the Fermi momentum, the Fermi velocity, and the chemical potential of species i = n, p, respectively.
2.3. Neutrino emissivity
The fastest reactions in NS cores are the direct Urca processes where l = e, μ. If these reactions are kinematically allowed, their emissivity is much greater than those produced by the modified Urca reactions (Yakovlev et al. 2001).
We write the neutrino emissivity and the net reaction rate of direct Urca 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 (11)where the term S_{l} 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 I_{D,ϵ} and I_{D,Γ}
are dimensionless phasespace integrals that contain the dependence of the emissivity and the net reaction rate, respectively, on the chemical imbalances and the energy gaps Δ_{n} and Δ_{p}.
To introduce these integrals, it is useful to define the usual dimensionless variables normalized by the thermal energy kT, as follows (12)which represent the energy of the nonsuperfluid degenerate particle j, the neutrino, and the chemical imbalance involving the lepton l, respectively, while for the superfluid nucleon i we write (13)where δ_{i} is defined in Eq. (5) in terms of Δ_{i}. Thus, in terms of these variables, (14)and (15)where f(·) is the Fermi function f(x) = 1/(1 + e^{x}), and the numerical factor in front of the integral is normalizes I_{D,ϵ} to 1 when the energy gaps and the chemical imbalances are zero.
In the nonsuperfluid case (i.e. δ_{n} = δ_{p} = 0), these integrals reduce to the polynomials calculated by Reisenegger (1995)The phasespace integrals I_{D,ϵ} and I_{D,Γ} cannot be solved analytically when one or two nucleon species are superfluid. Thus, Villain & Haensel (2005) and Petrovich & Reisenegger (2010) computed these integrals numerically and described the suppression produced by the Cooper pairs by means of the socalled reductions factors. The latter are defined as the ratio of the phasespace integrals with superfluid particle species to the nonsuperfluid integrals, i.e. R_{D,ϵ} = I_{D,ϵ}/F_{D}(ξ_{l}) and R_{D,Γ} = I_{D,Γ}/H_{D}(ξ_{l}). Hereafter, we calculate these integrals numerically by using the numerical method explained in Petrovich & Reisenegger (2010), which is based on the GaussLaguerre quadrature, taking advantage of the exponentially decaying behavior of the Fermi functions appearing in the integrand.
3. Results and discussion
3.1. Evolution
In this section, we show the peculiar behavior of rotochemical heating in MSPs with superfluid nucleons and direct Urca reactions. To do so, we numerically compute the evolution of Eqs. (1), (3), and (4). Hereafter, we omit the ∞ subscript in the temperature and the ∞ superscript in the chemical imbalances to simplify the notation.
The NS interior structure is modeled by the BPAL21 EOS of Prakash et al. (1988) for the core, supplemented with those of Pethick et al. (1995) and Haensel & Pichon (1994) for the inner and outer crust, respectively. This model opens electron direct Urca reactions in the core when the NS mass is greater than 1.67 M_{⊙} and muon direct Urca reactions when its mass is greater than 1.69 M_{⊙}, and reaches a maximum mass value of 1.71 M_{⊙}. Thus, the direct Urca reactions are allowed within the mass range 1.76 ± 0.2 M_{⊙} of the MSP J04374715 (Verbiest et al. 2008), the only MSP whose thermal radiation has been measured so far and, therefore, there is an estimate of its surface temperature (Kargaltsev et al. 2004).
For numerical calculations, the evolution of is computed by assuming magnetic dipole braking with no field decay, and relating the magnetic field on the magnetic equator, rotation period, and period derivative by the conventional formula B = 3.2 × 10^{19}(PṖ)^{1/2} G, where P is measured in seconds. The magnetic field is chosen to be B = 2.8 × 10^{8} G, so as to match the currently measured values of P and Ṗ of PSR J04374715.
In Fig. 1, we show the evolution for three different values of the neutron energy gap, with nonsuperfluid protons and no muons, illustrating the appearance of strong thermal oscillations at larger gaps. The left top panel, where Δ_{n} = 0.01 MeV, shows that rotochemical heating converges rapidly to a stable solution given by the quasisteady solution calculated from , as it happens when no superfluid nucleon species are considered (Fernández & Reisenegger 2005) or when superfluid nucleons are present, but only modified Urca reactions are allowed (Petrovich & Reisenegger 2010). A hint of an oscillation is seen before the solution converges.
In the left bottom panel of Fig. 1, we increase the value of the energy gap to Δ_{n} = 0.03 MeV. In this case, there is a clearly visible damped oscillation before the variables settle to their quasistate values.
Fig. 1 Evolution of the internal temperature T and the chemical imbalance η_{npe} for a 1.68 M_{⊙} star without muons built with the BPAL21 EOS (Prakash et al. 1988), with the initial conditions T = 10^{8} K and η_{npe} = 0. The spindown is assumed to be caused by magnetic dipole radiation, with dipolar field strength B = 2.8 × 10^{8} G and initial period P_{0} = 1 ms. The dashed lines are the quasisteady state solution calculated from . In all cases, the protons are taken to be nonsuperfluid, while the neutron superfluid energy gaps are Δ_{n} = 0.01 MeV (upper left panel), Δ_{n} = 0.03 MeV (lower left panel), and Δ_{n} = 0.05 MeV (upper and lower right panels). The lower right panel shows a closeup of the second peak in temperature of the evolution plot shown in the upper right panel (Δ_{n} = 0.05 MeV) with a linear scale in time. 

Open with DEXTER 
In the upper right panel of Fig. 1, we use Δ_{n} = 0.05 MeV to show strong oscillations without a clear damping effect. This behavior can be understood as follows. Once the chemical imbalance increases to a critical value ~Δ_{n}, the reaction rates increase dramatically, suddenly increasing the temperature and thus decreasing the imbalance. At this point, the imbalance is substantially below the energy gap and the direct Urca reactions are strongly inhibited (see Petrovich & Reisenegger 2010, for more details on the reduction factors). Thus, the star begins to cool almost exclusively by photon emission, and the spindown compression forces the chemical imbalance to grow almost unimpeded by the neutrino reactions until it reaches its critical value again. This process is repeated cyclically in a timescale given, essentially, by the time taken by the chemical imbalance to come back to its critical value ~Δ_{n}, where it erases any record of the previous cycles.
The lower right panel of Fig. 1, shows a closeup of one of the peaks depicting the evolution with Δ_{n} = 0.05 MeV in order to resolve and illustrate the rather sharp peaks. We also illustrate in this figure that the timescale involved in the rapid increase in temperature is on the order of ~10^{5} yr, while the decrease happens on a timescale of ~10^{6} yr.
This unexpected behavior will be the focus of our analysis.
3.2. Stability analysis
In this section, we analyze the local stability of the temperature and chemical imbalance by perturbing the evolution Eqs. (1) and (3) around the quasisteady state and computing the growth rates of these small displacements. For simplicity, we ignore the presence of muons in our analysis.
We can safely neglect the modified Urca reactions of the electrons because the direct Urca reactions are much more efficient and the superfluid suppression is similar in both cases. Thus, the evolution equations (with η_{npe} in units of temperature) can be written as We considered the fully accreted envelope model of Potekhin et al. (1997) to express the photon luminosity in terms of the internal temperature as , where is a numerical factor that depends on the stellar structure, and α = 2.42. Additionally, we defined the functions ℱ and to save notation in the subsequent analysis.
We write the perturbed solutions to these equations as (20)where T^{qs} and are the solutions of , which are unique for each combination of input parameters, i.e. energy gaps, spindown power, and structure constants. We define γ as the growth rate of the perturbations δT and δη_{npe}. Replacing these perturbed solutions in the differential equations leads to the usual eigenvalue problem for the growth rates of these perturbations (21)where I is the 2 × 2 identity matrix and J is the Jacobian matrix resulting from the linearization of Eqs. (18) and (19) and evaluated at T^{qs} and . Thus, we write the Jacobian matrix elements, replacing the quasisteady state and taking logarithmic derivatives, as We solve for the two growth rates γ_{ + } and γ_{ − } by simply calculating (26)where τ and Δ are the trace and the determinant of J, respectively. The stability of this system of differential equations, at least when the linear approximations is valid, i.e. close to quasisteady solutions, depends on these growth rates. Thus, if Re(γ_{ ± }) < 0, the fixed point solutions will be stable, as happens when we have only modified Urca reactions or no superfluid nucleons are considered. Otherwise, if Re(γ_{ ± }) > 0 the solutions will diverge. Moreover, the solution will oscillate if Im(γ_{ ± }) ≠ 0, and the oscillations become important during the evolution if Im(γ_{ ± }) are comparable to or greater than Re(γ_{ ± }).
Fig. 2 Real and imaginary parts of the growth rates γ_{ + } and γ_{ − } (upper panel) calculated from Eq. (26) and components of the Jacobian matrix (lower panel) in Eqs. (22)–(25), for different values of Δ_{n}. In both panels, the stellar parameters are the same as in Fig. 1 and we assume that . 

Open with DEXTER 
In the upper panel of Fig. 2, we show the real and imaginary parts of the two growth rates for different values of Δ_{n}. We fix the value of the spindown parameters to in order to compare our results to those of Fig. 1, where reaches this value at the age of 5 × 10^{6} yr and slowly changes until ~ 10^{8} yr. When the energy gap is Δ_{n} = 0.01 MeV (top panel of Fig. 1), the system is locally stable because Re(γ_{ ± }) < 0, and the oscillation is highly damped because  Re(γ_{ ± })  ≫  Im(γ_{ ± })  . Then, if we take Δ_{n} = 0.03 MeV (middle panel of Fig. 1), we can see that the system is still locally stable, but Re(γ_{ + }) and Re(γ_{−}) increase in such a way that  Re(γ_{ ± })  <  Im(γ_{ ± })  /(2π), and therefore the oscillations are less damped than those with lower energy gaps. Up to here, the linearized system predicts that the fixed point T^{qs} and behaves as a stable spiral(also known as an attractor or a sink) in the T − η_{npe} space, hence the linear analysis ensures that it is also correct for the nonlinear system close to this fixed point (Strogatz 1994).
Increasing the energy gap to a value higher than Δ_{n} = 0.03 MeV, the linearized system becomes unstable because Re(γ_{ ± }) changes from being negative to positive. Beyond this point, the nonlinear system acts as an unstable spiral (also known as a repeller or source) in the T − η_{npe} space. If the initial condition of the system of Eqs. (18) and (19) is close to the fixed point, it will diverge from it. Afterwards, the amplitude of the perturbations becomes large enough to make the nonlinear terms significant, and the linear analysis is no longer valid. We show in section Sect. 3.4 that, for high energy gaps, the highly nonlinear system converges to a stable limit cycle.
3.3. Superfluidity driving the oscillations
The oscillations are exclusively caused by the appearance of the superfluid energy gaps at the Fermi level of the nucleons. This drastically changes the dependence of the neutrino emissivity and the net reaction rate on the internal temperature and chemical imbalance, which causes the oscillations, as we show in this section.
In the lower panel of Fig. 2, we show the previous effect by means of the components of the Jacobian matrix in Eqs. (22)–(25), evaluated at the quasisteady state. From here, we can observe that close to the quasisteady state, the evolution represented by Eqs. (18) and (19) becomes unstable, essentially, because the component J_{11} = ∂ℱ/∂T changes from being negative to positive as we increase Δ_{n}, changing the sign of the trace of the Jacobian τ.
To understand why this happens, we analyze the dimensionless net heating function ξ_{e}I_{Γ} − I_{ϵ} that appears in Eq. (18). For the nonsuperfluid case, one finds from the polynomials given in Eqs. (16) and (17) that, in the quasisteady state, where T ≪ η, ξ_{e}I_{Γ} − I_{ϵ} ~ 21/(457π^{6})(η/T)^{6}, which gives ∂ln(ξ_{e}I_{Γ} − I_{ϵ})/∂lnT ~ − 6 , i.e. from Eq. (22) J remains negative. In contrast, for the superfluid case, we find that this function can increase with temperature for values of the energy gaps sufficiently greater than the temperature. We know that both of the functions I_{Γ} and I_{ϵ} increase with temperature since this increases the phasespace available to the particles. However, as Petrovich & Reisenegger (2010) showed by means of the reduction factors, superfluidity blocks the phasespace of the integrals I_{ϵ} and I_{Γ} quite differently, and for the quasisteady state we find that R_{ϵ} ≪ R_{Γ} (see Fig. 2 of Petrovich & Reisenegger 2010), which forces ξ_{e}I_{Γ} − I_{ϵ} to be an increasing function of the temperature. Finally, this happens until the temperature is high enough to make the reduction factors comparable and the function decreases with temperature in a similar way to its nonsuperfluid counterpart.
3.4. Limit cycle in η_{npe} – T: existence and Δdependence
In the upper panel of Fig. 3, we show the evolution of rotochemical heating in the η_{npe} − T space for different initial conditions that converge to the same the limit cycle. Moreover, we know that the trajectories inside the limit cycle are repelled by the fixed point of ℱ and , which in this case is an unstable spiral, as discussed in Sect. 3.2. Thus, we could construct a trapping region, excluding the fixed point, such as the vector that field points inward at the boundary of this region and, therefore, argue by means of the PoincaréBendixson theorem that there is a closed orbit, i.e. a stable limit cycle (see, e.g. Strogatz 1994, for details on limit cycles).
Fig. 3 Evolution of rotochemical heating in the η_{npe} − T space for a superfluid star with the same parameters of Fig. 1. Upper panel: evolution with Δ_{n} = 0.05 MeV for different initial conditions (dotteddashed lines) that converge to the same limit cycle (solid line). The quasisteady state is indicated (star). Lower panel: evolution of the limit cycle for Δ_{n} = 0.04, 0.05, 0.06 MeV. The arrows indicate the path of the cycles AB, BC, and CA for the star with Δ_{n} = 0.05 MeV. 

Open with DEXTER 
In the lower panel of Fig. 3, we show different limit cycles by changing the energy gap and label three sections of the evolution curve with Δ_{n} = 0.05 MeV, namely AB, BC, and CA. Along the path AB, the chemical imbalance is ~Δ_{n} and some reactions are opened, increasing the temperature until it reaches node B, without changing the chemical imbalance. Along the path BC, the temperature is high enough to open abruptly many reactions that tend to restore the beta equilibrium, and η_{npe} drops rapidly. Finally, along the path CA, the chemical imbalance is sufficiently below the energy gap to strongly inhibit the neutrino reactions. Thus, no heating mechanism is present and the star cools by photon emission. Moreover, no restoring mechanism for the beta equilibrium is present and the chemical imbalance grows by the spindown compression.
We observe from this figure that the amplitude of the limit cycles increases with the gap Δ_{n}, which we can easily understand from the previous analysis. For the path AB, the chemical imbalance must be larger for larger Δ_{n} to open reactions. Then, from B to C, we note from the linear analysis in Sect. 3.2 that for larger gaps the system becomes more unstable or, more precisely, the heating function ξ_{e}I_{Γ} − I_{ϵ} becomes more sensitive to small perturbations in temperature reaching higher temperatures, forcing the chemical imbalances to reach lower values in node C. Finally, from C to A, all the curves are parallel since the photon emission and spindown compression do not depend on Δ_{n}.
3.5. Limit cycle timescales and the MSP J0437
We can observe from the bottom panel of Fig. 1 that there are two very different scales governing the limit cycles. The first and shorter one involves the rapid increase in temperature and ends with a rapid decrease in the chemical imbalance, i.e. it goes from A to C in the lower part of Fig. 3. The second scale goes from C to A and is substantially dominant in the rotochemical heating evolution.
The first is quite stable and on the order of the timescales shown in Fig. 2 for the growth rates γ, i.e. a timescale of ~10^{5} yr. However, the second scale depends strongly on the energy gap, because the latter sets the amplitude of the drop in η_{npe}. It also depends on , since it accounts for the rate at which the chemical imbalance grows. For the examples displayed in Figs. 1 and 3, the longer timescale is ~10^{6}−10^{7} yr.
Considering the large difference between both timescales, we check that the fraction of the time in which the thermal emission would be detectable, i.e. the surface temperature T_{s} ≳ 10^{5} K, is in the range ~5–25%, mainly depending on the energy gaps and spindown parameters. Moreover, it is also unlikely to find it close to its quasisteady state and therefore are unable to draw further conclusions about the pairing gaps needed to explain observations as in Petrovich & Reisenegger (2010).
However, we can account for the maximum energy gap combination Δ_{thr} = Δ_{n} + Δ_{p} at which rotochemical heating evolves to a stable solution, i.e. its quasisteady state, and predict a surface temperature. Thus, for the spindown parameters of MSP J0437 and the set of EOSs from Prakash et al. (1988) that allow direct Urca reactions of electrons and muons for its mass range 1.76 ± 0.2 M_{⊙} (Deller et al. 2008), namely BPAL21, BPAL31, BPAL 32, and BPAL33, the maximum gap combination that predict a quasisteady state is in the range of Δ_{thr} ~ 0.01−0.1 MeV. Thus, the temperatures predicted when including superfluidity are a few times higher than those predicted by the nonsuperfluid case (Fernández & Reisenegger 2005). They are quite close to 10^{5} K and for some cases can explain the likely thermal emission of the MSP J0437 (Kargaltsev et al. 2004).
Finally, given the rapid increase in temperature, we could ask whether the finite thermal diffusion timescales are relevant to our study or make any difference to the validity of our assumption of an isothermal NS interior. We claim that this assumption is fairly safe, since the increase in temperature takes ~10^{5} yr, while the thermal diffusion time is just a few decades for these low temperatures.
Fig. 4 Evolution of the internal temperature T and the chemical imbalances η_{npe} and η_{npμ} with the initial condition T = 10^{8} K and null chemical imbalances. The spindown is assumed to be caused by the magnetic dipole radiation, with dipolar field strength B = 2.8 × 10^{8} G and initial period of P_{0} = 1 ms. The dashed lines indicate the energy gap Δ_{n} = 0.05 MeV. Upper panel: evolution of a 1.69 M_{⊙} star built with the BPAL21 EOS (Prakash et al. 1988), which allows muon direct Urca reactions, and we add the vortex creep reheating term of Eq. (27) with an excess of angular momentum J = 10^{43} erg s. Lower panel: evolution of a 1.68 M_{⊙} star built with the BPAL21 EOS (Prakash et al. 1988), which forbids direct Urca reactions with muons. 

Open with DEXTER 
3.6. Effect of muons and other reheating mechanisms.
The presence of muons does not introduce changes to the previous analysis, as long as direct Urca reactions with muons are allowed. On the other hand, if they react only by means of the modified Urca channels, the oscillations are damped and finally vanish.
In the upper panel of Fig. 4, we observe that the only effect to be considered is that the oscillations would be driven by the reactions involving muons, instead of those with electrons. This happens because the departure from beta equilibrium with muons, i.e. η_{npμ}, is more strongly affected by the spindown compression, and the direct Urca reactions with muons will lead the oscillations since they reach the critical value ~Δ_{n} first.
Moreover, including another reheating mechanism does not alter the existence of the oscillations as long as the temperature increase produced by this mechanism is not so high as to overcome the quasisteady value of the temperature. As Gonzalez & Reisenegger (2010) argued, the two most relevant mechanisms to reheat an old neutron star are rotochemical heating and vortex creep (Alpar et al. 1984; Shibazaki & Lamb 1989). Therefore we include the latter in the upper panel of Fig. 4 to show that the oscillations persist even when another reheating mechanism is included. The energy dissipation rate that must be included as a heating term in Eq. (1) is given by (Alpar et al. 1984) (27)where J is the excess of angular momentum. By choosing an intermediate value of J = 10^{43} erg s for the models in Shibazaki & Lamb (1989), we observe that the only effect of the extra heating term is to provide a higher floor to the internal temperature. This results in a straightforward way from the analysis in Sect. 3.2, since the Jacobian element in Eq. (22) is unaltered by the new term in Eq. (27). However, if this mechanism keeps the NS at substantially higher temperatures, it could displace the quasisteady state and change the behavior of the oscillations.
In the lower panel of Fig. 4, we show the evolution when direct Urca reactions with muons are forbidden, but modified Urca reactions involving muons are present. At first, the chemical imbalance of the muons reaches a value close to the energy gaps and reheats the star by means of modified Urca reactions. The chemical imbalance of the electrons then reaches this threshold value and the oscillations start to work normally until the chemical imbalance of the muons reaches its quasisteady state, which predicts higher temperatures than those without muons (see the bottom panel of Fig. 1). Hence, the oscillations vanish and the quasisteady state is mostly governed by the chemical imbalance of the muons.
Fig. 5 Evolution of the internal temperature T and the chemical imbalance η_{npe} for a 1.68 M_{⊙} star built with the EOS BPAL21 and no muons, with the initial condition T = 10^{9} K and null chemical imbalance. The spindown is assumed to be caused by the magnetic dipole radiation, with dipolar field strengths B = 10^{11} G (dashed line) and B = 10^{12} G (solid line) initial period of P_{0} = 5 ms. The dotteddashed line indicates the energy gap Δ_{n} = 0.05 MeV. 

Open with DEXTER 
3.7. Effect on classical pulsars
In this section, we show the effects of these oscillations on classical pulsars. In Fig. 5, we illustrate their unstable behavior by showing the evolution of rotochemical heating for a NS with an energy gap Δ_{n} = 0.05 MeV, electron direct Urca reactions, no muons, and two magnetic fields of B = 10^{11} G and B = 10^{12} G.
From this figure, we observe that the temperature and the chemical imbalance can oscillate only during half of the period for B = 10^{11} G or a full period for B = 10^{12} G since the rotation rate is insufficient for the imbalance to recover from its drop and, therefore, produce any other oscillation. Thus, after the sudden increase in temperature at ~10^{5} yr, the star cools almost exclusively by means of photon emission, while the chemical imbalance grows slowly, never reaching the threshold of the energy gap again.
If we changed the initial period to lower values (e.g. P_{0} = 1 ms), we could observe from the case with B = 10^{11} G one full oscillation period and another half, similar to the lower panel of Fig. 4. Similarly, lower values of the energy gap could produce more oscillations. In contrast, if we increase the initial period, the rotation energy to increase the chemical imbalance may not be enough to reach the energy gap threshold, and oscillation may not occur at all.
In conclusion, the oscillating behavior generally does not persist in classical pulsars, and the standard cooling governs their thermal evolution.
4. Conclusions
We have studied a new kind of oscillations of the temperature and chemical imbalances, produced during the evolution of rotochemical heating when the direct Urca reactions are allowed in superfluid MSPs with relatively high uniform and isotropic Cooper pairing gaps.
The oscillations work as follows. The direct Urca reactions are strongly blocked until the chemical imbalances grow up to a threshold value Δ_{thr} = Δ_{n} + Δ_{p}. Soon after this happens, strong reactions are suddenly turned on, increasing the temperature and decreasing the chemical imbalances. As this happens, the reactions are again strongly blocked and the star cools down by photon emission, while the chemical imbalances again increase until they reach the threshold value required to repeat the cycle after 10^{6 − 7} yr. The temperature stays high for only a
small fraction of the cycle, making it difficult to detect this effect and predict a NS temperatures at a given time. For gaps below ~ 0.05 MeV or with muons reacting only by modified Urca processes, the oscillations vanish and the system reaches a quasisteady state, as previously found in the nonsuperfluid or modified Urca cases.
Acknowledgments
We thank Denis González and Nicolás González 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), Proyecto Basal PFB06/2007, and Proyecto Límite VRAID No 15/2010.
References
 Alpar, M. A., Anderson, P. W., Pines, D., & Shaham, J. 1984, A&A, 276, 325 [NASA ADS] [CrossRef] (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)
 Deller, A., Verbiest, J., Tingay, S., & Bailes, M., 2008, ApJ, 685, L67 [NASA ADS] [CrossRef] (In the text)
 Fernández, R., & Reisenegger, A. 2005, ApJ, 625, 291 [NASA ADS] [CrossRef] (In the text)
 Flowers, E., Ruderman, M., & Sutherland, P. 1976, A&A, 205, 541 [NASA ADS] [CrossRef] (In the text)
 Gonzalez, D., & Reisenegger, A., 2010, A&A, 522, A16 [NASA ADS] [CrossRef] [EDP Sciences] (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, ApJ, 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. Notes Phys., 578, 30 [NASA ADS] [CrossRef] (In the text)
 Maxwell, O.V., 1979, A&A, 231, 201 [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)
 Petrovich, C., & Reisenegger, A., 2010, A&A, 521, A77 [NASA ADS] [CrossRef] [EDP Sciences] (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, Phys. Rev. Lett., 61, 2518 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Reisenegger, A. 1995, ApJ, 442, 749 [NASA ADS] [CrossRef] (In the text)
 Reisenegger, A. 1997, ApJ, 485, 313 [NASA ADS] [CrossRef] (In the text)
 Reisenegger, A., Jofré, P., Fernández, R., & Kantor, E. 2006, ApJ, 653, 568 [NASA ADS] [CrossRef] (In the text)
 Schaab, Ch., Sedrakian, A., Weber, F., & Weigel, M. K. 1999, A&A, 346, 465 [NASA ADS] (In the text)
 Shibazaki, N., & Lamb, F. K. 1989, ApJ, 346, 808 [NASA ADS] [CrossRef] (In the text)
 Strogatz, S. H. 1994, Nonlinear dynamics and chaos (AddisonWesley) (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, ApJ, 679, 675 [NASA ADS] [CrossRef] (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)
All Figures
Fig. 1 Evolution of the internal temperature T and the chemical imbalance η_{npe} for a 1.68 M_{⊙} star without muons built with the BPAL21 EOS (Prakash et al. 1988), with the initial conditions T = 10^{8} K and η_{npe} = 0. The spindown is assumed to be caused by magnetic dipole radiation, with dipolar field strength B = 2.8 × 10^{8} G and initial period P_{0} = 1 ms. The dashed lines are the quasisteady state solution calculated from . In all cases, the protons are taken to be nonsuperfluid, while the neutron superfluid energy gaps are Δ_{n} = 0.01 MeV (upper left panel), Δ_{n} = 0.03 MeV (lower left panel), and Δ_{n} = 0.05 MeV (upper and lower right panels). The lower right panel shows a closeup of the second peak in temperature of the evolution plot shown in the upper right panel (Δ_{n} = 0.05 MeV) with a linear scale in time. 

Open with DEXTER  
In the text 
Fig. 2 Real and imaginary parts of the growth rates γ_{ + } and γ_{ − } (upper panel) calculated from Eq. (26) and components of the Jacobian matrix (lower panel) in Eqs. (22)–(25), for different values of Δ_{n}. In both panels, the stellar parameters are the same as in Fig. 1 and we assume that . 

Open with DEXTER  
In the text 
Fig. 3 Evolution of rotochemical heating in the η_{npe} − T space for a superfluid star with the same parameters of Fig. 1. Upper panel: evolution with Δ_{n} = 0.05 MeV for different initial conditions (dotteddashed lines) that converge to the same limit cycle (solid line). The quasisteady state is indicated (star). Lower panel: evolution of the limit cycle for Δ_{n} = 0.04, 0.05, 0.06 MeV. The arrows indicate the path of the cycles AB, BC, and CA for the star with Δ_{n} = 0.05 MeV. 

Open with DEXTER  
In the text 
Fig. 4 Evolution of the internal temperature T and the chemical imbalances η_{npe} and η_{npμ} with the initial condition T = 10^{8} K and null chemical imbalances. The spindown is assumed to be caused by the magnetic dipole radiation, with dipolar field strength B = 2.8 × 10^{8} G and initial period of P_{0} = 1 ms. The dashed lines indicate the energy gap Δ_{n} = 0.05 MeV. Upper panel: evolution of a 1.69 M_{⊙} star built with the BPAL21 EOS (Prakash et al. 1988), which allows muon direct Urca reactions, and we add the vortex creep reheating term of Eq. (27) with an excess of angular momentum J = 10^{43} erg s. Lower panel: evolution of a 1.68 M_{⊙} star built with the BPAL21 EOS (Prakash et al. 1988), which forbids direct Urca reactions with muons. 

Open with DEXTER  
In the text 
Fig. 5 Evolution of the internal temperature T and the chemical imbalance η_{npe} for a 1.68 M_{⊙} star built with the EOS BPAL21 and no muons, with the initial condition T = 10^{9} K and null chemical imbalance. The spindown is assumed to be caused by the magnetic dipole radiation, with dipolar field strengths B = 10^{11} G (dashed line) and B = 10^{12} G (solid line) initial period of P_{0} = 5 ms. The dotteddashed line indicates the energy gap Δ_{n} = 0.05 MeV. 

Open with DEXTER  
In the text 