Evidence for a minimum ellipticity in millisecond pulsars

Neutron stars spin down over time due to a number of energy-loss processes. We provide tantalizing population-based evidence that millisecond pulsars (MSPs) have a minimum ellipticity of $\epsilon\approx10^{-9}$ around their spin axis and that, consequently, some spin down mostly through gravitational-wave emission. We discuss the implications of such a minimum ellipticity in terms of the internal magnetic field strengths and nuclear matter composition of neutron stars and show it would result in the Advanced LIGO and Virgo gravitational-wave detectors, or their upgrades, detecting gravitational waves from some known MSPs in the near future.


INTRODUCTION
Any rotating system generates gravitational waves if its mass is not arranged symmetrically around the axis of rotation, and the recent discovery of gravitational waves from the coalescing binary neutron star system GW170817 is a clear example of this general result (Abbott et al. 2017a). Even isolated, but rapidly-rotating, neutron stars (if sufficiently asymmetric) could generate gravitational-wave signals detectable on Earth. Indeed, approximately 451 of the 2 636 known radio and X-ray pulsars would generate continuous quadrupolar (l = 2) emission that falls within the frequency range of ground-based gravitational-wave detectors, prompting deep searches for these signals in data from both LIGO and Virgo (e.g., Aasi et al. 2014;Abbott et al. 2017b). 1 To date no signals have been detected, but these investigations have placed stringent upper limits on the mass quadrupole of many known pulsars, and in several instances limits on the fraction of overall lugraham.woan@glasgow.ac.uk 1 These numbers assume all pulsars with rotational frequencies greater than 10 Hz are within the sensitive frequency range of gravitational-wave detectors. They are taken from version 1.58 of the Australia Telescope National Facility (ATNF) Pulsar Catalog (Manchester et al. 2005). minosity that can be attributed to gravitational-wave emission (Abbott et al. 2017b).
The strength of the gravitational-wave emission from a neutron star depends on both the degree of asymmetry and the rotation rate. For a non-precessing triaxial star, rotating about its z principal axis, the asymmetry is dominated by the m = 2 mass quadrupole Q 22 and characterised by the moment of inertia ellipticity (see, e.g., Owen 2005) where I ii are the principal moments of inertia and I xx I yy I zz = I. A star of rotational period P and ellipticity will generate gravitational waves of period P/2 and with a gravitational luminosity of L GW = 2048π 6 5 G c 5 I 2 2 P −6 10 29 10 −9 2 1 ms P 6 W.
(2) The luminosity is proportional to the square of the third time derivative of the (reduced) moment of inertia tensor, giving it a strong period dependence. This loss of energy acts as a rotational brake on the neutron star and, for a pulsar, contributes to an observed rate of change of period, or "spin-down" rate.
However, gravitational radiation is usually expected to be a relatively small contribution to the overall spindown. Pulsars are thought to have strong external magnetic fields, with spin-down rates dominated by magnetic dipole, rather than gravitational quadrupole, radiation and with additional braking mechanisms present including wind-induced mass-loss.
A pulsar spinning at an angular frequency ω has a braking index n that satisfiesω ∝ −ω n (orṖ ∝ P −(n−2) ), and n can be measured for some pulsars, revealing a range of values (see, e.g., Table 1 of Lyne et al. 2015 and references therein).
The process of spin-down is clearly complicated, but an isolated, magnetically braked, rigid rotator would have a braking index of three. If gravitational emission is the dominant process we have a "gravitar" with a braking index of five (Palomba 2005).
Right from the start of pulsar studies, it was noted by Ferrari & Ruffini (1969) that the overall spin-down in even the simplest systems would have contributions from both, so that at the very leasṫ LIGO and Virgo observations have shown the gravitational contribution to the spin-down of the Crab pulsar to be tiny (Abbott et al. 2017b), but its spindown rate is high and it would need a possibly unphysical ∼ 10 −4 to be dominated by gravitationalwave emission. Conversely, gravitational observations have constrained the ellipticity of some millisecond pulsars (MSPs) to ∼ 10 −7 or less (Abbott et al. 2017b). Such ellipticities are well within the bounds set by likely neutron star equations of state (Owen 2005;Johnson-McDaniel & Owen 2013), and the short rotational period of MSPs would make them relatively luminous gravitational wave sources. The spin periods (P ) and spin-down rates "Pdot" (Ṗ ) of the pulsar population are traditionally displayed on a P-Pdot diagram (e.g., Figure 1.13 of Lorimer & Kramer 2004) showing how known pulsars cluster in particular regions of this parameter space, and new surveys have bolstered the number of sources considerably. The clustering is thought to result from a mixture of observational selection effects and underlying physics.
In this Letter, we discuss whether there is evidence for a new cutoff emerging in the diagram, at short-period and low-spin-down rate, caused by gravitational-wave emission and consistent with a minimum ellipticity for MSPs. Such a cutoff would correspond to a population of rapidly-rotating gravitars, sufficiently luminous to be detectable by current or future ground-based gravitational observatories.

