The High Time Resolution Universe Pulsar Survey X: Discovery of four millisecond pulsars and updated timing solutions of a further 12

Introduction

The bulk of the known pulsar population falls into two distinct groups when plotted on a period-period derivative diagram ($P$-$\dot{P}$ diagram). The normal or slow pulsars typically have spin periods between 0.1 and a few seconds, and have derived surface magnetic field strengths of $10^{11}$ to $10^{13}\rm\,G$. The second group has much lower magnetic field strengths of $10^{8}$ to $10^{9}\rm\,G$ and rapid spin periods measured in milliseconds. Members of this latter population are often referred to as the millisecond pulsars or MSPs. It is believed that MSPs are formed in binary systems in which the pulsar accretes matter from a companion star, gaining mass and angular momentum during the accretion process (Alpar et al. 1982; Tauris & van den Heuvel 2006, e.g.,). The pulsar is thus recycled and spun up to a very short spin period. At the same time the strength of its magnetic field is reduced, resulting in the typically small observed period derivative (Bhattacharya 2002, e.g.,). This formation scenario holds for most of the Galactic-field MSPs, whereas MSPs found in globular Clusters (GC) have more complicated histories, due to the significant probability of multiple exchange interactions with other cluster stars. In this paper we focus only on MSPs in the Galactic field.

MSPs are of particular interest mainly because of their typically high rotational stability, which combined with their short spin periods and narrow profile features, enable them to be timed precisely. This is in contrast to the younger group of normal or slow pulsars, which often show timing noise manifested as quasi-random variations in the rotational behaviour. MSPs are thus reliable and precise timing tools for a variety of astrophysical applications. For instance, MSPs have been employed in tests of gravity (Stairs 2003; Freire et al. 2012, e.g.,), for the detection of low frequency gravitational waves in Pulsar Timing Arrays (Yardley et al. 2011; van Haasteren et al. 2011, PTAs;), to provide measurements of neutron star masses to constrain the Equation of State (Demorest et al. 2010; Antoniadis et al. 2013), as precise clocks (Hobbs et al. 2012), in an array to constrain the Solar System ephemeris (Champion et al. 2010), and as an aid for the folding of gamma-ray photons to study the high-energy emission mechanism of pulsars (Abdo et al. 2009; Espinoza et al. 2013, e.g.,). At the same time, unique systems are constantly being discovered, including triple systems (Lynch et al. 2013), a highly-eccentric system (Champion et al. 2008) and the MSP J1719$-$1438 with an ultra-low mass companion (Bailes et al. 2011), challenging our theories of MSP formation and binary evolution.

To discover more MSPs and to improve our understanding of the MSP population as a whole, we began the High Time Resolution Universe (HTRU) survey in 2008. The HTRU is a blind pulsar survey of the Southern sky with the 64-m Parkes telescope (Keith et al. 2010) complemented by a twin survey in the north with the 100-m Effelsberg radio telescope (Barr et al. 2013, HTRU-North;). The surveys have benefited from recent advancements in technology and provide unprecedented time and frequency resolution, making them more sensitive to MSPs than previous efforts at these two telescopes. To date, the HTRU survey at Parkes has discovered more than 150 pulsars, of which 27 are MSPs.

The discovery of pulsars is however just a first step and, in fact, interesting science can usually only be revealed when a follow-up timing campaign is carried out. For MSPs, a coherent timing solution (i.e. when the number of rotations between every observation is well-determined) can be achieved typically within a few weeks of intense timing observations, providing preliminary determination of the rotational and orbital parameters, if any, of the MSP. This is the case for four newly-discovered MSPs presented in this paper. On the other hand, a timing campaign with a longer time baseline is necessary for improving the uncertainties of the timing solution and uncovering subtle details of each MSP system, such as proper motion, parallax, and possibly post-Keplerian binary parameters. This is demonstrated here by the further timing of 12 HTRU MSPs, the discoveries of which were first published two years ago (Bates et al. 2011; Keith et al. 2012; Bailes et al. 2011).

In this paper we describe the observations and analysis procedures used for obtaining the timing solutions (Section 2). We report on the discovery of four MSPs and present their initial timing solutions in Section 3, which include a discussion on the nature of the binary companions. We present the updated timing parameters for 12 further MSPs in Section 4, followed by a detailed discussion on various scientific implications arising from our measurements. Finally in Section 5 we present our conclusions.

OBSERVATIONS AND ANALYSIS

All 16 MSPs presented in this paper were discovered in 540-s-long integrations as part of the medium-latitude section ($-120°<\textit{l}<30°$, $|b|<15°$) of the HTRU survey. Survey observations were made using the 13-beam Multibeam receiver (Staveley-Smith et al. 1996) on the 64-m Parkes radio telescope. Full details of the survey parameters are given in Keith et al. 2010.

Follow-up timing observations were made at Parkes initially with a setup similar to that of the survey, employing the central beam of the same 13-beam Multibeam receiver at a centre frequency near 1.4 GHz and the Berkeley-Parkes-Swinburne Recorder (BPSR) with 1024 frequency channels incoherently dedispersed at a time resolution of 64$\mu$s. Later when the pulsar parameters were identified with sufficient accuracy, observations were carried out using the Digital Filterbank systems (DFBs) which are based on the implementation of a polyphase filter in FPGA processors with incoherent dedispersion. Coherently dedispersed data are collected by the CASPER Parkes Swinburne Recorder http://astronomy.swin.edu.au/pulsar/?topic=caspsr (CASPSR) and the ATNF Parkes Swinburne Recorder http://astronomy.swin.edu.au/pulsar/?topic=apsr (APSR). Pulsars with declination above $-35°$ are also being timed at the Jodrell Bank Observatory with the Lovell 76-m telescope, using a DFB backend and a ROACH backend. The latter is based on the ROACH FPGA processing board https://casper.berkeley.edu/wiki/ROACH and coherently dedisperses the data. Refer to Table 1 for the specifications of all observing systems employed.

Observations have also been taken at different frequencies at Parkes using the 10/50 cm receiver (Granet et al. 2005), to allow for precise dispersion measure (DM) measurements and to study any variations of pulsar profiles across frequencies. The various combinations of receivers and backends had central frequencies as listed in column 3 of Table 1. Note that predetermined offsets were applied to the observational data from Parkes to account for instrumental delay across observations with different backends in accordance with Manchester et al. 2013.

Timing observations of these 16 pulsars have first been made with an intense timing campaign within roughly their first year of discovery, and gradually decreased to weekly observation for the case of Jodrell Bank observations, whereas Parkes observations are more irregular with gaps ranging from days to months depending on telescope availability. Integration times vary from a few minutes to more than 2 hours, with longer observations for weaker pulsars to achieve adequate signal-to-noise (S/N) of at least 10.

We have used the psrchive data analysis package (Hotan et al. 2004) for data reduction. Each observation is corrected for dispersion and folded at the predicted topocentric pulse period, before finally summing over both frequency and time to produce an integrated profile. We align these profiles from each observation using an ephemeris created from the initial timing solution. This forms the basis of a noise-free analytic reference template, and we convolve the template with each individual profile to produce a Time of Arrival (TOA) (Taylor 1992). The DE421 Solar System ephemeris of the Jet Propulsion Laboratory (Folkner et al. 2009) was used to transform the TOAs to the Solar System barycentre. The tempo2 software package (Hobbs et al. 2006) was then used to fit a timing model to all TOAs, taking into account the astrometry, spin, and orbital motion of the pulsar. This process of cross-correlating a template with individual profiles can then be iterated to improve the quality of the model fit. We generate multiple TOAs per observation when possible, especially for the pulsars with small orbital periods. This is to make sure each TOA does not cover more than one tenth of an orbit, to avoid masking orbital information within a seemingly high S/N TOA. If simultaneous observations with different backends were taken, we include only one of the observations to avoid otherwise over-weighting duplicated TOAs.

