What the Timing of Millisecond Pulsars Can Teach us about Their Interior

What the Timing of Millisecond Pulsars Can Teach us about Their Interior

Mark G. Alford and Kai Schwenzer Department of Physics, Washington University, St. Louis, Missouri, 63130, USA

The cores of compact stars reach the highest densities in nature and therefore could consist of novel phases of matter. We demonstrate via a detailed analysis of pulsar evolution that precise pulsar timing data can constrain the star’s composition, through unstable global oscillations (r-modes) whose damping is determined by microscopic properties of the interior. If not efficiently damped, these modes emit gravitational waves that quickly spin down a millisecond pulsar. As a first application of this general method, we find that ungapped interacting quark matter is consistent with both the observed radio and x-ray data, whereas for ordinary nuclear matter some additional enhanced damping mechanism is required.

Pulsars are believed to be ultradense compact objects that may consist of nuclear matter Lattimer:2004pg () or may include more exotic material such as quark matter Itoh:1970uw (); Witten:1984rs (); Alford:2007xm (). There is a wealth of very precise radio—and increasingly also high-energy TheFermi-LAT:2013ssa ()—pulsar timing data Manchester:2004bp (), showing that they are extremely stable systems with known frequency and spindown (SD) rate. The goal of this paper is to show that these data can be used to constrain hypotheses about the interior composition of the star. Our approach relies on r-modes (RMs) Andersson:1997xt (); Andersson:2000mf (), global oscillations which are unstable via the Friedman-Schutz mechanism Friedman:1978hf (). If they are not effectively damped, r-modes grow spontaneously and emit gravitational waves (GWs), spinning the star down. Different possible phases of dense matter have different viscosities, and hence differ in their ability to damp r-modes. Therefore observations of high-spin pulsars indicate that sufficiently strong damping must be present, constraining the possible phases of matter in the star.

The macroscopic state of the star is specified by its angular velocity , its core temperature , the amplitude of the r-mode (which is unobservable). The evolution is determined by conservation equations which involve energy loss rates, namely the power fed into the r-mode by radiating gravitational waves, the dissipated power that heats the star and the thermal luminosity that cools it. If not globally so at least over certain temperature ranges, they follow power laws which for the r-mode read


where at small , are logarithmic correction factors that arise from non-Fermi liquid (NFL) effects in certain forms of quark matter Schwenzer:2012ga () and we neglect an amplitude dependence of the luminosity since it is only relevant for the unrealistic case Alford:2012yn (). The prefactors in (1) are given in tab. 2. They are given by a few dimensionless parameters that encode the relevant properties of the star, see tab. 1. We will first review and refine constraints on the composition from measurements of and , concluding that currently known damping mechanisms have difficulty explaining the pulsar data within minimal hadronic matter models of neutron stars, which include only viscous damping which can be reliably estimated. We then show how this conclusion is confirmed and enhanced by measurements of and (“timing data”).

compact star
SS (NFL) - -
Table 1: Parameters characterizing the strange star (SS) and neutron star (NS) considered in this work Alford:2010fd (). The exponents , , (,,) arise in the parameterizations eq. (1) for mechanisms in tab. 2 along with the corresponding constants , , .
GW luminosity
Shear viscosity (SV)
Bulk viscosity (BV, low )
Ekman layer (EL)
Neutrino luminosity
Photon luminosity
Table 2: Parameters in the general parameterization eq. (1) for the energy loss rates, in terms of the star’s mass and radius , the gravitational constant , generic scales and . For Ekman damping parameters see Lindblom:2000gu (). The dimensionless constants , and encode properties of the stellar interior Alford:2012yn () and are evaluated in table 1.

The left panel of fig. 1 shows - data for low mass x-ray binaries (LMXBs) Haskell:2012 (), which are being heated and potentially spun up by accretion from a companion. is the core temperature, inferred from X-ray spectra using a model of the envelope Gudmundsson:1983ApJ (). These involve uncertainties (estimated by the error-bars) or provide only upper limits (left-pointing arrows). The figure also shows static instability boundaries Andersson:2000mf () for a few hypothesized star compositions. The boundaries are determined by and explicitly given for the individual segments with a given dominant damping mechanism by Alford:2010fd (); Reisenegger:2003cq ()