DATA AND MODEL
To investigate this apparent cutoff we will consider known pulsars with periods P < 10 ms and period derivativesṖ < 10 −18 s s −1 . The current ATNF pulsar catalogue (Manchester et al. 2005) contains 199 pulsars that fulfil these criteria. The observed values of their period derivatives will, to varying degrees, be contaminated by radial accelerations due to proper motion (the Shklovskii effect), differential Galactic rotation (Damour & Taylor 1991) and, for pulsars in globular clusters, local forces (Freire et al. 2017). We therefore need to consider these effects and, as far as possible, work with the true (i.e., intrinsic) period derivatives.
We begin by excluding all 59 globular cluster pulsars from our sample. We also exclude PSR J1801−3210, which shows no measurable period derivative but which is thought to be affected by Galactic acceleration (Ng et al. 2014), and J1400−1431, for which there is only an upper limit on Pdot (Swiggum et al. 2017); this leaves 128 MSPs. For 28 of these we use intrinsic period derivatives from the literature, already corrected for Shklovskii and differential Galactic rotation effects using parallaxbased distance estimates (and in many cases also corrected for Lutz-Kelker bias) (Deller et al. 2012;Ng et al. 2014;Desvignes et al. 2016;Spiewak et al. 2018). For the remaining pulsars we calculate the corrections using the model of Damour and Taylor but with a Galactic radius of 8.3 kpc. We use parallax-derived distances when available, either from the literature (corrected for Lutz-Kelker bias; Desvignes et al. 2016;Reardon et al. 2016) or using the values given by the ATNF pulsar catalog (Manchester et al. 2005). If no parallax distance is known we use the best-estimate distance given in the ATNF catalogue, which by default uses the measured dispersion measure and the Galactic electron density model of Yao et al. (2017). The P-Pdot diagrams using intrinsic (circles) and observed (stars) period derivatives are shown in Figure 1.
We can cast Equation (3) into a standard form by assuming a neutron star of radius R, with magnetic spindown due to vacuum dipole radiation and surface magnetic field intensity B s : Using canonical values for the radius (10 km) and moment of inertia (10 38 kg m 2 ) of the neutron star this becomes Ṗ 10 −20 s/s =0.98 1 ms P B s 10 8 Gauss 2 + 2.7 1 ms P 3 10 −9 2 .
Canonical gravitars (i.e., neutron stars obeying Equation (5) but with negligible magnetic field) would fall on the straight orange lines in Figure 1. The blue curves show where canonical pulsars with = 10 −9 would be located for different values of surface magnetic field. Pulsars with a range of ellipticities and magnetic fields would also fall in this region, consistent with Equation (5). The number of pulsars involved is small, but there is some evidence for a cutoff in the population below the gravitar = 10 −9 line in the limit of low magnetic field, consistent with the notion that MSPs have a residual ellipticity that does not tend to fall below this level. Although such weakly magnetized pulsars would constantly radiate the equivalent of several solar-luminosities in gravitational waves, their reservoir of rotational kinetic energy is huge (∼ 10 43 J) and their gravitational spin-down age, P/(4Ṗ ), is 10 8 -10 10 years. These limiting neutron stars must still have sufficient magnetic field to be seen as pulsars, but their dominant spin-down mechanism would be gravitational. The two pulsars below the = 10 −9 line are PSR J2322−2650 (Spiewak et al. 2018), a recently discovered MSP with low radio luminosity and a low-density planetary-mass companion, and PSR J1017−7156 (Ng et al. 2014).