All 16 MSPs in this work are in binary systems. The Damour-Deruelle (DD) timing model (Damour & Deruelle 1986) in tempo2 is a theory-independent description for eccentric binary orbits. However, for binaries with small eccentricities the location of periastron is not well-defined and using the DD timing model results in a high covariance between the longitude of periastron $(\omega)$ and the epoch of periastron $\rm(T_{0})$. A useful quantity to help choosing the best timing model is $xe^{2}$, where $e$ is the eccentricity and $x$ is the projected semi-major axis of the pulsar orbit as defined by:

with $a_{\rm{p}}$ being the semi-major axis, $i$ the orbital inclination and $c$ the speed of light. For pulsars with $xe^{2}$ smaller than the timing precision as represented by the RMS, we use the ELL1 timing model (Lange et al. 2001) alternatively. The ELL1 timing model avoids the covariance by using the Laplace-Lagrange parameters ($\epsilon_{1}=e\sin\omega$ and $\epsilon_{2}=e\cos\omega$) and the time of ascending node passage $(\rm{T_{asc}})$ instead of $\rm{T_{0}}$ as in the DD timing model.

Towards the end of the timing analysis procedure when the respective reduced $\rm\chi^{2}$ comes close to one, we can then assume a reliable fit is achieved which is only influenced by the presence of radiometer noise in the template. As a last step, we compensate for these systematic effects by calculating dataset-specific calibration coefficients (also known as ‘EFAC’ in tempo2). These coefficients are applied to scale the TOA uncertainties such that each final respective reduced $\rm\chi^{2}$ is unity.

In addition, full flux density and polarisation calibration are implemented for the four newly-discovered MSPs, in order to study their polarisation profiles. This analysis is not repeated for the rest of the 12 MSPs in this paper since their polarisation properties are already presented in Keith et al. 2012. With the only exception of PSR J1017$-$7156, a high-precision timing pulsar which is noticeably polarised in both linear and circular sense, we have fully calibrated the data to correctly assess the uncertainties on the TOAs.

To carry out the calibration we make use of Parkes DFB observations which record the four Stokes parameters in each frequency channel. We calibrate each observation for the differential gain and phase between the feed with an observation of the noise diode coupled to the receptors in the feeds. This calibration observation triggers a square-wave signal which is used to retrieve the true Stokes parameters, and it is important that this calibration is taken adjacent to the targeted pulsar observations. In addition, we correct for the non-orthogonality of the receptors in the Multibeam receiver by computing a model of the Jones matrix for the receiver using an averaged observation of the bright pulsar J0437$-$4715, in accordance with the ‘measurement equation modelling’ technique described in van Straten 2004, and we calibrate the flux density by using an averaged observation of Hydra A.

DISCOVERY OF FOUR MILLISECOND PULSARS

We present the discoveries of four MSPs in the HTRU survey, namely PSRs J1056$-$7117, J1525$-$5545, J1528$-$3828 and J1755$-$3716. They all have observations spanning more than one year, and their coherent timing solutions are shown in Table 2. All four are in binary systems.

PSR J1528$-$3828 and PSR J1056$-$7117 are likely to be formed from wide-orbit low mass X-ray binaries (LMXBs), leading to the formation of classic MSPs with Helium white dwarf (He-WD) companions. According to Tauris 2011, wide-orbit LMXBs with $P_{\rm{orb}}\geq 1$ day lead to He-WDs with masses between about $0.15$ to $0.46\,M_{☉}$.

PSR J1755$-$3716 has a relatively high median companion mass of $0.35\,M_{☉}$. Although this would fit in the above classification, the fact that PSR J1755$-$3716 has a spin period of 12.8 ms implies that the system is only mildly recycled. This, combined with its $P_{\rm{orb}}$ of just 11.5 days (which is too short for LMXB evolution to produce a $0.35\,M_{☉}$ WD, Tauris & Savonije 1999), indicates that its evolutionary track is more likely to have started from an intermediate mass X-ray binary pulsar (IMXB) accreting via early Case B Roche-lobe overflow (RLO) (Tauris 2011). The companion of PSR J1755$-$3716 is probably a CO-WD.

PSR J1525$-$5545 has a solar mass companion with a median mass of 0.99 $M_{☉}$ and an $P_{\rm{orb}}$ of 0.99 days. These fit the typical characteristics of binary evolution from a wide-orbit IMXB via Case C RLO and common envelope evolution (Tauris 2011). The companion is likely to be a massive CO-WD, or an ONeMg-WD if the orbital inclination is low.

2 Polarisation Profiles

Figure 1 shows the integrated polarisation profiles of the four MSPs in total intensity, linear and circular polarisation. We measure the Faraday rotation observed towards each pulsar by fitting the position angle (P.A.) variations across the 256 MHz band centred at 1369 MHz, and the plots shown here have their rotation measure (RM) corrected with the respective RMs as listed in Table 2. Multi-frequency data are included only if the S/N ratio is high enough, and are plotted here with an arbitrary alignment. None of the four MSPs are detectable at 3100 MHz with at least 1 hour of observation, except a tentative detection of PSR J1755$-$3716. At 732 MHz only PSRs J1056$-$7117 and J1525$-$5545 are detectable, both with limited S/N. Although pulsars typically have steep spectral indexes and thus higher flux at lower observing frequencies, our receiver system at 50 cm has a reduced sensitivity due to its higher system temperature and narrow bandwidth (Table 1). Hence we cannot comment if there is any profile evolution across frequency.

PSR J1056$-$7117 has a profile comprising three components. The emission of the middle component changes handedness in circular polarisation, whereas the S/N of the other two components are not sufficient for identifying the polarisation fraction. Linear polarisation is present in the middle component, although noisy. PSR J1525$-$5545 has a simple, single peak profile. It is almost completely unpolarised, and such low polarisation profile is typically associated with aligned gamma-ray and radio profiles (Espinoza et al. 2013). Although no Fermi Gamma-ray Space Telescope (Fermi) association has been reported for PSR J1525$-$5545 yet, it is worth following-up as the radio ephemeris improves with longer timing baseline. PSR J1528$-$3828 has a broad single peak profile with a hint of interpulse, and the P.A. is relatively flat over the profile. PSR J1755$-$3716 also has a profile formed of three components with some degree of linear polarisation in the middle component which is narrower compared to the total intensity, and the P.A. seems to show an ‘S-shaped’ swing.

UPDATED TIMING OF 12 HTRU MILLISECOND PULSARS

We have achieved considerable improvement in the timing accuracy for 12 HTRU MSPs compared with results published in their respective discovery papers (Bates et al. 2011; Keith et al. 2012; Bailes et al. 2011). This is thanks to the now longer timing baseline of more than three years in all cases, and only slightly less for PSR J1337$-$6423 which has 2.7 years of timing data. The timing parameters resulting from the best fits to the expanded set of TOAs are presented in Tables 3 and 4 for pulsars fitted with the ELL1 timing model and in Table 5 for pulsars fitted with the DD timing model.

In the following we discuss the physical implications arising from our timing measurements, including DM variations (Section 4.1), proper motion and transverse velocities (Section 4.2), intrinsic period derivatives (Section 4.3), binary companions and mass functions (Section 4.4), Galactic height distributions (Section 4.5), orbital eccentricities (Section 4.6), change in projected semi-major axis (Section 4.7), orbital period variation (Section 4.8), variation in the longitude of periastron (Section 4.9), and gamma-ray associations (Section 4.10).