The region above a boundary is where dissipation is insufficient to damp the r-modes. The dissipation arises from shear viscosity Shternin:2008es (); Heiselberg:1993cr (), bulk viscosity Sawyer:1989dp (); Alford:2010gw (); Schwenzer:2012ga () or another mechanism like surface rubbing in a viscous boundary layer at a solid crust Lindblom:2000gu ().

The solid line is the instability boundary for a model of interacting ungapped quark matter Schwenzer:2012ga () which is compatible with the data via the no r-mode scenario, where r-modes are completely damped; this is due to the resonant enhancement of bulk viscosity Madsen:1998qb (); Madsen:1999ci (); Schwenzer:2012ga () which creates a large stability window at . In contrast, a model of non-interacting quark matter (short-dashed) does not explain the data.

The long-dashed line is the instability boundary for stars made of hadronic matter, taking into account viscous damping only. Most of the data points lie above this line, indicating that this model would leave r-modes unsuppressed. Even if we add maximum viscous damping at the crust-core boundary Lindblom:2000gu (), requiring an implausibly thin (cm-size) Ekman layer, using the improved shear viscosity result Shternin:2008es () we still get an instability line (dotted) that is below some points 111Magnetic fields can enhance Ekman damping, but have a minor impact for in old ms-pulsars Mendell:2001nz ().. Therefore, the hadronic matter model is only compatible with the data if there is some additional damping mechanism, or in a saturated r-mode scenario, where non-linear damping () limits r-modes to a tiny amplitude , determined by the condition . We use a general power-law parametrization of the saturation amplitude as realized for proposed mechanisms Bondarescu:2013xwa (); Haskell:2013hja (); Lindblom:2000az ().

Figure 1: Boundaries of the r-mode instability regions for different star compositions compared to pulsar data. Left: Standard static instability boundary compared to x-ray data Haskell:2012 (); Tomsick:2004pf () with error estimates from different envelope models Gudmundsson:1983ApJ (); Potekhin:1997 (). Right: Dynamic instability boundary in timing parameter space compared to radio data Manchester:2004bp () (all data points are upper limits for the r-mode component of the spindown). The curves represent: neutron star with standard viscous damping Shternin:2008es (); Sawyer:1989dp () (long-dashed) and with additional boundary layer rubbing Lindblom:2000gu () at a rigid crust (dotted) as well as strange star Madsen:1992sx (); Heiselberg:1993cr () (short-dashed) and same with long-ranged NFL interactions causing enhanced damping Schwenzer:2012ga () (using the strong coupling ) (solid)—more massive stars are not qualitatively different. The thin curves show for the neutron star exemplarily the analytic approximation for the individual segments. The encircled points denote the only ms-radio pulsar J0437−4715 with a temperature estimate and the only LMXB IGR J00291+5934 that has been observed to spin down during quiesence.

To see how small in hadronic matter has to be, we need to calculate the spindown evolution. It is crucial to note that the thermal evolution is always faster than the spindown Alford:2012yn (), so the temperature reaches a steady state where cooling matches heating Bondarescu:2013xwa (), , giving

Figure 2: The thermal steady-state spindown curves, along which the star evolves, for several - and -independent saturation amplitudes. Shown are also the boundaries of the instability regions for neutron stars with different damping sources as discussed in fig. 1. Left: Static instability boundary compared to x-ray data Haskell:2012 (); Tomsick:2004pf (). Evolution curves are shown for a neutron star with modified Urca cooling. The vertical line gives the temperature below which photon emission replaces neutrino emission as the dominant cooling mechanism. Right: Dynamic instability boundary in timing parameter space compared to radio data Manchester:2004bp ().

In Fig. 2 (left panel) we plot the same LMXB observations along with the spindown curves eq. (3) for hadronic matter for a range of values of . We assume photon and modified Urca cooling, with independent of and Owen:1998xg (). For surface luminosity we take the surface temperature to be related to core temperature via the unaccreted envelope model Gudmundsson:1983ApJ (), , with , . The data points are for LMXBs which are heated by accretion, so they can lie to the right of the spindown curves. We conclude that for these sources and similar bounds were obtained in Mahmoodifar:2013quw (). Moreover, it is expected that the saturation mechanism is insensitive to the detailed star configuration (mass, radius, magnetic field, …) in which case the lower bound should approximately hold for all sources. No saturation mechanism proposed so far gives such a low Lindblom:2000az (); Bondarescu:2013xwa (); Haskell:2013hja (), so a new mechanism would be required to make the data compatible, via the saturated r-mode scenario, with the interior of the star being hadronic matter. Modifying our assumptions about the saturation and cooling mechanisms does not qualitatively change this conclusion. Presently proposed saturation mechanisms allow to depend on and to negative powers Alford:2012yn (), which makes the curves steeper but the intersection with the boundary of the instability region is invariant Alford:2012yn (), so the constraints on are only slightly weakened. Direct Urca cooling Yakovlev:2004iq () gives a slightly more restrictive limit. The crust model, e.g. with accreted envelope Potekhin:1997 (), has a minor impact on the results.