Fitting the P-Pdot distribution to the model
The model described by Equation (5) assumes the same moment of inertia for all neutron stars. In reality, the moment of inertia of a star depends on its mass via the equation of state, and particularly massive or light neutron stars may have larger (or smaller respectively) values of Pdot than those predicted by Equation (5) of up to a factor ≈ 2 (Worley et al. 2008), leading to additional scatter in the P-Pdot diagram. The significance of the apparent cutoff therefore needs to be considered more carefully, with these effects in mind.
We therefore construct a simple model of the distribution of pulsars over the P-Pdot plane. We assume a priori that spin-down rates are distributed uniformly in log-space, with a lower cutoff that follows a power law of the formṖ = kP −(n−2) , corresponding to a braking index of n. A special case of this cutoff process would be gravitational radiation (Equation (4) with B s = 0) from a common ellipticity. The lower-right corner of the P-Pdot plane is also largely free of sources and is delimited by the "death line," below which neutron stars are not observable as pulsars. We use the death line defined in Equation (3) of Zhang et al. (2000) to exclude the region belowṖ = L d (P/1 s) 11/4 s/s, where L d = 10 −14.62 , and find the common cutoff process that best explains the observed distribution. There are also few MSPs witḣ P 10 −19 s s −1 , which will depress the model evidence values but not significantly affect comparisons between parameterizations of the lower cutoff.
We assign a Gaussian likelihood for the true period derivative of the jth pulsarṖ j , given our measurement of its intrinsic Pdot and its uncertainty, so that where the values of d j ≡ {Ṗ int j , σ j } are set by the intrinsic Pdot estimation procedure described in Section 2, and shown in Figure 1. If uncertainties for the intrinsic Pdot are available from the literature, then we use those (see Section 2). Otherwise we combine, in quadrature, the measurement uncertainties on Pdot, the Shklovskii correction, and the Galactic correction, and these comprise the error bars in Figure 1. We assign a conservative 50% uncertainty on any distances derived from dispersion measure. We incorporate the cutoff process and death line into the log-uniform prior onṖ j as a common thresholdṖ th j (P j ), such that whereṖ max is the maximum value ofṖ int in our sample.
The exclusion region delimited byṖ th (P ) depends on both our new cutoff line and the death line, and we takė P th (P ) to be the maximum ofṖ (P ) defined by the cutoff line (parameterized by k and n) and the death line and incorporate the moment of inertia of each pulsar into the cutoff process by setting k j = I j k . We can then set a Gaussian prior on I j , with a common mean µ I and standard deviation σ I over all pulsars and choose µ I = 2×10 −38 kg m 2 and σ I = 3×10 37 kg m 2 using the ranges shown in Figures 4 & 7 of Worley et al. (2008), with hard bounds of [1×10 38 , 3×10 38 ] kg m 2 . We also incorporate uncertainty in the position of the death line by allowing the coefficient L d to be a free parameter with a prior that is uniform in log-space between [10 −15.62 , 10 −14.62 ], which spans most of the range of death lines in Zhang et al. (2000). Therefore, for a given pulsar,Ṗ th j is defined by P j , k , I j , n and L d .
For each pulsar, we first marginalize overṖ j and I j to give We then form a joint likelihood over all N pulsars in our sample and marginalize over L d and k to give the marginal likelihood, or evidence, for the model as a function of the cutoff process braking index n, giving where d is the combined data {d i }, and the prior on k , p(k ), is log-uniform over a range for which ln (k max /k min ) ≈ 44.
Using all the data from pulsars for which we have confident measurements of intrinsic Pdot (shown in Figure 1) we can therefore determine Bayes factors, comparing models with and without a common cutoff process other than the death line. If we take the cutoff process to be due to an n = 5 process, then we find that B is hugely in favour of a distribution with a nonzero ellipticity (a factor of ∼ 6 400), corresponding to an ellipticity of ∼ 5.3 +0.4 −0.7 × 10 −10 . We note here that as we have taken a distribution for I centred at twice the canonical value, this ellipticity corresponds to the same mass quadrupole, and therefore gravitational wave amplitude, as choosing an ellipticity of ∼ 10 −9 and using the canonical moment of inertia.
Comparing the evidence for an n = 5 process to an n = 3 process we find that the former is favoured by a factor ∼ 35. In fact, we find that a process with n ≈ 5.6 +1.1 −0.4 is the most favoured case, being ∼ 10, 300 times more likely than having no cutoff, and ∼ 55 times more likely than an n = 3 process. We conclude, therefore, that the observed fall-off in the number of sources at the bottom-left of Figure 1 is highly significant and much more likely due to a common ellipticity amongst MSPs than due to a common minimum surface magnetic field. Although this simple model explains the shortperiod/low Pdot cutoff nicely, there could be other processes and selection effects involved that we are not modeling. Very low intrinsic period derivatives are difficult to measure, but the most important corrections (applied above) are thought to be well-understood. As is clear in the case of globular cluster pulsars, significant but uncorrected additive errors in Pdot tend to scatter pulsars well away from the small region of (linear) parameter space containing the cutoff boundary, and such pulsars would not influence the analysis. It should also be noted that a relatively long time baseline is required to measure low Pdot values, but this becomes somewhat easier for shorter period pulsars, so these should not be preferentially excluded by this requirement.
Although our Bayes factors are naturally modulated by prior assumptions, they strongly support the apparent power-law cutoff one sees in the P-Pdot plot, and show it to be consistent with a limiting braking index of n = 5, rather than n = 3. Evidence for or against a simple cutoff would be strengthened by a larger sample of pulsars with periods below ∼ 3 ms, where the death line has less of an influence. Similarly, a larger number of pulsars close to the cutoff line might indicate a form more complex than a power law. However, with the current limited sample size, more complex models are naturally disfavored by the Occam factor.