Temporal variations in DM, due to turbulence in the ionised interstellar medium (ISM) and the changing line-of-sight to the pulsar, are in theory present in the TOAs of every pulsar (Petroff et al. 2013, see e.g.,). However this is typically not observable in slow pulsars since they have limited timing precision. In contrast, for MSPs such variations in DM can become significant and thus require special data treatment (You et al. 2007).

Indeed for the high-precision timing of PSR J1017$-$7156 we identified significant temporal variations in its DM measurement, implying changes in the electron density in the ISM along the line-of-sight over a time scale of a few months. We have attempted to model this variation via three correction methods, firstly by fitting DM variations across short ranges of TOAs while holding fixed all other parameters, secondly by including higher order DM derivatives and thirdly by the DM model described in Keith et al. 2013. In Figure 2 we plot the manually identified values of DM across every few TOAs in black. We plot the best-fit curve from the timing solution of tempo2, employing up to eight DM-derivatives as the green dashed line. We plot the DM model derived using the method outlined in Keith et al. 2013 as red crosses, and the red solid line joining them shows the resulting DM model. It can be seen that the DM derivatives provide a smooth fit to the DM variations, however there are still small scale variations that are not properly accounted for. On the other hand, the DM model essentially creates a linear interpolation between DM offsets identified at specific epochs (note that here we have adopted a gap of 50 days between successive DM offsets), and hence can be tailor-made to follow more closely variations on all scales. We conclude that the DM model of Keith et al. 2013 gives a more successful fit and hence have adopted this for the timing solution of PSR J1017$-$7156.

2 Proper motion and transverse velocities

The proper motion $(\mu)$ of a pulsar introduces a positional offset over time and is measurable from pulsar timing data. Within our sample of 12 MSPs with extended timing solutions, we have measured five new proper motions with significances greater than 3-$\sigma$, for PSRs J1017$-$7156, J1125$-$5825, J1446$-$4701, J1708$-$3506 and J1719$-$1438. PSR J1811$-$2405 is very close to the ecliptic plane with $(\lambda,\beta)=(272.586°,-0.675°)$ which means its proper motion in ecliptic latitude ($\mu_{\beta}$) cannot be well-constrained. With a $\lambda$ so close to 270$°$, the translation from ecliptic frame to equatorial frame would have almost no rotation. This implies that the large uncertainty associated with $\beta$ is only inherited in the declination, $\delta$, without also contaminating the right ascension, $\alpha$. Hence for PSR J1811$-$2405 we can choose to continue using the equatorial coordinates and we fixed $\mu_{\delta}$ at zero for the rest of the analysis. For the four newly-discovered MSPs in this paper, their time spans are not yet long enough for proper motion to be detected with significance.

From $\mu$ and their respective pulsar distances, $d$, we can derive their corresponding transverse velocities, $V_{\rm{T}}$, with the following equation,

In this work we have calculated pulsar distances based on the NE2001 electron density model (Cordes & Lazio 2002) and we assume an associated uncertainty of 25$\%$ for each DM-derived distance. MSP proper motion measurements are relatively rare and hence there are not many derived velocities (only about 40 currently published values in the literature), making it difficult to place constraints on MSP velocity distribution models. The latest MSP velocity discussions can be found in Toscano et al. 1999 and Hobbs et al. 2005, proposing an average velocity for recycled MSPs of $85\pm 13\rm\,km\,s^{-1}$ and $87\pm 13\rm\,km\,s^{-1}$ respectively. Hobbs et al. 2005 also quoted a median velocity for recycled MSPs of $73\rm\,km\,s^{-1}$. Our new $V_{\rm{T}}$ measurements largely agree with these previous results (refer to Table 3 to 5). Note that we believe the high VT of 670 and 350$\rm\,km\,s^{-1}$ reported for PSRs J1708$-$3506 and J1731$-$1847 in Bates et al. 2011 should in fact be corrected to more modest values of $70\pm 20$ and $80\pm 40\rm\,km\,s^{-1}$ respectively.

3 Observed and inferred intrinsic period derivatives

The vast majority of pulsars are rotation-powered objects and hence their respective period derivatives $(\dot{P})$ are fundamental to their identities. The observed period derivatives $(\dot{P}_{\rm{obs}})$ however contain a contribution from kinematic effects (Shklovskii 1970) and acceleration due to the Galactic potential (Damour & Taylor 1991). Determination of the intrinsic period derivative is important for properly placing pulsars in the $P$-$\dot{P}$ diagram from which physical conclusions (such as magnetic field strength, characteristic ages) may be drawn. To obtain the intrinsic period derivative $(\dot{P}_{\rm{int}})$ we employed the following equation,

The term $\dot{P}_{\rm{shk}}$ accounts for the apparent acceleration that arises from the transverse motion of the pulsar. It is related to the pulsar spin period, $P$, the proper motion, $\mu$, and the pulsar distance, $d$, by the following equation from Shklovskii 1970,

The term $\dot{P}_{\rm{gal}}$ accounts for difference in the line-of-sight components of the acceleration of the pulsar and the Solar System under the influence of the Galactic gravitational potential. There exist several Galactic potential models in the literature, and we have chosen the one described in Paczynski 1990. This model reproduces a flat rotation curve and uses a Solar Galactocentric distance $R_{0}$ of 8 kpc and a Solar Galactic rotation velocity of $220\rm\,km\,s^{-1}$.

Table 6 lists the $\dot{P}$ contributions as calculated for the 12 MSPs with updated timing solutions in our sample. Monte Carlo simulations with 1,000,000 runs per pulsar have been used to estimate the associated error. Note that the errors in $\dot{P}_{\rm{shk}}$ and $\dot{P}_{\rm{gal}}$ do not reflect the effect of errors in the distance estimates. The results are illustrated in Figure 3, which is a $P$-$\dot{P}$ diagram around the region where MSPs are located. The $\dot{P}_{\rm{obs}}$ and the corrected $\dot{P}_{\rm{int}}$ of the 12 MSPs studied in this paper are plotted, together with other known pulsars in this region.

Some of the results (noticeably those of PSR J1337$-$6423 and J1502$-$6752) have large associated errors and should be considered with caution. One reason is that $\dot{P}_{\rm{shk}}$ relies on the square of $V_{\rm{T}}$, which is in turn dependent on proper motion as seen from Equation (2). Hence $\dot{P}_{\rm{shk}}$ is only meaningful for MSPs with well-constrained proper motion measurements. Additionally, $\dot{P}_{\rm{gal}}$ is dependent on the distance of the pulsar, $d$. As mentioned in Section 4.2, the DM-derived distance is thought to have $\sim 25\%$ error, and can be much larger for individual pulsars.

Disregarding these two unconstrained measurements, PSR J1017$-$7156 stands out with one of the smallest inferred intrinsic $\dot{P}$ at a value of $1.2\times 10^{-21}$. We are aware that if red noise is present in the data this could potentially also contaminate our $\dot{P}$ measurement. However if we include the frequency second derivative in the model fit in an attempt to whiten the data with a quadratic component, the $\dot{P}$ measurement remains statistically consistent. PSR J1017$-$7156 is thus located at the bottom left of the $P$-$\dot{P}$ diagram, which yields a characteristic age, $\tau_{\rm{c}}\equiv$ $P/(2\dot{P})$, of 31 Gyr, i.e. larger than the Hubble age. Note that $\tau_{\rm{c}}$ is by no means a reliable age indicator for MSPs, since it is only applicable for pulsars which have a braking index $n=3$ and an initial spin period ($P_{0}$) much less than the current spin period, which is not thought to be the case for MSPs. However for MSPs with such small $\dot{P}_{\rm{obs}}$ like that of PSR J1017$-$7156, we can deduce that the MSP was probably born with small initial period derivative and must not have moved very far from its current location on the $P$-$\dot{P}$ diagram since its birth (Tauris et al. 2012). The derived surface magnetic field strength at the equator Note throughout the paper we differentiate between the derived surface magnetic field at the polar region ($B_{\rm{p}}$), and that of the equatorial region ($B_{\rm{eq}}$) which is only half the strength comparing to the polar region. of PSR J1017$-$7156 is also at one of the lowest known at $5.3\times 10^{7}$ G.