We now turn to the timing data, a much larger data set of and for millisecond (ms) radio pulsars whose temperatures are generally unknown. The r-mode spindown rate is Owen:1998xg (); Alford:2012yn (), where is the moment of inertia of the star. Along the thermal steady state eq. (3) this yields the effective spindown equation Alford:2012yn () in terms of the effective braking index . Inverting it we find the evolution path in a --plot


This equation is valid even if other spindown mechanisms—like magnetic braking—are present, since for them the lost rotational energy does not heat the star. By analyzing where the evolution leaves the static instability region we obtain novel dynamic instability boundaries in --space. The result is in the case


Note that this expression depends on the cooling behavior, but is, like eq. (2), completely independent of the saturation mechanism and amplitude. These analytic expressions exhibit the complete dependence on the underlying physics and allow us to make quantitative predictions with control over the uncertainties.

The timing data are plotted in the right panel of fig. 1. The horizontal axis shows the amount of spindown due to r-modes, so every point is an upper limit, since other mechanisms might contribute to the spindown. The dynamic instability boundaries eq. (5) are also plotted for the previously considered star compositions 222The envelope model Gudmundsson:1983ApJ () is used both for neutron stars and for strange quark stars (which are assumed to have a hadronic crust suspended by electrostatic forces Alcock:1986hz ()).. We see that the ungapped quark matter model of ref. Schwenzer:2012ga () (solid (blue) line) is compatible with the timing data, since r-modes cannot be the dominant spindown mechanism. As in the - plot, this is thanks to a large stability window. However, the radio data cannot be explained by the no r-mode scenario, because once the accretion stops theses sources quickly cool until they reach the boundary of the instability region. If the saturated thermal steady state curve eq. (3) is to the right of the boundary eq. (2), the stars will periodically be heated out of the instability region and cool back in again Andersson:2001ev (). In this boundary-straddling scenario the r-mode amplitude is eventually dynamically set by thermal balance at the boundary Reisenegger:2003cq () without the need for a non-linear saturation mechanism. We therefore expect that non-accreting sources such as ms-pulsars will cluster along the low-temperature boundary of the stability window in the - plane. Since there are generally no temperature estimates for such sources we cannot yet test this prediction. However, the timing data impose a constraint: these sources should be found on or to the right (higher ) of the stability boundary in the - plane because other spindown mechanisms may also be operating. This is exactly what the data (fig. 1) show. The interacting quark matter model is therefore compatible with the timing data.

Concerning the hadronic matter model (long-dashed line) we first note that the clustering of young radio sources Manchester:2004bp () just below the hadronic matter instability boundary at means that the hadronic model can—if in young pulsars is sufficiently large—explain why they do not spin faster Alford:2012yn (). However, the large column of data points around , consisting of old ms-pulsars, is a problem for the minimal hadronic matter model, since these upper limits are within its dynamic instability region. In principle this might be explained in a saturated r-mode scenario, but as we see from the thermal steady state curves in the right panel of fig. 2, this scenario requires , a similar value to that obtained from the - data, which is not provided by current mechanisms Lindblom:2000az (); Bondarescu:2013xwa (); Haskell:2013hja ().

It might seem that a no-r-mode scenario is also possible, since each point is just an upper limit, so in the - plot these stars could be at  K, outside the hadronic matter instability region, so their is really zero. To see that this is not possible, consider the evolutionary history of these sources: they are recycled pulsars, previously spun up by accretion in a LMXB Brown:1998ch (). During accretion they are heated to by nuclear reactions in the crust (left panel of fig. 1). If the hadronic matter model were correct, these stars would then become trapped in the instability region, since once the accretion stops the stars would quickly cool until they reach the spindown line (determined by ) at which cooling balances r-mode heating. They then spin down along the curve very slowly: at the frequency would change only by a few tens of Hz in a billion years. The only way that the fastest spinning pulsars could escape the instability region is if , allowing them to cool to  K in a few million years without crossing the r-mode spindown line. So this scenario requires an even lower saturation amplitude than the saturated r-mode scenario, much lower than is predicted by any proposed r-mode saturation mechanism. This account is not contradicted by the only ms-radio source J0437-4715, for which a temperature estimate is available Durant:2011je (). Its spin frequency is low enough that it could have cooled out of the instability region without crossing an r-mode spindown curve for any . The figures also show the first LMXB (IGR J00291+5934), whose (temporary) spindown during quiescence has been observed Patruno:2010qz ().