DISCUSSION
It is reasonable to ask what could cause the minimal ellipticity described in Section 2. Magnetic field evolution in MSPs is not well-understood; however, these stars have relatively small external dipole magnetic fields (B s ∼ 10 8 G) compared to the younger population of radio pulsars (10 11 B s /G 10 13 ).
One possibility is that the external magnetic field becomes buried while the system is undergoing accretion from the binary companion. If this is the case, one may expect an internal magnetic field on the order of 10 11 G (Vigelius & Melatos 2009). MSPs are old, cold stars, and the protons in their cores are expected to form a type II superconductor. For such stars, the ellipticity is linear in the internal magnetic field strength. The exact value is model-dependent, but generally of the order of (Lander et al. 2012;Lander 2014) where H c is the lower critical field for superconductivity. An ellipticity of 10 −9 would therefore be consistent with a buried field of B ∼ 10 11 G, itself consistent with the field strengths observed in the general pulsar population. Another possible explanation for such a minimum ellipticity could be strain in the elastic crust; a value of ≈ 10 −9 can be accommodated with dimensionless strain levels of around 10 −4 (Ushomirsky et al. 2000), comfortably below the breaking strains indicated by molecular simulations (Horowitz & Kadau 2009). Such strains might be produced by asymmetries in the accretion process through which these MSPs were originally spun-up (Bildsten 1998;Ushomirsky et al. 2000), or else by asymmetric crustal fracture during the accretion spin-up (Fattoyev et al. 2018) or the subsequent post-accretion spin-down (Baym & Pines 1971).
At this point it is interesting to ask whether there is any evidence for such a minimum ellipticity in the population of accreting millisecond X-ray pulsars in low-mass X-ray binaries, which are thought to be the progenitors of millisecond radio pulsars. They have also long been studied as a source of gravitational waves because their spin periods appear to cluster around P ≈ 2 ms -well below the theoretical Keplerian breakup frequency to which a neutron star could be spun-up . Gravitational waves could be providing an additional spin-down torque that balances the spin-up torque due to accretion. This leads to equilibrium at the observed periods for ellipticities of = 10 −8 in most of the population (Bildsten 1998;Ushomirsky et al. 2000). There are however, two systems close to the minimum observed period that are transitional pulsars; i.e., are making the transition from being accretion powered Xray pulsars to being rotationally powered millisecond radio pulsars: J1023+0038 and J1227−4853. Both are accreting at low rates, so that the accretion torque would be balanced for ≈ 2 × 10 −10 . Nevertheless, these systems cannot be in spinequilibrium as they are observed to spin-down in radio, and J1023+0038 not only appears in our sample close to the cutoff, but it is also observed to spin-down ≈ 30% faster in X-ray (Jaodand et al. 2016). Intriguingly, this effect could be explained by gravitational-wave emission , suggesting that these systems may be the progenitors of the pulsars close to the cutoff in the P-Pdot diagram.