3.2 PSR J1801−-3210

There is one peculiar case, PSR J1801$-$3210, for which no significant period derivative has been measured, even with more than 4 years of timing data. The best-fit solution in tempo2 shows a $\dot{P}_{\rm{obs}}$ of $-4\pm 4\times 10^{-23}$, an extremely small number compared to that of typical MSPs ($\dot{P}_{\rm{obs}}$ of the order of $10^{-19}$ to $10^{-20}$). A 2.7-$\sigma$ $\dot{P}$ value of $0.265(97)\times 10^{-20}$ was presented in the initial discovery paper by Bates et al. 2011 which at that time had just over one year of timing data, however this value is inconsistent with our current longer time baseline TOAs. Shortening our data span to the same epoch as that in Bates et al. 2011 results in an unconstrained $\dot{P}_{\rm{obs}}$ measurement of 0.2(20)$\times 10^{-20}$, eliminating the possibility of an actual change in period derivative over time.

Referring again to Equation (3), proper-motion-induced $\dot{P}_{\rm{shk}}$ has an always positive contribution to $\dot{P}_{\rm{obs}}$, so that $\dot{P}_{\rm{int}}$ will be even smaller. $\dot{P}_{\rm{gal}}$ however could have a positive or negative contribution depending on the relative location of the pulsar in the Galaxy with respect to the Earth.

PSR J1801$-$3210 has a proper motion measurement of $15(7)\rm\,mas\rm\,yr^{-1}$, corresponding to a positive $\dot{P}_{\rm{shk}}$ of the order of $10^{-20}$. The Paczynski 1990 Galactic potential model shows that at the NE2001 DM-derived distance of 4 kpc PSR J1801$-$3210 would be accelerated away from the Sun, giving a positive $\dot{P}_{\rm{gal}}$ of the order of $10^{-21}$ (Table 6) to further decrease the already negative $\dot{P}_{\rm{obs}}$. Even if we assume the proper motion to be zero to get the smallest possible contribution from $\dot{P}_{\rm{shk}}$, we still cannot overcome this apparent negative $\dot{P}_{\rm{int}}$ at the given DM distance of 4 kpc, since the $\dot{P}_{\rm{gal}}$ is positive and dominates the tiny $\dot{P}_{\rm{obs}}$ of $10^{-23}$. We acknowledge that the Paczynski 1990 model consists basically of only three elements: a bulge, a disk and the surrounding halo. However this is considered a valid approximation, and for example the effect of spiral arm structure should not significantly skew the model.

In the following we consider other potential explanations to this apparent negative $\dot{P}_{\rm{int}}$, i.e. effects that would have contributed to the $\dot{P}_{\rm{obs}}$ but are not yet accounted for in Equation (3). We discuss the cases of (a) acceleration due to local stars; (b) acceleration due to giant molecular clouds (GMCs); and (c) acceleration due to a third orbiting object if PSR J1801$-$3210 is in a triple system.

If there exists a third body (with mass $M_{3}$) located near the pulsar, in a direction towards the Earth and close to the line-of-sight, it will potentially cause a radial acceleration of PSR J1801$-$3210 towards the Earth. We can express the mass of the third body required to produce a $\dot{P}$ contribution of $\dot{P}_{\rm{M_{3}}}$ as,

where $r$ is the distance between the third body and the pulsar, $G$ is Newton’s gravitational constant and $\theta$ is the angle between the direction from the pulsar to the third body and the direction from the pulsar to the Sun. We imagine the scenario of $\theta\approx 0°$ where the line-of-sight acceleration induced on the pulsar is the largest, and we first examine the potential contribution from stars located near the pulsar. The probability distribution of fluctuation in Galactic acceleration due to local clustering centres has been studied in the literature (Holtsmark 1919, see e.g.,), and based on Equation (3.1) and (3.5b) in Damour & Taylor 1991 one finds for PSR J1801$-$3210 at 1-$\sigma$ confidence level,

where $\dot{P}_{*}$ is the potential period derivative contribution from nearby stars, $\hat{M}$ is the average of mass taken over the mass spectrum of the attracting centres and we use the same value of 1 $M_{☉}$ as in Damour & Taylor 1991. The local stellar-mass density, $\rho_{☉}$, has a value of 0.06$\,M_{☉}\,\rm{pc}^{-3}$ according to Mihalas & Binney 1981 and the stellar-mass density, $\rho$, can be extrapolated by,

where $R_{0}$ is the aforementioned Solar Galactocentric distance at 8 kpc. $L_{\rm{disk}}$ is the stellar disk scale length and $z_{\rm{h}}$ is the scale height of the stellar disk component, which from the most recent literature by Bovy & Rix 2013 $L_{\rm{disk}}=2.15\rm\,{kpc}$ and $z_{\rm{h}}=0.4\rm\,{kpc}$. $R$ is the distance of the pulsar from the Galactic Centre, and for the case of PSR J1801$-$3210 it is approximately 4 kpc as derived from the NE2001 model. This corresponds to a Galactic height, $z$, of 0.32 kpc. Substituting these into Equations (7) and (6) gives $\rho=0.17\,M_{☉}\,\rm{pc}^{-3}$ and a tiny $\dot{P}_{*}$ of the order of $10^{-24}$ which is unlikely to have led to the negative $\dot{P}_{\rm{int}}$. To appreciate the improbability of this scenario we can also hypothesise a nominal $\dot{P}_{*}$ of the order of $-10^{-21}$. From Equation (6) this would require $\rho$ to be more than 300$\,M_{☉}\,\rm{pc}^{-3}$ and nowhere along the line-of-sight direction of PSR J1801$-$3210 has such high stellar-mass density.

Alternatively let us consider the contribution from GMCs, and again we assume that there exists such an acceleration acting upon the pulsar towards the Earth which induces a nominal $\dot{P}_{\rm{GMC}}$ of the order of $-10^{-21}$. GMCs typically have masses between $10^{3}$ to $10^{7}\,M_{☉}$; substituting this into Equation (5) corresponds to a distance, $r$, of about 2 to 190 pc from the pulsar. No GMCs are known to exist near PSR J1801$-$3210, but not all GMCs have necessarily been detected, so this possibility cannot be ruled out. It may also be that multiple smaller molecular clouds (also known as Bok globules) act together to accelerate PSR J1801$-$3210 in our direction.

Another possible candidate of this third body could be a tertiary star or an exoplanet orbiting PSR J1801$-$3210 in a weakly-bounded hierarchical triple orbit. This third component would accelerate the pulsar system towards it, and hence if the third component happened to provide a net acceleration on PSR J1801$-$3210 towards the Earth it would lead to the negative $\dot{P}_{\rm{int}}$ like in the case of a GMC as mentioned above. We can achieve the same $\dot{P}_{\rm{exo}}$ of the order of $-10^{-21}$, for example with an Earth-sized exoplanet at distance of $\sim$20 AU in an orbit of $\sim$70 years around PSR J1801$-$3210, or a Jupiter-sized exoplanet at a distance of $\sim$400 AU in a large orbit of $\sim$6000 years, assuming circular orbit.