Finally consider the thermal state of old radio sources. In the no r-mode scenario they should have cooled to very low temperatures (potentially set by other heating sources Durant:2011je ()), whereas in the boundary straddling scenario the temperature would be independent of the spindown rate and determined by eq. (2). Yet, if the saturated r-mode scenario is realized then, on each side of fig. 2, ms-pulsars sit on the spindown curve for the physical value of . If the star is in steady state balance between r-mode heating and neutrino/photon cooling then its temperature is a function of its spin and r-mode spindown rate . Using eqs. (3) and (4) this r-mode steady-state temperature is


where is the moment of inertia. The striking feature of this simple expression is that it is independent of the saturation physics. This is because it is determined by rotational energy being transformed into gravitational wave, neutrino and/or photon energy, irrespective of the r-mode physics that accomplishes this. However, eq. (6) depends on the cooling behavior which differs for various forms of matter. The corresponding r-mode steady-state temperatures of radio pulsars that would be observed on earth are shown in fig. 3 for the two extreme cases of standard modified Urca cooling and fast direct Urca cooling. If the spindown of a star is dominated by r-mode gravitational emission then eq. (6) tells us its core temperature. If only a fraction of the spindown rate is due to r-modes then it provides only an upper bound. For the actual temperature to be significantly below eq. (6) would require that only a fraction of the observed spindown rate is due to r-modes, which would require an even smaller value of than the bounds obtained from fig. 2, and we do not know of any mechanism that could accomplish this. We conclude that if radio pulsars are, as the hadronic model requires, undergoing r-mode spindown, even tiny amplitude r-modes would have a big impact on the thermal evolution. In this case radio pulsars should have observable surface temperatures as measured at infinity of  K, which is significantly hotter than standard cooling estimates suggest Yakovlev:2004iq ().

Figure 3: The r-mode temperatures as observed at infinity of radio pulsars with known timing data as well as the corresponding temperatures of LMXBs (dots) compared to r-mode instability boundaries. R-modes temperatures are given for sources with assuming neutron stars with modified Urca (stars) and direct Urca cooling (diamonds).

We conclude that the novel dynamic instability regions and the r-mode temperature show us how timing data from the large population of radio pulsars can constrain their interior constitution. Beyond the illustrative examples discussed here, there are various options, like superfluid pairing (including the important effect of mutual friction) Haskell:2012 (), magnetic fields Mendell:2001nz (), hyperonic matter Reisenegger:2003cq () or color-superconducting phases Alford:2007xm (), that might be responsible for enhanced damping. To come to definite conclusions will require both observational and theoretical progress. On the theoretical side, we need saturation amplitudes, dynamic instability regions and r-mode temperatures for all hypothesized forms of dense matter with distinct damping and cooling properties. Observationally, it would be particularly useful to obtain temperature measurements or bounds for nearby ms-radio pulsars that spin with frequencies above . The comparison with the theoretical r-mode stability boundary could reveal whether the saturated or boundary-straddling r-mode scenarios can be realized. If they are so cool as to lie outside the boundary this would be inexplicable in the minimal hadronic model. This is just one example of how the combination of radio, x-ray and future gravitational wave data will allow us to discriminate the no-, saturated- or boundary-straddling r-mode scenarios and eventually different phases of dense matter.