IMPLICATIONS FOR FUTURE GRAVITATIONAL-WAVE DETECTIONS
We now determine the corresponding gravitationalwave signal strength from our selected MSPs (using, e.g., Equation (3) of Aasi et al. 2014) under two models: (i) that they each have an ellipticity of 10 −9 as considered above and (ii) that they are all maximal gravitars, with signal strengths set by their spin-down limits.
In the first case we are able to re-include those pulsars in globular clusters, as knowledge of their true period derivatives is not required. We use the canonical moment of inertia and distances obtained as described in Section 2, keeping in mind that there could be factortwo uncertainties in both these values. In the gravitar case we calculate spin-down limits only for those pulsars not associated with globular clusters and for which we have confident measurements of their intrinsic Pdot, as discussed in Section 2 and shown in Figure 1.
Under these models we have estimated the signal-tonoise ratios (see Equation (2) of Pitkin 2011) that pulsars with rotation periods less than 10 ms would have in future networks of gravitational-wave detectors, assuming fully-coherent targeted searches. Their signalto-noise ratios depend on the (unknown) inclination of each pulsar's rotation axis with respect to the line-ofsight. We have calculated signal-to-noise ratios for the angle-averaged case, a factor of ∼ 1.69 times below the best-case, assuming a uniform prior on orientation (Pitkin 2011). The networks we have considered are: (a) two advanced LIGO (aLIGO) detectors, in Hanford and Livingston, operating at design sensitivity (LIGO Scientific Collaboration 2010; Aasi et al. 2015) with the advanced Virgo (AdV) detector also operating at design sensitivity (Virgo Collaboration 2009;Acernese et al. 2015); (b) this configuration combined with an equivalent LIGO detector in India (Unnikrishnan 2013) and the KAGRA detector (Aso et al. 2013) at design sensitivity; (c) two upgraded aLIGO+ (LIGO Scientific Collaboration 2018) detectors, together with AdV; (d) Cosmic Explorer (Abbott et al. 2017c); (e) the Einstein Telescope in its ET-D configuration (Hild et al. 2011;Sathyaprakash et al. 2012), assuming three co-located detectors. In all cases we use a one-year observation pe- riod and coherently combine data from all detectors in a network. The distributions of signal-to-noise ratios for all these scenarios are shown in Figure 2, where the filled histograms represent sources all having ellipticities of 10 −9 and the canonical moment of inertia, and the unfilled histograms show the distribution when the MSPs are gravitars. If our minimum-ellipticity hypothesis is correct, we would expect the true distribution to be somewhere between these two extremes, modulo the uncertainties in distance and moment of inertia. The pulsar that would give the highest signal-to-noise ratio with an ellipticity of 10 −9 in all these scenarios is PSR J1643−1224, at a distance of 0.79 kpc and rotating at 216.4 Hz. The pulsar that would have the highest signal-to-noise ratio as a pure gravitar is PSR J0711−6830, with a rotation frequency of 182.1 Hz and a distance of 0.11 kpc. Figure 2 indicates that marginal detections from this MSP population could be possible with the network of second-generation detectors, and more confident detections possible with future third-generation detectors.
It is worth noting that, on the time scale of upgrades to aLIGO, the Square Kilometre Array plans to start operation, bringing with it the possibility of finding ∼ 6 000 more MSPs (Smits et al. 2009). Pulsars that occupy the current gap at the bottom-left of the P-Pdot diagram, i.e., with short periods and very low spin-down rates, would provide strong evidence against our hypothesis. Of course the hypothesis would be strongly supported if the braking index of a boundary MSP was measured to be close to n = 5. However many decades or even centuries of observations, would most likely be required to confidently determine the second period derivative necessary for this.
We thank the Institute for Nuclear Theory at the University of Washington, the Department of Energy, and the organizers of the "Astro-Solids, Dense Matter, and Gravitational Waves" workshop for their hospitality and for partial support during completion of this work. M.P. is funded by the UK Science & Technology Facilities Council (STFC) under grant ST/N005422/1. P.D.L. is supported through Australian Research Council (ARC) Future Fellowship FT160100112 and ARC Discovery Project DP180103155.
B.H. acknowledges support from the Polish National Science Centre (SONATA BIS 2015/18/E/ST9/00577) and from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 702713. D.I.J. acknowledges funding from STFC through grant No. ST/M000931/1. Partial support comes from PHAROS, COST Action CA16214. The work has been made possible using data from the ATNF Pulsar Catalog (Manchester et al. 2005). This work has been assigned LIGO document number LIGO-P1800125 and INT preprint number INT-PUB-18-031.