The relative motion between the pulsar system and the exoplanet would have induced variations in the acceleration, as well as variations in the second derivative of spin frequency, $\ddot{\nu}$ (Backer et al. 1993). We do not have a significant measurement of $\ddot{\nu}$ except a 2-$\sigma$ upper limit of $8\times 10^{-26}\rm\,s^{-3}$. This thus excludes the existence of a nearby exoplanet and favours the case of a further-out heavier object. However at the same time, for a third orbiting object to stay bound with the pulsar system, a very strict limit on the post-supernova (SN) recoil velocity of the inner binary is required (Hills 1983). Precisely, the recoil velocity has to be no more than $30\rm\,km\rm\,s^{-1}$ and 7$\rm\,km\rm\,s^{-1}$ for the case of an Earth-mass and a Jupiter-mass exoplanet respectively. According to simulations by Tauris & Bailes 1996, the recoil velocity of any surviving binary is expected to be larger than 20$\rm\,km\rm\,s^{-1}$, even for a symmetric SN explosion, unless the pulsar formed via an accretion-induced collapse of a white dwarf (Nomoto et al. 1979). Hence, we are inclined to exclude a very distant third body with a Jupiter mass, and notice that a closer Earth-mass object would require quite some fine-tuning in the SN event to remain bound. To summarise, we conclude that this scenario of an exoplanet is possible but unlikely.

Finally, we consider the possibility that the NE2001 DM-derived distance of 4 kpc is significantly wrong, hence locating PSR J1801$-$3210 in a different quadrant of the Galaxy which would reverse the direction of the Galactic potential and the sign of $\dot{P}_{\rm{gal}}$. In Figure 4 we plot the various $\dot{P}$ contribution as a function of distance along the line-of-sight of PSR J1801$-$3210. It can be seen that in the limiting case of $\dot{P}_{\rm{shk}}$ being zero, we can achieve a positive period derivative beyond a distance of 8 kpc, and can reach an upper limit of $\dot{P}_{\rm{int}}$ of $3\times 10^{-20}$ at a distance of 8.5 kpc. At a distance of 8 kpc, the NE2001 model requires a corresponding DM of 326.1$\rm\,cm^{-3}\rm\,pc$ which is inconsistent with the well-constrained DM measurement of PSR J1801$-$3210 of only 177.713(4)$\rm\,cm^{-3}\rm\,pc$. However other electron density models give very different results. For example the TC93 model (Taylor & Cordes 1993) requires a corresponding DM of only 227.0$\rm\,cm^{-3}\rm\,pc$, whereas including a thick disk component to the TC93 model (Schnitzeler 2012) predicts an even smaller corresponding DM of 185.5$\rm\,cm^{-3}\rm\,pc$, which is only a factor of 1.07 from our measured value. These large discrepancies between various models reflect uncertainties in the electron density distribution along this line-of-sight, and thus it seems plausible that the DM-derived distances of PSR J1801$-$3210 have been underestimated. PSR J1801$-$3210 is located at $(l,b)=(358.922°,-4.577°)$, a distance of at least 8 kpc in this direction would put PSR J1801$-$3210 just beyond the Galactic Centre, hence reversing the direction of $\dot{P}_{\rm{gal}}$. In any case, we suggest that PSR J1801$-$3210 would serve as an important test pulsar for improving future electron density models.

Otherwise, if PSR J1801$-$3210 has indeed an extremely small $\dot{P}_{\rm{int}}$ it would imply an exceptionally small surface magnetic field. Popular theories on the pulsar emission mechanism require electron-positron pair production, and the longer the spin period of the pulsar, the larger the potential needed to power the particle acceleration (Beskin et al. 1988, see for example). The following implication is known as the ‘pulsar death line’, which predicts for a particular pulsar spin period, there exists a lower limit of period derivative and surface magnetic field for which radio emission can be produced. Therefore, we can derive a lower limit of $\dot{P}_{\rm{int}}$ for PSR J1801$-$3210 to stay above the pulsar death line. We adopt the theoretical study from Chen & Ruderman 1993 which described three possible death lines also plotted in Figure 3. If we take the lowest limiting case imposed by death line B, we derive a lower limit of $\dot{P}_{\rm{int}}=7.9\times 10^{-24}$ and a corresponding surface magnetic field at the equator ($B_{\rm{eq}}$) of $7.8\times 10^{6}$ G. We note that this derivation assumes a contribution only from a model of a vacuum magnetic dipole. However as discussed by Tauris et al. 2012, if the spin-down torque caused by the plasma current in the magnetosphere (Spitkovsky 2006) is also taken into account, the realistic surface magnetic field would even be lower, by at least a factor of $\sqrt{3}$.

4 Binary companions and mass functions

A plot of mass function versus orbital period is a standard way of distinguishing different types of binary systems and can be used to gain insight into the nature of the binary companion, as shown in Figure 5. Indeed it can be seen immediately that PSR J1719$-$1438 occupies an otherwise empty region in the bottom left corner of this figure, as a result of its uniquely light, planet-mass companion. This has been extensively discussed in the literature (Bailes et al. 2011; van Haaften et al. 2012, e.g.,) so will not be further elaborated in this paper.

A cluster of pulsars can be seen in the left side of Figure 5, with $P_{\rm{orb}}$ $\leq$ 1 day and mass functions between $10^{-7}$ to $10^{-4}\,M_{☉}$. They are considered descendants of close LMXB systems, resulting in the formation of a binary with an ultra-light companion (Tauris 2011), also known as the ‘very low-mass binary pulsars’ (VLMBPs). In our sample we have three MSPs that fit into this category, namely PSRs J1446$-$4701, J1502$-$6752 and J1731$-$1847.

Some of the VLMBPs exhibit eclipses and are typically referred to as black widow pulsars (Roberts 2013, BW;). Eclipses have already been reported for PSR J1731$-$1847 by Bates et al. 2011, but not for PSRs J1446$-$4701 nor J1502$-$6752. Freire 2005 proposed a correlation between the possibility of observing eclipses and orbital inclination for these VLMBPs in GC. The essence of the idea is that the companions of these VLMBPs have a narrow intrinsic mass distribution, and subsequently whether a VLMBP shows eclipses or not, becomes exclusively dependent on its orbital inclination. In other words, a VLMBP viewed relatively face-on (low inclination) is less likely to be observed as an eclipsing system and will also have a smaller mass function, and vice versa.

While this hypothesis seems to work well in GCs, there has not yet been a similar study on the non-GC associated VLMBP population. We have compiled all related literature, and colour-coded in Figure 5 the known ‘eclipsers’ as red circles and the ‘non-eclipsers’ as blue diamonds. Two distinct groups composed of ‘eclipsers’ and ‘non-eclipsers’ do seem to exist, with only one outlier, PSR J1311$-$3430, which is the tightest binary pulsar known with a $P_{\rm{orb}}$ of just 93 minutes (Romani et al. 2012; Pletsch et al. 2012; Ray et al. 2013). But this pulsar may have evolved from an ultra-compact X-ray binary (UCXB), hence belonging to a different population (van Haaften et al. 2012) and might not be applicable to the hypothesis as mentioned above. Disregarding this system, it is striking to see a bimodal distribution. Particularly interesting is that there is no non-eclipsing system found within the red cluster of ‘eclipsers’, although from a pulsar searching point-of-view these kinds of systems should in fact be easier to detect due to their non-eclipsing nature.

Plotted as a dotted line in the zoomed-in panel of Figure 5 is our nominal split between the ‘eclipsers’ and the ‘non-eclipsers’, representing a dividing mass function of $6.7\times 10^{-6}\,M_{☉}$. We assume a pulsar mass of $1.7\,M_{☉}$, and an orbital inclination of 70$°$ to postulate a lower limit on inclination which eclipses can be observed. This dividing mass function would then correspond to $m_{\rm{c}}=0.029\,M_{☉}$, which is also within the range of typical companion masses of BWs as shown in Chen et al. 2013. Indeed orbital eclipses are observed for PSR J1731$-$1847 which has a median $m_{\rm{c}}$ of 0.0385 $M_{☉}$ and lies above the dotted line, whereas no eclipse is observed for PSRs J1446$-$4701 and J1502$-$6752 with lower companion masses (median $m_{\rm{c}}$ of 0.022 $M_{☉}$ and 0.025 $M_{☉}$ respectively) located below this line. These measurements are in agreement with Freire 2005.