We are grateful to Simin Mahmoodifar and Tod Strohmayer for helpful discussions. This research was supported by the Offices of Nuclear and High Energy Physics of the U.S. Department of Energy under contracts #DE-FG02-91ER40628, #DE-FG02-05ER41375.


  • (1) J. Lattimer and M. Prakash, Science 304, 536 (2004), astro-ph/0405262.
  • (2) N. Itoh, Prog. Theor. Phys. 44, 291 (1970).
  • (3) E. Witten, Phys. Rev. D30, 272 (1984).
  • (4) M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schafer, Rev. Mod. Phys. 80, 1455 (2008), 0709.4635.
  • (5) Fermi-LAT collaboration, A. Abdo et al., Astrophys.J.Suppl. 208, 17 (2013), 1305.4385.
  • (6) R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, Astron. J. 129, 1993 (2005), astro-ph/0412641.
  • (7) N. Andersson, Astrophys. J. 502, 708 (1998), gr-qc/9706075.
  • (8) N. Andersson and K. D. Kokkotas, Int. J. Mod. Phys. D10, 381 (2001), gr-qc/0010102.
  • (9) J. L. Friedman and B. F. Schutz, Astrophys. J. 222, 281 (1978).
  • (10) K. Schwenzer, (2012), 1212.5242.
  • (11) M. G. Alford and K. Schwenzer, Astrophys.J. 781, 26 (2014), 1210.6091.
  • (12) M. Alford, S. Mahmoodifar, and K. Schwenzer, Phys.Rev. D85, 024007 (2012), 1012.4883.
  • (13) L. Lindblom, B. J. Owen, and G. Ushomirsky, Phys.Rev. D62, 084030 (2000), astro-ph/0006242.
  • (14) B. Haskell, N. Degenaar, and W. C. G. Ho, Mon. Not. Roy. Astron. Soc. 424, 93 (2012), 1201.2101.
  • (15) E. H. Gudmundsson, C. J. Pethick, and R. I. Epstein, Astrophys.J. 272, 286 (1983).
  • (16) A. Reisenegger and A. A. Bonacic, Phys. Rev. Lett. 91, 201103 (2003), astro-ph/0303375.
  • (17) P. S. Shternin and D. G. Yakovlev, Phys. Rev. D78, 063006 (2008), 0808.2018.
  • (18) H. Heiselberg and C. J. Pethick, Phys. Rev. D48, 2916 (1993).
  • (19) R. F. Sawyer, Phys. Rev. D39, 3804 (1989).
  • (20) M. G. Alford, S. Mahmoodifar, and K. Schwenzer, J. Phys. G37, 125202 (2010), 1005.3769.
  • (21) J. Madsen, Phys. Rev. Lett. 81, 3311 (1998), astro-ph/9806032.
  • (22) J. Madsen, Phys. Rev. Lett. 85, 10 (2000), astro-ph/9912418.
  • (23) R. Bondarescu and I. Wasserman, Astrophys.J. 778, 9 (2013), 1305.2335.
  • (24) B. Haskell, K. Glampedakis, and N. Andersson, (2013), 1307.0985.
  • (25) L. Lindblom, J. E. Tohline, and M. Vallisneri, Phys. Rev. Lett. 86, 1152 (2001), astro-ph/0010653.
  • (26) J. A. Tomsick, D. M. Gelino, J. P. Halpern, and P. Kaaret, Astrophys.J. 610, 933 (2004), astro-ph/0404287.
  • (27) A. Y. Potekhin, G. Chabrier, and D. G. Yakovlev, Astron. Astrophys. 323, 415 (1997), arXiv:astro-ph/9706148.
  • (28) J. Madsen, Phys. Rev. D46, 3290 (1992).
  • (29) B. J. Owen et al., Phys. Rev. D58, 084020 (1998), gr-qc/9804044.
  • (30) S. Mahmoodifar and T. Strohmayer, Astrophys.J. 773, 140 (2013), 1302.1204.
  • (31) D. G. Yakovlev and C. J. Pethick, Ann. Rev. Astron. Astrophys. 42, 169 (2004), astro-ph/0402143.
  • (32) N. Andersson, D. I. Jones, and K. D. Kokkotas, Mon. Not. Roy. Astron. Soc. 337, 1224 (2002), astro-ph/0111582.
  • (33) E. F. Brown, L. Bildsten, and R. E. Rutledge, Astrophys.J. 504, L95 (1998), arXiv:astro-ph/9807179.
  • (34) M. Durant et al., Astrophys.J. 746, 6 (2012), 1111.2346.
  • (35) A. Patruno, Astrophys.J. 722, 909 (2010), 1006.0815.
  • (36) G. Mendell, Phys.Rev. D64, 044009 (2001), gr-qc/0102042.
  • (37) C. Alcock, E. Farhi, and A. Olinto, Astrophys.J. 310, 261 (1986).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description