#  Evidence for a minimum ellipticity in millisecond pulsars

G. Woan  M. D. Pitkin  SUPA, School of Physics & Astronomy, University of
Glasgow, Glasgow G12 8QQ, United Kingdom<EMAIL_ADDRESS>B. Haskell
Nicolaus Copernicus Astronomical Center of the Polish Academy of Sciences,
Ulica Bartycka 18, 00-716 Warszawa, Poland  D. I. Jones  Mathematical Sciences
and STAG Research Centre, University of Southampton, Southampton SO17 1BJ,
United Kingdom  P. D. Lasky  School of Physics and Astronomy, Monash
University, Vic 3800, Australia  OzGrav: The ARC Centre of Excellence for
Gravitational Wave Discovery, Monash University, Vic 3800, Australia

###### Abstract

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\approx 10^{-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.

pulsars: general — stars: neutron — stars: magnetic field — gravitational
waves

††software: Jupyter (Kluyver et al., 2016), NumPy (Oliphant, 2006), SciPy
(Jones et al., 2001–), Matplotlib (Hunter, 2007), PSRQpy (Pitkin, 2018)

arXiv:1806.02822 LIGO-P1800125 INT-PUB-18-031

##  1 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).111These 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). 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 luminosity 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 $\epsilon$ (see, e.g., Owen, 2005)

$\epsilon=\frac{I_{xx}-I_{yy}}{I_{zz}}=\frac{Q_{22}}{I}\sqrt{\frac{8\pi}{15}},$
(1)

where $I_{ii}$ are the principal moments of inertia and $I_{xx}\simeq
I_{yy}\simeq I_{zz}=I$. A star of rotational period $P$ and ellipticity
$\epsilon$ will generate gravitational waves of period $P/2$ and with a
gravitational luminosity of

$L_{\text{GW}}=\frac{2048\pi^{6}}{5}\frac{G}{c^{5}}I^{2}\epsilon^{2}P^{-6}\simeq
10^{29}\left(\frac{\epsilon}{10^{-9}}\right)^{2}\left(\frac{1\,\text{ms}}{P}\right)^{6}\\!\\!\text{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 spin-down. 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 $\omega$ has a braking index $n$
that satisfies $\dot{\omega}\propto-\omega^{n}$ (or $\dot{P}\propto
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 least

$\dot{\omega}=-\alpha\omega^{3}-\beta\omega^{5}.$ (3)

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 spin-
down rate is high and it would need a possibly unphysical $\epsilon\sim
10^{-4}$ to be dominated by gravitational-wave emission. Conversely,
gravitational observations have constrained the ellipticity of some
millisecond pulsars (MSPs) to $\epsilon\sim 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” ($\dot{P}$) 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.

##  2 Data and model

To investigate this apparent cutoff we will consider known pulsars with
periods $P<10$ ms and period derivatives $\dot{P}<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 parallax-based 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 spin-down due to vacuum dipole radiation and surface
magnetic field intensity $B_{\text{s}}$:

$\dot{P}=\frac{32\pi^{3}R^{6}}{3Ic^{3}\mu_{0}}B_{\text{s}}^{2}P^{-1}+\frac{512\pi^{4}GI}{5c^{5}}\epsilon^{2}P^{-3}.$
(4)

Using canonical values for the radius (10 km) and moment of inertia ($10^{38}$
kg m2) of the neutron star this becomes

$\begin{split}\left(\frac{\dot{P}}{10^{-20}\,\text{s/s}}\right)=&0.98\left(\frac{1\,\text{ms}}{P}\right)\left(\frac{B_{\text{s}}}{10^{8}\,\text{Gauss}}\right)^{2}\\\
&+2.7\left(\frac{1\,\text{ms}}{P}\right)^{3}\left(\frac{\epsilon}{10^{-9}}\right)^{2}.\end{split}$
(5)

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 $\epsilon=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
$\epsilon=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 ($\sim 10^{43}$ J) and
their gravitational spin-down age, $P/(4\dot{P})$, 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 $\epsilon=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).

Figure 1:  Observed (black stars) and intrinsic (blue circles) period
derivatives of the MSPs in our sample vs. period (i.e., the bottom-left corner
of the standard P-Pdot diagram) excluding pulsars in globular clusters. The
intrinsic period derivatives include corrections for differential Galactic
rotation and Shklovskii effects. Straight (orange) lines show where canonical
gravitars with ellipticities of $10^{-8}$ and $10^{-9}$ would sit. Note the
sudden fall in source numbers below $\epsilon\simeq 10^{-9}$. Also shown are
curves (blue) of constant surface magnetic field $B_{\text{s}}$, assuming
spin-down is in accordance with Equation (5) and that all the pulsars have a
common ellipticity of $10^{-9}$. The shaded area corresponds to the death line
exclusion region of Zhang et al. (2000).

###  2.1 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 $\approx 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 $\dot{P}=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_{\rm 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 $\dot{P}=L_{\text{d}}(P/1\,\text{s})^{11/4}$ s/s, where
$L_{\text{d}}=10^{-14.62}$, and find the common cutoff process that best
explains the observed distribution. There are also few MSPs with
$\dot{P}\gtrapprox 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 $j$th
pulsar $\dot{P}_{j}$, given our measurement of its intrinsic Pdot and its
uncertainty, so that

$p(d_{j}|\dot{P}_{j})\propto\exp\left[-\frac{1}{2\sigma_{j}^{2}}\left(\dot{P}_{j}-\dot{P}^{\text{int}}_{j}\right)^{2}\right],$
(6)

where the values of $d_{j}\equiv\\{\dot{P}^{\text{int}}_{j},\sigma_{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 $\dot{P}_{j}$ as a common threshold $\dot{P}_{j}^{\text{th}}(P_{j})$,
such that

$p(\dot{P}_{j}|\dot{P}_{j}^{\rm
th})=\begin{cases}\left[\dot{P}_{j}\ln\frac{{\dot{P}^{\rm
max}}}{\dot{P}^{\text{th}}_{j}}\right]^{-1}&\text{if }\dot{P}_{j}^{\rm
th}\leqslant\dot{P}_{j}\leqslant\dot{P}^{\rm max},\\\
0&\text{otherwise},\end{cases}$ (7)

where $\dot{P}^{\text{max}}$ is the maximum value of $\dot{P}^{\text{int}}$ in
our sample.

The exclusion region delimited by $\dot{P}^{\text{th}}(P)$ depends on both our
new cutoff line and the death line, and we take $\dot{P}^{\rm th}(P)$ to be
the maximum of $\dot{P}(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^{\prime}$. We can then
set a Gaussian prior on $I_{j}$, with a common mean $\mu_{I}$ and standard
deviation $\sigma_{I}$ over all pulsars and choose
$\mu_{I}=2\\!\times\\!10^{-38}\,{\rm kg}\,{\rm m}^{2}$ and
$\sigma_{I}=3\\!\times\\!10^{37}\,{\rm kg}\,{\rm m}^{2}$ using the ranges
shown in Figures 4 & 7 of Worley et al. (2008), with hard bounds of
$[1\\!\times\\!10^{38},3\\!\times\\!10^{38}]\,{\rm kg}\,{\rm m}^{2}$. We also
incorporate uncertainty in the position of the death line by allowing the
coefficient $L_{\text{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,
$\dot{P}_{j}^{\text{th}}$ is defined by $P_{j}$, $k^{\prime}$, $I_{j}$, $n$
and $L_{\text{d}}$.

For each pulsar, we first marginalize over $\dot{P}_{j}$ and $I_{j}$ to give

$p(d_{i}|L_{\text{d}},k^{\prime},n)=\iint
p(d_{i}|\dot{P}_{j})p(\dot{P}_{j}|\dot{P}^{\text{th}}_{j})p(I_{j})\,\text{d}\dot{P}_{j}\,\text{d}I_{j}.$
(8)

We then form a joint likelihood over all $N$ pulsars in our sample and
marginalize over $L_{\text{d}}$ and $k^{\prime}$ to give the marginal
likelihood, or evidence, for the model as a function of the cutoff process
braking index $n$, giving

$p(\bm{d}|n)=\iint\prod_{j=1}^{N}p(d_{j}|L_{\text{d}},k^{\prime},n)p(L_{\text{d}})p(k^{\prime})\,\text{d}L_{\text{d}}\,\text{d}k^{\prime},$
(9)

where $\bm{d}$ is the combined data $\\{d_{i}\\}$, and the prior on
$k^{\prime}$, $p(k^{\prime})$, is log-uniform over a range for which
$\ln{(k^{\prime}_{\text{max}}/k^{\prime}_{\text{min}})}\approx 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,

${\cal B}(n)=\frac{p(\bm{d}|n)}{p(\bm{d}|k^{\prime}=0)},$ (10)

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 $\cal B$ is hugely in favour of a distribution with a non-zero
ellipticity (a factor of $\sim 6\,400$), corresponding to an ellipticity of
$\sim 5.3^{+0.4}_{-0.7}\\!\times\\!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 $\sim 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 $\sim 35$. In fact, we find that a process
with $n\approx 5.6_{-0.4}^{+1.1}$ is the most favoured case, being $\sim
10,300$ times more likely than having no cutoff, and $\sim 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 short-period/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 $\sim 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.

##  3 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_{\text{s}}\sim 10^{8}$ G) compared to the younger population of radio
pulsars ($10^{11}\lesssim B_{\text{s}}/{\rm G}\lesssim 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)

$\epsilon\sim 10^{-8}\left(\frac{\langle
B_{\text{i}}\rangle}{10^{12}\,\text{G}}\right)\left(\frac{\langle
H_{\text{c}}\rangle}{10^{15}\,\text{G}}\right),$ (11)

where $H_{\text{c}}$ is the lower critical field for superconductivity. An
ellipticity of $10^{-9}$ would therefore be consistent with a buried field of
$B\sim 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 $\epsilon\approx 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\approx 2$ ms — well below the theoretical Keplerian breakup frequency to
which a neutron star could be spun-up (Patruno et al., 2017). 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 $\epsilon=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 X-ray 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 $\epsilon\approx 2\times 10^{-10}$ (Patruno et al., 2017).
Nevertheless, these systems cannot be in spin-equilibrium 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 $\approx 30\%$ faster in
X-ray (Jaodand et al., 2016). Intriguingly, this effect could be explained by
gravitational-wave emission (Haskell & Patruno, 2017), suggesting that these
systems may be the progenitors of the pulsars close to the cutoff in the
P-Pdot diagram.

##  4 Implications for future gravitational-wave detections

We now determine the corresponding gravitational-wave 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 factor-two
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-to-noise 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 signal-to-noise ratios depend on the
(unknown) inclination of each pulsar’s rotation axis with respect to the line-
of-sight. We have calculated signal-to-noise ratios for the angle-averaged
case, a factor of $\sim 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 period and
coherently combine data from all detectors in a network.

Figure 2: TDistribution in signal-to-noise ratios of a selection of MSPs,
using a variety of future gravitational-wave observatory networks with one
year of coherent observations. The shaded histograms give the distribution if
we assume that all MSPs have ellipticities of $10^{-9}$ and the canonical
moment of inertia, and the unfilled histograms are the values if they are all
gravitars emitting at their spin-down limits.

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 $\sim 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.

## References

  * Aasi et al. (2014) Aasi, J., Abadie, J., Abbott, B. P., et al. 2014, ApJ, 785, 119, doi: 10.1088/0004-637X/785/2/119
  * Aasi et al. (2015) Aasi, J., Abbott, B. P., Abbott, R., et al. 2015, CQGra, 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, PhRvL, 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, ApJ, 839, 12, doi: 10.3847/1538-4357/aa677f
  * Abbott et al. (2017c) —. 2017c, CQGra, 34, 044001, doi: 10.1088/1361-6382/aa51f4
  * Acernese et al. (2015) Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, CQGra, 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Aso et al. (2013) Aso, Y., Michimura, Y., Somiya, K., et al. 2013, PhRvD, 88, 043007, doi: 10.1103/PhysRevD.88.043007
  * Baym & Pines (1971) Baym, G., & Pines, D. 1971, AnPhy, 66, 816, doi: 10.1016/0003-4916(71)90084-4
  * Bildsten (1998) Bildsten, L. 1998, ApJ, 501, L89, doi: 10.1086/311440
  * Damour & Taylor (1991) Damour, T., & Taylor, J. H. 1991, ApJ, 366, 501, doi: 10.1086/169585
  * Deller et al. (2012) Deller, A. T., Archibald, A. M., Brisken, W. F., et al. 2012, ApJ, 756, L25, doi: 10.1088/2041-8205/756/2/L25
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341, doi: 10.1093/mnras/stw483
  * Fattoyev et al. (2018) Fattoyev, F. J., Horowitz, C. J., & Lu, H. 2018, arXiv:1804.04952
  * Ferrari & Ruffini (1969) Ferrari, A., & Ruffini, R. 1969, ApJ, 158, L71, doi: 10.1086/180435
  * Freire et al. (2017) Freire, P. C. C., Ridolfi, A., Kramer, M., et al. 2017, MNRAS, 471, 857, doi: 10.1093/mnras/stx1533
  * Haskell & Patruno (2017) Haskell, B., & Patruno, A. 2017, PhRvL, 119, 161103, doi: 10.1103/PhysRevLett.119.161103
  * Hild et al. (2011) Hild, S., Abernathy, M., Acernese, F., et al. 2011, CQGra, 28, 094013, doi: 10.1088/0264-9381/28/9/094013
  * Horowitz & Kadau (2009) Horowitz, C. J., & Kadau, K. 2009, PhRvL, 102, 191102, doi: 10.1103/PhysRevLett.102.191102
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jaodand et al. (2016) Jaodand, A., Archibald, A. M., Hessels, J. W. T., et al. 2016, ApJ, 830, 122, doi: 10.3847/0004-637X/830/2/122
  * Johnson-McDaniel & Owen (2013) Johnson-McDaniel, N. K., & Owen, B. J. 2013, PhRvD, 88, 044004, doi: 10.1103/PhysRevD.88.044004
  * Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas: Proceedings of the 20th International Conference on Electronic Publishing, IOS Press, Amsterdam, 87
  * Lander (2014) Lander, S. K. 2014, MNRAS, 437, 424, doi: 10.1093/mnras/stt1894
  * Lander et al. (2012) Lander, S. K., Andersson, N., & Glampedakis, K. 2012, MNRAS, 419, 732, doi: 10.1111/j.1365-2966.2011.19720.x
  * LIGO Scientific Collaboration (2010) LIGO Scientific Collaboration. 2010, Advanced LIGO anticipated sensitivity curves, Tech. Rep. LIGO-T0900288-v3
  * LIGO Scientific Collaboration (2018) —. 2018, The A+ design curve, Tech. Rep. LIGO-T180042-v5, LIGO Project
  * Lorimer & Kramer (2004) Lorimer, D. R., & Kramer, M. 2004, Cambridge Observing Handbooks for Research Astronomers, Vol. 4, Handbook of Pulsar Astronomy (Cambridge, UK: Cambridge University Press)
  * Lyne et al. (2015) Lyne, A. G., Jordan, C. A., Graham-Smith, F., et al. 2015, MNRAS, 446, 857, doi: 10.1093/mnras/stu2118
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993, doi: 10.1086/428488
  * Ng et al. (2014) Ng, C., Bailes, M., Bates, S. D., et al. 2014, MNRAS, 439, 1865, doi: 10.1093/mnras/stu067
  * Oliphant (2006) Oliphant, T. E. 2006, A Guide to NumPy (USA: Trelgol Publishing)
  * Owen (2005) Owen, B. J. 2005, PhRvL, 95, 211101, doi: 10.1103/PhysRevLett.95.211101
  * Palomba (2005) Palomba, C. 2005, MNRAS, 359, 1150, doi: 10.1111/j.1365-2966.2005.08975.x
  * Patruno et al. (2017) Patruno, A., Haskell, B., & Andersson, N. 2017, ApJ, 850, 106, doi: 10.3847/1538-4357/aa927a
  * Pitkin (2011) Pitkin, M. 2011, MNRAS, 415, 1849, doi: 10.1111/j.1365-2966.2011.18818.x
  * Pitkin (2018) —. 2018, JOSS, 3, 538, doi: 10.21105/joss.00538
  * Reardon et al. (2016) Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751, doi: 10.1093/mnras/stv2395
  * Sathyaprakash et al. (2012) Sathyaprakash, B., Abernathy, M., Acernese, F., et al. 2012, CQGra, 29, 124013, doi: 10.1088/0264-9381/29/12/124013
  * Smits et al. (2009) Smits, R., Kramer, M., Stappers, B., et al. 2009, A&A, 493, 1161, doi: 10.1051/0004-6361:200810383
  * Spiewak et al. (2018) Spiewak, R., Bailes, M., Barr, E. D., et al. 2018, MNRAS, 475, 469, doi: 10.1093/mnras/stx3157
  * Swiggum et al. (2017) Swiggum, J. K., Kaplan, D. L., McLaughlin, M. A., et al. 2017, ApJ, 847, 25, doi: 10.3847/1538-4357/aa8994
  * Unnikrishnan (2013) Unnikrishnan, C. S. 2013, IJMPD, 22, 1341010, doi: 10.1142/S0218271813410101
  * Ushomirsky et al. (2000) Ushomirsky, G., Cutler, C., & Bildsten, L. 2000, MNRAS, 319, 902, doi: 10.1046/j.1365-8711.2000.03938.x
  * Vigelius & Melatos (2009) Vigelius, M., & Melatos, A. 2009, MNRAS, 395, 1985, doi: 10.1111/j.1365-2966.2009.14698.x
  * Virgo Collaboration (2009) Virgo Collaboration. 2009, Advanced Virgo Baseline Design, Tech. Rep. VIR-027A-09
  * Worley et al. (2008) Worley, A., Krastev, P. G., & Li, B.-A. 2008, ApJ, 685, 390, doi: 10.1086/589823
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * Zhang et al. (2000) Zhang, B., Harding, A. K., & Muslimov, A. G. 2000, ApJ, 531, L135, doi: 10.1086/312542


This would rule out a number of interesting models and so it seems
unsatisfactory. An alternative is to define the iteration operator by:
\[
\begin{array}{rcl}
\textbf{iter}\ \A{Z}\ z\ s &\reducesTo & z\\
\textbf{iter}\ (\A{S}\ n)\ z\ s &\reducesTo & (\textbf{iter}\ n\ z\ s)
+ \epsilon(s\ (\textbf{iter}\ n\ z\ s))
\end{array}
\]
Fixed $z, s$, and defining the map $\mu(n) \defeq \textbf{iter}\ n\ z\ s$, its
derivative $\D[\mu](n, \A{S}\ \A{Z})$ is precisely $s(\mu(n))$. Or, in other
words, the function $\mu : \NN \to A$ is a ``curve'' which starts at $z$
and whose derivative at a given point $n$ is $s(\mu(n))$ -- this boils down to
stating that the curve $\mu$ is an integral curve for the vector field $s$
satisfying the initial condition $\mu(\A{Z}) = z$! Hence it may be possible to
understand iteration as a discrete counterpart of the Picard-Lindel\"of theorem,
which states that such integral curves always exist (locally).

It would be of great interest to extend $\lambda_\epsilon$ with an interation
operator and give its semantics in terms of differential (or difference)
equations. Studying recurrence equations using the language of differential
equations is a very useful tool in discrete analysis; for example, one can treat
the recursive definition of the Fibonacci sequence as a discrete ODE and use
differential equation methods to find a closed-form solution. We believe that in
a language which frames iteration in such terms may be amenable to optimisation
by similar analytic methods.

\bibliographystyle{spmpsci}
\begin{thebibliography}{10}
\providecommand{\url}[1]{{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{DOI~\discretionary{}{}{}#1}\else
\providecommand{\doi}{DOI~\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi

[1]
Alvarez-Picallo, M.: Change actions: from incremental computation to discrete
derivatives (2020)

[2]
Alvarez-Picallo, M., Lemay, J.S.P.: Cartesian difference categories.
\newblock In: International Conference on Foundations of Software Science and
Computation Structures, p. to appear. Springer (2020)

[3]
Alvarez-Picallo, M., Lemay, J.S.P.: Cartesian difference categories: Extended
report (2020)

[4]
Alvarez-Picallo, M., Ong, C.H.L.: Change actions: models of generalised
\newblock In: International Conference on Foundations of Software Science and
Computation Structures, pp. 45--61. Springer (2019)

[5]
Alvarez-Picallo, M., Peyton-Jones, M., Eyers-Taylor, A., Ong, C.H.L.: Fixing
incremental computation.
\newblock In: European Symposium on Programming. Springer (2019).
\newblock In press

[6]
Blute, R.F., Cockett, J.R.B., Seely, R.A.G.: Cartesian differential categories.
\newblock Theory and Applications of Categories \textbf{22}(23), 622--672

[7]
Bucciarelli, A., Ehrhard, T., Manzonetto, G.: Categorical models for simply
typed resource calculi.
\newblock Electronic Notes in Theoretical Computer Science \textbf{265},
213--230 (2010)

[8]
Cockett, J.R.B., Lemay, J.S.: Integral categories and calculus categories.
\newblock Mathematical Structures in Computer Science \textbf{29}(2), 243--308

[9]
Ehrhard, T.: An introduction to differential linear logic: proof-nets, models
and antiderivatives.
\newblock Mathematical Structures in Computer Science \textbf{28}(7), 995--1060

[10]
Ehrhard, T., Regnier, L.: The differential lambda-calculus.
\newblock Theoretical Computer Science \textbf{309}(1-3), 1--41 (2003).
\newblock \doi{10.1016/S0304-3975(03)00392-X}

[11]
Ehrhard, T., Regnier, L.: {Uniformity and the Taylor expansion of ordinary
\newblock Theoretical Computer Science \textbf{403}(2), 347--372 (2008)

[12]
Lavendhomme, R.: Basic concepts of synthetic differential geometry, vol.~13.
\newblock Springer Science \& Business Media (2013)

[13]
Lemay, J.S.P.: Exponential functions in cartesian differential categories.
\newblock arXiv preprint arXiv:1911.04790  (2019)

[14]
Manzonetto, G.: {What is a Categorical Model of the Differential and the
Resource $\lambda$-Calculi?}
\newblock Mathematical Structures in Computer Science \textbf{22}(03), 451--520

[15]
Vaux, L.: {$\lambda$-calculus in an algebraic setting}.
\newblock unpublished note  (2006)

[16]
Vaux, L.: {The algebraic lambda calculus}.
\newblock Mathematical Structures in Computer Science \textbf{19}(05),
1029--1059 (2009)

\end{thebibliography}

\end{document}