5 Galactic height distribution

Based on theoretical grounds we expect an anti-correlation between the absolute Galactic height and the inferred mass function of binary pulsars. The reason is the following: assuming that the momentum kick imparted to a newborn neutron star during the SN explosion is independent of exterior parameters, such as the mass of the companion star, the resulting systemic recoil velocity is larger for systems with smaller companion star masses (and thus smaller mass functions) as a simple consequence of conservation of momentum. Since the acquired amplitude of the Galactic motion of the system only depends on the systemic recoil velocity, we therefore expect the above mentioned anti-correlation between the distribution of observed Galactic heights and the measured mass functions of pulsar binaries. Some theoretical studies (Tauris & Bailes 1996, e.g.,) have suggested the possibility of a weak relation between orbital period and systemic recoil velocity of pulsar binaries. However, Gonzalez et al. 2011 found no observational evidence for such a relation based on the 2D velocities of binary MSPs. Thus we disregard orbital periods in the following discussion.

Our sample of 16 MSPs has a wide distribution of mass functions, from PSR J1719$-$1438 with an ultra-low mass companion and a mass function of 7.8$\times 10^{-10}\,M_{☉}$ to PSR J1525$-$5545 with a massive WD companion and a mass function of 0.11$\,M_{☉}$. With the addition of these systems, we investigate whether there exists a correlation between mass function and vertical distance from the Galactic plane $(|d\sin{b}|)$. We have taken our sample of MSPs from the ATNF Pulsar Catalogue http://www.atnf.csiro.au/people/pulsar/psrcat/ (Manchester et al. 2005) and an online MSP catalogue maintained by Lorimer http://astro.phys.wvu.edu/GalacticMSPs/GalacticMSPs.txt. We have included the 16 MSPs in this work and also six additional newly-discovered HTRU MSPs (Ng et al. in prep; Thornton et al. in prep). All recycled MSPs in binary systems are considered, provided that they are not associated with a GC or extragalactic, which amounts to 164 MSPs in total. We continue to use the Cordes & Lazio 2002 model of Galactic electron density to derive the distances of all known pulsars in order to calculate their respective Galactic heights. Independent distance measurements are available for 19 binary systems and we use these, instead of the DM distances, when calculating their Galactic heights. In Figure 6 we plot the absolute Galactic heights against mass functions, and we classify the nature of each of the binary companions in accordance with the description in Tauris et al. 2012. This results in five binary groups, namely those with ultra-light (UL) companions, with He-WD companions, with massive CO or ONeMg-WD companions, neutron-star$-$neutron-star (NS-NS) systems and those with main-sequence star (MS) companions. For the rest of the discussion we set aside the nine NS-NS systems, since they were born with two SN explosions (hence received two kicks) and would complicate our discussion.

Table 7 summarises the statistical distribution of Galactic height for each of the binary groups mentioned above, from which we draw two main interpretations. Firstly, the heavier systems tend to stay closer to the plane, as seen for example from the MS systems with a mean Galactic height of only 0.06 kpc, whereas the lightest UL systems tend to be found at a higher Galactic height with a mean of 0.52 kpc. Secondly there is a larger scatter in the height distribution of the lighter systems, whereas the heaviest MS system are found almost exclusively within the Galactic plane. We note that a potential caveat here is that the ages of the MSPs might also have an influence on the Galactic height scattering. For example the fully-recycled He-WD binaries are generally older and hence might have more time to scatter away from the Galactic plane, whereas the less-recycled binaries with heavier companions tend to be younger. In addition, there is also a longer time interval between the SN explosion and the formation of the MSP for systems with UL and with He-WD companions, because their low-mass progenitors have much longer nuclear evolution timescales. Nonetheless, this does not change the outcome of the overall picture in Figure 6, explicitly that the distribution of the total mass of binary systems is inversely-related to the Galactic height distribution.

We are aware that the MSP distribution depicted in Figure 6 is skewed by another observational bias. That is from a pulsar searching point-of-view, pulsars with shorter spin periods, meaning the more recycled UL and He-WD systems, are more difficult to be discovered at higher DM regions, for example deep in the Galactic plane. This is because short spin period pulsars are more vulnerable to dispersion smearing and interstellar scattering. However, the less recycled massive WD and MS systems have longer spin periods, and we should have a relatively more uniform ability to detect them whether they are in the Galactic plane where DM is high or out of the plane.

This leads to two further implications. The first is that the smaller Galactic heights of the heavier systems are genuine, since if massive WD or MS systems exist at high Galactic heights we would have been more likely to have discovered them, given that we have detected the in-theory more difficult He-WD at those Galactic heights. The second is that this gives an explanation to the lack of light systems at small Galactic heights in the Galactic plane, resulting in the sparsely populated region in Figure 6 below 0.2 kpc and for mass function less than $10^{-3}\,M_{☉}$. Indeed a large number of the UL systems at high Galactic heights are only discovered thanks to the Fermi Large Area Telescope (Atwood et al. 2009, LAT;), which has much less ability to detect pulsars in the Galactic plane due to confusion with background emission.

These results show that the observed MSP distribution is not as isotropic as previously thought (Johnston & Bailes 1991, see for example) prior to the latest generation of pulsar surveys with improved backends, which have allowed us to probe a much bigger volume within the Galaxy. Conventional MSP population synthesis using the scale factor method typically takes into consideration only the pulsar luminosities (Levin et al. 2013, see for example), and we suggest that including the mass function as an extra parameter could be a potential improvement for future population studies.

6 Orbital eccentricity

We have measured initial eccentricities for the four newly-discovered binary MSPs and improved precision for the eccentricities of the 12 previously published MSPs, except for PSRs J1502$-$6752 where only upper limits can be achieved. Figure 7 shows a plot of orbital period versus orbital eccentricity and the 16 MSPs in this work are marked together with 1- and 2-$\sigma$ uncertainties of their eccentricities. The dotted lines denote the eccentricity predicted by the convective fluctuation-dissipation theory of Phinney 1992, applicable to binary systems formed by stable mass transfer from a Roche-lobe filling red giant. It can be seen that our MSPs with He-WD companions (plotted as red stars in Figure 7) largely agree with the predictions of Phinney 1992. Within the 2-$\sigma$ eccentricity measurement uncertainties, only PSRs J1017$-$7156, J1811$-$2405 and J1801$-$3210 lie outside the 95% confidence-level range (the first one above and the latter two below). However as seen in Figure 7 they have the same scatter as the rest of the MSP population. In addition, these three pulsars have typical He-WD companions and their spin periods indicate highly recycled systems. Therefore, we find little evidence for unusual evolutionary scenarios for these three pulsars.

The low eccentricity of $e=2.1\pm 1.1\times 10^{-6}$ of PSR J1801$-$3210 combined with its large orbital period of $P_{\rm{orb}}=21$ days makes it a ‘wide-orbit binary millisecond pulsar’ (WBMSP), and an interesting object to be employed for tests of the strong equivalence principle (SEP) as described in Damour & Schäfer 1991; Stairs et al. 2005; Freire et al. 2012. The basic idea being that in the case of SEP violation, the extreme difference between the gravitational binding energy of the heavy neutron star and its much less compact companion star implies that they would experience different accelerations in the presence of an external gravitational field (Nordvedt effect). This translates to an observable effect, most prominent in systems with small eccentricity and wide orbits, that the eccentricity would oscillate between the minimum and maximum value. The dashed line overplotted on Figure 7 indicates $e\propto{P_{\rm{orb}}}^{2}$, a figure-of-merit for a SEP test. With a ${P_{\rm{orb}}}^{2}/e$ ratio of $2.1\times 10^{8}\rm\,day^{2}$, PSR J1801$-$3210 thus provides the best test for SEP together with PSRs J1835$+$1303 and J0407$+$1607 as detailed in Gonzalez et al. 2011. Note that although PSR J1711$-$4322 appears to lie close to the figure-of-merit in Figure 7, it is in fact not usable for this SEP test (Kehl & Krieger 2012).

7 Change in projected semi-major axis, x˙\dot{x}

For PSR J1017$-$7156 we determine a change in projected semi-major axis $(\dot{x})$ of $9.1\pm 1.7\times 10^{-15}$. The projected semi-major axis, $x$, is related to the semi-major axis, $a_{\rm{p}}$, and the inclination, $i$, by Equation (1). Hence a measurement of $\dot{x}$ could be due either to a physical change of the intrinsic orbit size as measured by $a_{\rm{p}}$, or to a change in $i$, or both.

In the case of an actual change in $a_{\rm{p}}$ due to gravitational wave emission, we would expect this to also be reflected in a detection of $\dot{P}_{\rm{orb}}$ (Peters 1964). From this we can predict the corresponding observable change in $a_{\rm{p}}\sin{i}/c$ to be of the order of $10^{-21}$ for PSR J1017$-$7156, which is many orders of magnitude too small to be observed. So we conclude that the observed $\dot{x}$ is most likely due to an apparent change in the orbital inclination as a result of proper motion affecting the viewing geometry. This effect has been first proposed by Arzoumanian et al. 1996 and Kopeikin 1996 using,

In this equation proper motion has a total magnitude of $\mu$ and a position angle of $\Theta$, whereas $\Omega$ is the position angle of the line of nodes.

To assess if any physical constraints of the orientation of the line of nodes in relation to the direction of the proper motion (i.e. $\Theta\,-\,\Omega$) can be subsequently drawn, one must compare the uncertainty of the measured $\dot{x}$ with the product of $\mu$ and $x$. For PSR J1017$-$7156 we have $\mu{x}=7.6\times 10^{-15}$, which is indeed in the same order of magnitude as compared to our $\dot{x}$ measurement, and can already provide constraints to the possible ranges of $\Omega$. Future improved timing precision and additional information, such as constraints on or detection of a Shapiro delay, will allow us to extract more information on the binary systems, including mass measurements. None of the other MSPs in this paper have a detectable $\dot{x}$ yet and are unlikely to be measurable in the near future. With the possible exceptions of PSRs J1125$-$5825 and J1708$-$3506, which both have $\mu$$x$ of the order of $10^{-14}$, we can quote a marginal $\dot{x}$ limit of $1.6\pm 2.0\times 10^{-14}$ and $-9\pm 6\times 10^{-14}$ respectively. Hence they might achieve reliable $\dot{x}$ measurements with additional timing data.

8 Orbital period variation, P˙orb\dot{P}_{\rm{orb}}

We measure an orbital period variation $(\dot{P}_{\rm{orb}})$ in PSR J1731$-$1847. However rather than due to gravitational-wave damping, the $\dot{P}_{\rm{orb}}$ observed in this case is more likely due to the eclipsing nature of PSR J1731$-$1847, a BW system, inducing orbital interaction. We refer to Lazaridis et al. 2011 for a detailed discussion of such orbital period variations caused by changes in the gravitational quadrupole moment of a tidally interacting BW system. For the case of PSR J1731$-$1847, a straight-forward fit of $\dot{P}_{\rm{orb}}$ is not adequate, since the orbital period exhibits quadratic changes over the last three years. We have achieved the best fit using the BTX model (Nice, D., unpublished) implemented in tempo2, taking into account the orbital frequency changes up to the second order term (i.e. ${n_{\rm{b}}}$, $\dot{n}_{\rm{b}}$, $\ddot{n}_{\rm{b}}$). The phase $(\phi)$ of the orbit is thus a function of the binomial expansion of the $n_{\rm{b}}^{(k)}$ terms, where $k$ denotes the ${k^{\rm{th}}}$ derivative with respect to time. At any particular time, $t$, the phase $\phi$ can be represented by,

To get a better visualisation of the change of the orbit over time, we express this phase shift as the shift of the epoch of periastron ($T_{0}$). One can consider that a positive phase shift corresponds to an earlier arrival of the observed periastron, $T_{\rm{0,obs}}$, as compared to the predicted arrival of the periastron, $T_{\rm{0,pre}}$. The result is a negative $\Delta T_{0}$, which also symbolises a decrease in $\dot{P}_{\rm{orb}}$,

Figure 8 shows this $\Delta T_{0}$ as derived from the $n_{\rm{b}}^{(k)}$ terms of the BTX model fit in tempo2. It can be seen that the orbit of PSR J1731$-$1847 shrinks until approximately MJD 55800 but gets wider after. We identified manually a value of $T_{0}$ across every few TOAs, while holding fixed all other parameters (shown by black points in Figure 8). The BTX model results in a close agreement. We remark, however, that this model has no predictive power for the orbital period variations outside of the current TOA timeline.

9 Variation in the longitude of periastron, ω˙\dot{\omega}

We measured a marginally significant variation in the longitude of periastron $(\dot{\omega})$ for PSR J1017$-$7156 with a value of $0.022\pm 0.009\,°\rm\,yr^{-1}$. If we assume a typical pulsar mass of 1.4$\,M_{☉}$ and an orbital inclination of 60$\,°$, using Equation (2) of Weisberg & Taylor 1981 we obtain a predicted $\dot{\omega}$ in general relativity of 0.012$\,°\rm\,yr^{-1}$, which agrees with our measured value within 1.1-$\sigma$. In general $\dot{\omega}$ is a useful Post-Keplerian (PK) parameter as it can be used to calculate the total mass of the binary system, from which a measurement of the pulsar mass may be extracted. The variation in $\dot{\omega}$ is the easiest to measure for orbits with significant eccentricities. In the case of PSR J1017$-$7156 with $e=1.4\times 10^{-4}$ and an already good timing residual RMS of 1.3 $\rm\mu$s, we expect its $\dot{\omega}$ measurement to be much improved with another 5 years of timing data.

10 Gamma-ray pulsation searches

Among the pulsars in our sample, PSRs J1125$-$5825 and J1446$-$4701 have been observed to emit $>0.1$ GeV pulsations by Keith et al. 2012, through the analysis of data taken by the Fermi Large Area Telescope (Atwood et al. 2009, LAT;), with post-trial significances just under 5-$\sigma$. High confidence detections of these two MSPs in gamma rays were later presented in Abdo et al. 2013.

To determine whether other MSPs in our sample also emit gamma-ray pulsations, we analysed LAT photons recorded between 2008 August 4 and 2013 May 1, with energies from 0.1 to 100 GeV, and belonging to the ‘Source’ class of the reprocessed P7REP data, a version of Pass7 data See Bregeon et al. 2013 and http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Pass7REP_usage.html for more information. reprocessed with improved calibrations. Events with zenith angles larger than $100°{}$ were excluded, to reject atmospheric gamma rays from the Earth’s limb. In addition, events recorded when the instrument was not operating in nominal science operations mode, when the limb of the Earth infringed upon the regions of interest (see below for the definition of these regions), or when the data were not flagged as good were excluded. These cuts were made using the Fermi Science Tools http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/overview.html (STs) v9r32p5, and the selected photons were assigned pulse phases using the ephemerides listed in Tables 3 and 5 and the Fermi plug-in for tempo2 (Ray et al. 2011; Hobbs et al. 2006).

The spectral parameters listed in Table 8 for PSRs J1125$-$5825 and J1446$-$4701 are consistent with those reported in Abdo et al. 2013 to within uncertainties. The parameters for PSRs J1543$-$5149 and J1811$-$2405 are only weakly constrained at present, but are reminiscent of those of known gamma-ray MSPs (Abdo et al. 2013). Also listed in Table 8 are the gamma-ray luminosities $L_{\gamma}$ deduced from the energy flux measurements, and the efficiencies of conversion of spin-down power into gamma-ray emission, $\eta=L_{\gamma}/\dot{E}$, calculated using the Shklovskii-corrected $\dot{E}$ values and the DM distances given in Tables 3 and 5. The uncertainties reported in Table 8 are statistical. Studies of systematic uncertainties in the effective area suggest a 10% uncertainty at 100 MeV, decreasing linearly in Log(E) to 5% in the range between 316 MeV and 10 GeV and increasing linearly in Log(E) up to 10% at 1 TeV see http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html.

The two newly-identified gamma-ray pulsars, PSRs J1543$-$5149 and J1811$-$2405, bring the total number of MSPs with detected gamma-ray pulsations to 53 objects. It is unlikely that the Fermi LAT will detect many of the remaining MSPs presented in this paper. Assuming an average gamma-ray efficiency for MSPs of 0.245 following the prescription of Johnson et al. 2013, we derive expected energy fluxes for these pulsars much smaller than the lowest value reported in Abdo et al. 2013 for an MSP, because of the generally large distance values; with the notable exception of PSR J1731$-$1438. The latter MSP may be inefficient at converting its spin-down power into gamma-ray emission, or its gamma-ray beams may not cross the Earth’s line-of-sight. The high $\dot{E}$ but distant MSPs in this sample could contribute to the diffuse emission seen by the Fermi LAT around the Galactic plane.

Conclusion

The High Time Resolution Universe survey for pulsars and fast transients has discovered 27 MSPs to date, of which four are announced in this work. All four MSPs have phase-coherent timing solution with RMS already of the order of tens of $\mu$s. We have presented their pulse profiles and polarimetric properties at different frequencies. PSRs J1528$-$3828 and J1056$-$7117 are likely to be formed from wide-orbit LMXBs, leading to the formation of classic MSPs with He-WD companions. PSR J1755$-$3716 is likely to have evolved from an IMXB and possesses a CO-WD companion. PSR J1525$-$5545 is likely to have a massive CO-WD companion, or an ONeMg-WD companion if the orbital inclination angle is low.

In addition, we present updated timing solutions for 12 previously published HTRU MSPs, as compared with results in their respective discovery papers (Bates et al. 2011, Keith et al. 2012, Bailes et al. 2011), thanks to the now longer timing baseline of over three years in all cases, except one with 2.7 years of timing data.

We measure 5 new proper motions with significance greater than 3-$\sigma$, from PSRs 1017$-$7156, J1125$-$5825, J1446$-$4701, J1708$-$3506 and J1719$-$1438. Their derived transverse velocities are all consistent with previous MSP velocity distribution studies. In turn, with the proper motion measurements, we are able to constrain the period derivative contribution from the Shklovskii effect. In addition, we take into account the acceleration due to the Galactic potentials and correct for the intrinsic period derivatives for the 12 MSPs in this work. PSR J1017$-$7156 has one of the smallest inferred intrinsic period derivatives at 1.2$\times 10^{-21}$, hence also one of the lowest derived surface magnetic field strength within the known MSP population at a value of $5.4\times 10^{7}$ G.

We further discuss the case of PSR J1801$-$3210 for which no significant period derivative can be measured, even with more than 4 years of timing data. The best-fit solution in tempo2 shows a $\dot{P}_{\rm{obs}}$ of $-4\pm 4\times 10^{-23}$, an extremely small number comparing to that of a typical MSP. The both positive $\dot{P}_{\rm{shk}}$ and $\dot{P}_{\rm{gal}}$ of the order of $10^{-20}$ and $10^{-21}$, respectively, act to further decrease the already negative period derivative. It seems unlikely that the DM-derived distance is significantly wrong and hence reversing the direction of Galactic potential. Alternatively we consider the presence of a third body near PSR J1801$-$3210 which might be accelerating the pulsar towards Earth. Giant molecular clouds seem to be a plausible scenario, whereas an exoplanet orbiting in a large hierarchical orbit seems unlikely due to the small probability of surviving the SN, as well as the fact that we do not measure any significant second derivatives of spin frequency. Based on radio emission theory, we derive a theoretical lower limit of period derivative of $7.9\times 10^{-24}$ and a corresponding surface magnetic field strength at the equator of $7.8\times 10^{6}$ G for PSR J1801$-$3210, in order for it to stay above the pulsar death line. We also highlight the potential of PSR J1801$-$3210 to be employed in the SEP test due to its wide and circular orbit.

We have undertaken a comparison study between MSPs in our sample and the complete known pulsar population. We point to a strong dependence on inclination for eclipses to be observed in VLMBPs, as indicated by an apparent bimodal distribution of eclipsing and non-eclipsing systems separated by a companion mass of about $0.027\,M_{☉}$. We also suggest that the distribution of the total mass of binary systems is inversely-related to the Galactic height distribution. In other words, MSPs with the heaviest companions have larger tendencies to stay close to the Galactic plane, whereas lighter systems with smaller mass functions show larger mean value and larger scatter in the Galactic height distribution.

A change in the projected semi-major axis ($\dot{x}$) is observed in PSR J1017$-$7156 at $9.1\pm 1.7\times 10^{-15}$. Rather than due to gravitational wave emission, this $\dot{x}$ is likely due to an apparent change in the orbital inclination as a result of proper motion affecting the viewing geometry. We also report an $\dot{\omega}$ of 0.022(9)$\,°\rm\,yr^{-1}$, and we highlight the potential of measuring more relativistic orbital parameters with PSR J1017$-$7156. Together with its small period derivative and the corresponding low derived magnetic field as mentioned above, this makes PSR J1017$-$7156 a very interesting pulsar to be closely followed-up with further timing campaign, and indeed it is already being monitored by the Parkes Pulsar Timing Array Project (Manchester et al. 2013). Although we stress the importance of a careful dispersion measure variation treatment as discussed in Section 4.1 and a proper polarisation calibration to correctly assess the uncertainties on the TOAs for the high-precision timing required for PSR J1017$-$7156. Furthermore, orbital period variations are observed in the BW system PSR J1731$-$1847. We present the timing solution with the BTX timing model which demonstrate the quadratic changes in orbital period over the last three years.

We detected highly significant gamma-ray pulsations from PSRs J1125$-$5825 and J1446$-$4701, confirming the results of Keith et al. 2012 and Abdo et al. 2013. PSRs J1543$-$5149 and J1811$-$2405 were identified for the first time as gamma-ray pulsars: after folding the Fermi LAT photons with radio timing ephemerides, we obtained $>$ 5-$\sigma$ detections of these two MSPs, bringing the total number of MSPs with detected gamma-ray pulsations to 53 objects.

acknowledgements

The Parkes Observatory is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. A large amount of the timing data for PSR J1017$-$7156 has been taken as part of the Parkes Pulsar Timing Array Project P456.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France.

The authors would like to thank Patrick Lazarus for providing an automated RFI cleaning routine, Dominic Schnitzeler for useful discussions on pulsar distances, Joris Verbiest for his expertise on the use of tempo2, Charlotte Sobey and Gregory Desvignes for advices on data calibration, Kejia Lee for his help in plotting Figure 7, Chris Flynn for sharing his knowledge on Galactic stellar density, Megan DeCesar from the Fermi collaboration for reviewing the paper and providing many constructive suggestions, and Lijing Shao for carefully reading the manuscripts. CN was supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne.

References