A millisecond pulsar in an extremely wide binary system

A millisecond pulsar in an extremely wide binary system

C. G. Bassa, G. H. Janssen, B. W. Stappers, T. M. Tauris, T. Wevers, P. G. Jonker,L. Lentati, J. P. W. Verbiest, G. Desvignes, E. Graikou, L. Guillemot, P. C. C. Freire,P. Lazarus, R. N. Caballero, D. J. Champion, I. Cognard, A. Jessner, C. Jordan,R. Karuppusamy, M. Kramer, K. Lazaridis, K. J. Lee, K. Liu, A. G. Lyne,J. McKee, S. Osłowski, D. Perrodin, S. Sanidas, G. Shaifullah, R. Smits,G. Theureau, C. Tiburzi and W. W. Zhu
ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA, Dwingeloo, the Netherlands
Jodrell Bank Centre for Astrophysics, The University of Manchester, Manchester, M13 9PL, United Kingdom
Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany
Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany
Department of Astrophysics/IMAPP, Radboud University Nijmegen, P.O. box 9010, 6500 GL Nijmegen, The Netherlands
SRON, Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, the Netherlands
Institute of Astronomy / Battcock Centre for Astrophysics, University of Cambridge, Madingley Road, Cambridge CB3 0HA, United Kingdom
Fakultät für Physik, Universität Bielefeld, Postfach 100131, 33501 Bielefeld, Germany
Laboratoire de Physique et Chimie de l’Environnement et de l’Espace LPC2E CNRS-Université d’Orléans, F-45071 Orléans, France
Station de radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU F-18330 Nançay, France
Kavli institute for astronomy and astrophysics, Peking University, Beijing 100871, P.R. China
INAF - Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Laboratoire Univers et Théories, Observatoire de Paris, CNRS/INSU, Université Paris Diderot, 5 place Jules Janssen, 92190 Meudon, France
email: bassa@astron.nl
Accepted July 12, 2019. Received July 12, 2019; in original form July 12, 2019

We report on 22 yrs of radio timing observations of the millisecond pulsar J10240719 by the telescopes participating in the European Pulsar Timing Array (EPTA). These observations reveal a significant second derivative of the pulsar spin frequency and confirm the discrepancy between the parallax and Shklovskii distances that has been reported earlier. We also present optical astrometry, photometry and spectroscopy of 2MASS J102438690719190. We find that it is a low-metallicity main-sequence star (K7V spectral type, ,  K) and that its position, proper motion and distance are consistent with those of PSR J10240719. We conclude that PSR J10240719 and 2MASS J102438690719190 form a common proper motion pair and are gravitationally bound. The gravitational interaction between the main-sequence star and the pulsar accounts for the spin frequency derivatives, which in turn resolves the distance discrepancy. Our observations suggest that the pulsar and main-sequence star are in an extremely wide ( yr) orbit. Combining the radial velocity of the companion and proper motion of the pulsar, we find that the binary system has a high spatial velocity of  km s with respect to the local standard of rest and has a Galactic orbit consistent with halo objects. Since the observed main-sequence companion star cannot have recycled the pulsar to millisecond spin periods, an exotic formation scenario is required. We demonstrate that this extremely wide-orbit binary could have evolved from a triple system that underwent an asymmetric supernova explosion, though find that significant fine-tuning during the explosion is required. Finally, we discuss the implications of the long period orbit on the timing stability of PSR J10240719 in light of its inclusion in pulsar timing arrays.

stars: individual: PSR J10240719 – stars: neutron – binaries: general – superovae: general
pubyear: 2016pagerange: A millisecond pulsar in an extremely wide binary systemReferences

1 Introduction

The recent direct detection of gravitational waves by ground-based interferometers has opened a new window on the Universe (Abbott et al., 2016). Besides ground-based interferometers, gravitational waves are also predicted to be measurable using a pulsar timing array (PTA), in which an ensemble of radio millisecond pulsars are used as extremely stable, celestial clocks representing the arms of a Galactic gravitational wave detector (Detweiler, 1979; Hellings & Downs, 1983).

At present, three PTAs are in operation; the European Pulsar Timing Array in Europe (EPTA; Desvignes et al. 2016), the Parkes Pulsar Timing Array (PPTA) in Australia (Manchester et al., 2013) and NANOGrav in North-America (Demorest et al., 2013). The three PTAs have been in operation for about a decade, and have so far provided limits on the stochastic gravitational wave background at nano-hertz frequencies (e.g. Lentati et al. 2015; Shannon et al. 2013; Arzoumanian et al. 2015). The three PTAs collaborate together as the International Pulsar Timing Array (IPTA; Verbiest et al. 2016). As a result of improvements in instrumentation, analysis software, and higher observing cadence, the timing precision has increased with time, which leads previously unmodelled effects to become important.

The first data release by the EPTA (EPTA DR1.0) has recently been published by Desvignes et al. (2016) and provides high-precision pulse times-of-arrival (TOAs) and timing ephemerides for an ensemble of 42 radio millisecond pulsars spanning baselines of 7 to 18 yrs in time. While preparing this data release, we rediscovered the apparent discrepancy in the distance of the isolated millisecond pulsar J10240719 which was reported by Espinoza et al. (2013) and also noticed by Matthews et al. (2016). This discrepancy arises due to the high proper motion of the pulsar, which gives rise to an apparent positive radial acceleration and hence a change in the spin period known as the Shklovskii effect (Shklovskii, 1970). In the case of PSR J10240719, the Shklovskii effect exceeds the observed spin period derivative for distances larger than 0.43 kpc. Pulsar timing yields parallax distances beyond that limit, which would require an unphysical negative intrinsic spin period derivative for PSR J10240719.

This discrepancy could be resolved through a negative radial acceleration due to a previously unknown binary companion in a very wide orbit. In a search for non-thermal emission from isolated millisecond pulsars, Sutaria et al. (2003) studied the field of PSR J10240719 in optical bands. Their observations revealed the presence of two stars near the position of the pulsar. The broadband colors and spectrum of the bright star, 2MASS J102438690719190 (, hereafter star B) were consistent with those of a main-sequence star of spectral type K5, while the other, fainter star (, hereafter star F) had broadband colors which Sutaria et al. (2003) suggested could be consistent with non-thermal emission of the neutron star in PSR J10240719. However, the astrometric uncertainties in both the optical and radio did not allow Sutaria et al. (2003) to make a firm conclusion on the association of either stars B or F and PSR J10240719.

Here, we present radio observations of PSR J10240719 and optical observations of stars B and F (§2) that allow us to improve the timing ephemeris of PSR J10240719 and the astrometry, photometry and spectroscopy of stars B and/or F (§3). In §4 we present our results which show that star B is in an extremely wide orbit around PSR J10240719. An evolutionary formation scenario for this wide orbit is presented in § 5. We discuss our findings and conclude in §6. This research is the result of the common effort to directly detect gravitational waves using pulsar timing, known as the European Pulsar Timing Array (Desvignes et al., 2016)111http://www.epta.eu.org.

While preparing this paper, we became aware that a different group had used independent radio timing observations and largely independent optical observations to reach the same conclusion about the binarity of PSR J10240719 (Kaplan et al., 2016). The analysis and results presented in this paper agree very well with those of Kaplan et al. (2016), though these authors suggest an alternate formation scenario for PSR J10240719 (see §5.4).

2 Observations

2.1 Radio

We present radio timing observations of PSR J10240719 that were obtained with telescopes participating in the EPTA over the past 22 years. The EPTA DR1.0 data presented by Desvignes et al. (2016) contains timing observations of PSR J10240719 from the ’historical’ pulsar instrumentation at Effelsberg, Jodrell Bank, Nançay and Westerbork and contains 561 times-of-arrival (TOAs) over a 17.3 yr time span from January 1997 to April 2014.

The DR1.0 data set was extended back to March 1994 with observations obtained at 410, 610 and 1400 MHz with the Jodrell Bank analog filterbank (AFB; Shemar & Lyne 1996). The first year of this data set has been published in the PSR J10240719 discovery paper by Bailes et al. (1997). We also included observations obtained with the new generation of baseband recording/coherent dedispersing pulsar-timing instruments presently in operation within the EPTA. These instruments are PuMa II at Westerbork (Karuppusamy et al., 2008), PSRIX at Effelsberg (Lazarus et al., 2016), the ROACH at Jodrell Bank (Bassa et al., 2016) and BON and NUPPI at Nançay (Cognard & Theureau, 2006; Cognard et al., 2013). A subset of the Nançay TOAs on PSR J10240719 were presented in Guillemot et al. (2016). Combining these data sets yields 2249 TOAs spanning 22 yrs since the discovery of PSR J10240719 (Table 1 and Fig. 1).

2.2 Optical

We retrieved imaging observations of the field of PSR J10240719 from the ESO archive. These observations consist of  min -band,  min -band, and  min -band exposures, which were obtained with the FORS1 instrument (Appenzeller et al., 1998) at the ESO VLT at Cerro Paranal on March 28, 2001 under clear conditions with seeing. The high resolution collimator was used, providing a  pix pixel scale and a field-of-view. Photometric standards from the Landolt (1992) SA109 field were observed with the same filters and instrument setup. Analysis of this data has been presented by Sutaria et al. (2003).

The field of PSR J10240719 was also observed with OmegaCAM (Kuijken, 2011) at the VLT Survey Telescope. OmegaCAM consists of 32 4k2k CCDs with a  pix pixel scale. SDSS observations were obtained between February 24 and March 2, 2012 as part of the ATLAS survey (Shanks et al., 2015). Here we specifically use the 45 s SDSS -band exposure obtained on February 27, 2012. The seeing of that exposure was .

Follow-up observations of the field of PSR J10240719 were taken on June 10 and 11, 2015 with FORS2 at the ESO VLT at Cerro Paranal in Chile. A series of 8 dithered 5 s exposures in the -band filtered were obtained on June 10th, under clear conditions with mediocre seeing of . The standard resolution collimator was used with  binning, providing a pixel scale of  pix with a field-of-view.

Two long-slit spectra were also obtained with FORS2, one with the 600RI grism on June 10th and one with the 600z grism on June 11th. These grisms cover wavelength ranges of 5580 to 8310 Å and 7750 to 10430 Å, respectively. Both exposures were 1800 s in length and used the slit. The CCDs were read out using binning, providing a resolution of 6.3 Å, sampled at 1.60 Å pix for both the 600RI and 600z grism. The seeing during these exposures was and , respectively. The spectrophotometric standard LTT 3864 was observed with the same grisms and a slit. Arc-lamp exposures for all grism and slit combinations were obtained as part of the normal VLT calibration programme.

Fit and data-set
MJD range 49416.8—57437.0
Data span (yr) 21.96
Number of TOAs 2249
Rms timing residual (s) 2.81
Weighted fit Y
Reduced value 0.98
Reference epoch (MJD) 55000
Measured Quantities
Right ascension,
Proper motion in R.A., (mas yr)
Proper motion in decl., (mas yr)
Parallax, (mas) 0.77(11)
Pulse frequency, (s) 193.7156834485468(7)
First derivative of , (s)
Second derivative of , (s)
Third derivative of , ()
Fourth derivative of , ()
Dispersion measure, DM (pc cm) 6.4888(8)
First derivative of DM (pc cm yr)
Second derivative of DM (pc cm yr)
Clock correction procedure TT(BIPM2011)
Solar system ephemeris model DE421
Units TCB
Table 1: Parameters for PSR J10240719. All reported uncertainties have been multiplied by the square root of the reduced .

3 Analysis

Figure 1: Radio timing observations at different observing frequencies are shown in the top panel as a function of time. Observations from different observatories are shown in different colors. The middle panel shows timing residuals of PSR J10240719 based on our timing ephemeris (Table 1). The bottom panel shows the dispersion measure (DM) values as modelled by the power law model. The dashed lines indicate the uncertainty in the DM.

3.1 Pulsar timing

We use the standard approach, as detailed in Edwards et al. (2006), for converting the topocentric TOAs to the solar-system barycenter through the DE421 ephemeris (Folkner et al., 2009) and placing them on the Terrestial Time standard (BIPM2011; Petit 2010). The TOAs from the different combinations of telescope and instrument were combined in the Tempo2 pulsar timing software package (Hobbs et al., 2006).

The best-fit timing solution to the TOAs was found through the standard Tempo2 iterative least-squares minimization. The free parameters in the fit consist of the astrometric parameters (celestial position, proper motion and parallax), a polynomial describing the spin frequency as a function of time, and a polynomial to describe variations in the pulsar dispersion measure (DM) as a function of time. We also fitted for offsets between the telescope and instrument combinations to take into account difference in pulse templates and instrument delays.

To perform model selection between models that include increasing numbers of spin frequency derivatives when incorporating increasingly complex noise models, we use the Bayesian pulsar timing package TempoNest (Lentati et al., 2014), which operates as a plugin to Tempo2. In all cases we include parameters to model the white noise that scale, or add in quadrature with, the formal TOA uncertainties, and define these parameters for each of the 20 observing systems included in the data set. We also include combinations of additional time correlated processes, including DM variations and system-dependent noise. For this final term we use the approach described in Lentati et al. (2016). Apart from the spin frequency, derivatives we marginalize over the timing model analytically. The spin frequency derivatives are included in the analysis numerically with priors that are uniform in the amplitude of the parameter in all cases.

When including only white noise parameters, and limiting our model for DM variations to a quadratic in time, we find the evidence supports a total of four frequency-derivatives. However, we find there is significant support from the data for higher order DM variations. In particular, the difference in the log evidence for models with and without an additional power law DM model besides the quadratic model is 36 in favor of including the additional stochastic term. We find the spectral index of the DM variations () to be consistent with that expected for a Kolmogorov turbulent medium. When including these additional parameters to describe higher-order DM variations, we find that the data only support the first and second frequency-derivatives. The significant detections of the third and fourth frequency derivatives observed in the simpler model can therefore not be confidently interpreted as arising from a binary orbit, as they are highly covariant with the power-law variations in DM, thereby reducing their significance. Finally, we find no support for any system dependent time-correlated noise in this data set. The timing ephemeris using a DM power law model is tabulated in Table 1 and timing residuals and DM values as a function of time are plotted in Fig. 1.

3.2 Photometry

The FORS1 and FORS2 images were bias-corrected and flatfielded using sky flats. Instrumental magnitudes were determined through point spread function (PSF) fitting with the DAOPHOT II package (Stetson, 1987). The instrumental PSF magnitudes of the 2001 observations were calibrated against 17 photometric standards from the SA109 field, using calibrated values from Stetson (2000) and fitting for zeropoint offsets and color terms, but using the tabulated ESO extinction coefficients of 0.113, 0.109 and 0.087 mag airmass for , respectively. The residuals of the fit were 0.03 mag in , 0.06 mag in and 0.07 mag in . We find that star B has , and , while star F has , and .

The instrumental magnitudes of the 8 FORS2 -band images from 2015 were calibrated against the 2001 FORS1 observations. As the FORS2 -band (R_SPECIAL) filter has a different transmission curve compared to the FORS1 Bessell -band (R_BESS) filter, we fitted -band offsets as a function of the color, finding a color term of and residuals of 0.05 mag. We find that star B has , i.e. consistent with the FORS1 measurement.

3.3 Astrometry

We use the OmegaCAM -band exposure to obtain absolute astrometry of the field of PSR J10240719 as the individual CCDs from OmegaCAM have the largest field-of-view each (). The image containing the field of PSR J10240719 was calibrated against astrometric standards from the fourth U.S. Naval Observatory CCD Astrograph Catalog (UCAC4, Zacharias et al. 2013), which provides proper motions through comparison with other catalogs. A total of 14 UCAC4 astrometric standards overlapped with the image. The UCAC4 positions were corrected for proper motion from the reference epoch (2000.0) to the epoch of the OmegaCAM image (2012.15). The centroids of these stars were measured and used to compute an astrometric calibration fitting for offset, scale and position angle. Rejecting two outliers with residuals larger than , the calibration has rms residuals of in right ascension and in declination.

The position of star B on the OmegaCAM -band image is , , where the uncertainties are the quadratic sum of the uncertainty in the astrometric calibration and the positional uncertainty on the image ( in both coordinates). Star F is not detected on the OmegaCAM images. Hence, we transferred the OmegaCAM calibration to the median combined FORS1 -band image using 7 stars in common to both (excluding star B). The residuals of this transformation are in both coordinates. Star F is at , . Here, the uncertainties also include the uncertainty in the calibration between the OmegaCAM and FORS1 image.

The 14.2 yr time baseline between the FORS1 and FORS2 images allows us to determine the proper motions of stars B and F. To prevent pollution by non-random proper motions of stellar objects, we determined the absolute proper motion with respect to back ground galaxies. We selected 7 objects on the FORS1 and FORS2 images, located within of PSR J10240719, which were clearly extended in comparison to the PSF of stars. The 5 s FORS2 -band exposures were not deep enough to accurately measure the positions of the objects, so instead we used the 40 s FORS2 acquisition exposure for the long-slit spectroscopic observations. This image was taken with the OG590 order sorting filter, which cuts light blue-wards of 6000 Å.

Using the 7 extra-galactic objects, we determined the transformation between each of the three FORS1 -band images and the FORS2 OG590 image. The residuals of the transformations were typically in each coordinate. The transformations were then used to compute the positional offsets of stars B and F between the FORS1 and FORS2 images. Averaging the positional offsets and taking into account the 14.2 yr time baseline yielded the proper motion of both stars. We find that star B had  yr and  yr, while star F has  yr and  yr. Table 2 lists the magnitudes, position and proper motion of stars B and F.

Object (yr) (yr) ref.
PSR J10240719 1
PSR J10240719 2
PSR J10240719 3
PSR J10240719 4
Star B (PPMXL) 5
Star B (APOP) 6
Star B 1
Star F 1

References: (1) this work; (2) Desvignes et al. (2016); (3) Reardon et al. (2016); (4) Matthews et al. (2016);
(5) Roeser et al. (2010); (6) Qi et al. (2015)

Table 2: The positions and proper motions of stars B and F and PSR J10240719 are referenced to epoch MJD 55984.135. The uncertainties in the magnitudes of stars B and F are instrumental. The zeropoint uncertainties of 0.03 mag in , 0.06 mag in and 0.07 mag in should be added in quadrature to obtain absolute uncertainties.

3.4 Spectroscopy

The spectroscopic observations were bias corrected and flatfielded using lamp-flats. Spectra of the counterpart and the spectrophotometric standard were extracted using the optimal extraction method by Horne (1986), and wavelength calibrated using the arc-lamp exposures. The residuals of the wavelength calibration were  Å or better. The extracted instrumental fluxes were corrected for slit losses by using the wavelength dependent spatial profile to estimate the fraction of the flux being masked by the finite slit width (about 28% for the 1 slit and 1% for the 5 slit). The instrumental response of both instrument setups was determined from the standard observations and a comparison with calibrated spectra from Hamuy et al. (1992, 1994). The response was then used to flux calibrate the spectra of star B.

4 Results

4.1 Pulsar timing

As the data set presented here is an extension of the EPTA DR1.0 data set, the parameters in Table 1 have improved with respect to the Desvignes et al. (2016) ephemeris. The parameters are generally consistent with the exception of the polynomial describing the DM variations, primarily due to the inclusion of the power law model. The inclusion of the multi-frequency observations from the Jodrell Bank AFB and PuMa II at WSRT have significantly improved the measurement of DM variations.

Our timing ephemeris is also generally consistent with the PPTA and NANOGrav results. In Table 2 we list the positions and proper motion of PSR J10240719 from the timing ephemerides by Desvignes et al. (2016); Reardon et al. (2016) and Matthews et al. (2016) propagated to the epoch of the OmegaCAM observations. We find that the positions and proper motions are consistent within the uncertainties. The observed parallax  mas corresponds to a distance of  kpc after correcting for Lutz-Kelker bias (Lutz & Kelker, 1973; Verbiest et al., 2010), and is consistent with the values found by Desvignes et al. (2016); Reardon et al. (2016); Matthews et al. (2016).

The improvements in the timing ephemeris allow us to measure a significant second derivative of the spin frequency ; . This higher order derivative likely accounts for the steep red noise spectrum seen in the PPTA observations of PSR J10240719 (Reardon et al., 2016). Our measurement of translates to  s. For comparison, Matthews et al. (2016) use the NANOGrav data and the measurement by Verbiest et al. (2009) to set a limit on  s. Guillemot et al. (2016) measured  s from Nançay observations of PSR J10240719. This value is inconsistent with our value; we attribute this difference to unmodelled DM variations in the Guillemot et al. (2016) timing solution.

The observed spin frequency derivative contains contributions from the intrinsic pulsar spindown, the Shklovskii effect due to non-zero proper motion (Shklovskii, 1970), differential Galactic rotation and Galactic acceleration (e.g. Nice & Taylor 1995) and any contribution due to radial acceleration. These effects are more commonly described as derivatives of the spin period such that .

Figure 2: The top panel shows the spin period derivative due to the non-zero proper motion (the Shklovskii effect) as a function of distance as the dashed line. For distances larger than 0.43 kpc, it exceeds the observed spin period derivative (solid horizontal line). The histogram on the right-hand side of the top panel shows the distribution of intrinsic spin period derivatives for all radio millisecond pulsars in the Galactic field (Manchester et al., 2005). The bottom panel shows the radial acceleration required to explain the observed spin period derivative of PSR J10240719 given the Shklovskii contribution assuming four different values for the intrinsic spin period derivative of the pulsar. The solid and dashed vertical lines denote the Lutz & Kelker (1973) bias corrected distance of PSR J10240719 derived from our parallax measurement.

Figure 2 shows the observed spin period derivative and the Shklovskii contribution as a function of distance. As noted by Desvignes et al. (2016), Matthews et al. (2016) and Guillemot et al. (2016), the large proper motion of PSR J10240719 leads the Shklovskii contribution to exceed the observed spin period derivative  s s for distances larger than 0.43 kpc. Our results confirm this. The histogram on the right of Fig. 2 shows the distribution of intrinsic spin period derivatives of known millisecond pulsars ( ms; Manchester et al. 2005). These range between and  s s, and PSR J10240719 is expected to have a in this range. Desvignes et al. (2016) and Matthews et al. (2016) estimate that for PSR J10240719 the contributions due to differential Galactic rotation and Galactic acceleration only account for  s s, which is negligible with respect to the other terms making up . Depending on the value of , the remaining spin period derivative due to radial acceleration required to explain the observed spin period derivative sets the radial acceleration between and  cm s.

4.2 Association with PSR J10240719

In Table 2, we list the proper motion and pulsar position propagated to the epoch of the OmegaCAM observations (MJD 55984.135) for several timing ephemerides of PSR J10240719. Within the uncertainties, the pulsar position and proper motion as measured by the three pulsar timing arrays are consistent with each other. We find that star B is offset from PSR J10240719 by and , corresponding to a total offset of . Here, the uncertainty is dominated by the astrometric calibration against the UCAC4 catalog. For star F, the total offset is . Moreover, the proper motions determined from the 14.2 year baseline between the FORS1 and FORS2 observations show that, within , star B has a proper motion consistent with the pulsar.

Figure 3: The position of PSR J10240719 (black), star B (blue) and star F (grey) between 1950 and 2025. The positional uncertainties of stars B and F are indicated by the ellipses at the OmegaCAM observation of star B (epoch 2012.15) and the FORS1 observation of star F (epoch 2001.24). The lines expanding from the ellipse illustrate the positional uncertainty due to position and proper motion as a function of time. On the scale of this plot, the uncertainties on the timing position and proper motion of PSR J10240719 are negligible. Also plotted are measurements of the position of star B in the USNO-A2 (Monet et al., 1998), USNO-B1 (Monet et al., 2003), 2MASS (Cutri et al., 2003; Skrutskie et al., 2006) and PPMXL (Roeser et al., 2010) catalogs. The epochs of these measurements are given in brackets. The USNO-A2 catalog does not provide positional uncertainties; we conservatively estimate uncertainties on the position.
Figure 4: Proper motion measurements of PSR J10240719 and star B. The proper motion of PSR J10240719 is denoted by the thick black dot. The proper motion of star B as determined in this work, the PPMXL catalog (Roeser et al., 2010) and the APOP catalog (Qi et al., 2015) are shown with the triangle, circle and square, respectively. The small points and histograms at the top and right of the figure represent proper motion measurements from stars in the APOP catalog, selecting 14754 stars within a radius of around PSR J10240719.

The similarity in position and proper motion between star B and PSR J10240719 is independently confirmed by several surveys. At , star B is bright enough to have been recorded on historic photographic plates and hence is included in the USNO-A2 (Monet et al., 1998) and USNO-B1 (Monet et al., 2003) astrometric catalogs. It is also detected in the near-IR in the 2MASS survey (Cutri et al., 2003; Skrutskie et al., 2006). Star B is present in the PPMXL catalog by Roeser et al. (2010), which combines the USNO-B1 and 2MASS astrometry to determine proper motions, and also the APOP catalog by Qi et al. (2015), which uses STScI digitized Schmidt survey plates originally utilized for the creation of the GSC II catalog (Lasker et al., 2008) to derive absolute proper motions. The position and proper motion of star B in these catalogs is plotted in Fig. 3 and Fig. 4 and listed in Table 2.

The probability that an unrelated star has, within the uncertainties, a position and proper motion consistent with PSR J10240719 is minuscule. The APOP catalog and the FORS2 photometry yield a stellar density for stars with (equal or brighter than star B) of 1 to 2 stars per square arcminute, while within from PSR J10240719, only 8 out of 7300 APOP stars with have a proper motion in right ascension and declination that is within 10 mas yr of that of PSR J10240719 ( mas yr). Based on these numbers, we estimate that the chance probability of a star having a similar position and proper motion to PSR J10240719 is about for stars equal or brighter than star B. Such a low probability confirms that star B is associated with PSR J10240719 and that both objects form a common proper motion pair.

At or above the brightness level of star F, there are about 10 objects per square arcminute in the FORS1 -band image, suggesting that there is a probability of about 7% of finding an object as bright and close as star F with respect to the pulsar position. Hence, we consider star F as a field star, not related to PSR J10240719.

Figure 5: The FORS2 optical spectrum of star B, formed from the combination of the single exposures obtained with the 600RI and 600z grisms, is shown in black. The spectrum has been shifted to zero velocity. The and -band magnitudes have been converted to fluxes using the zeropoints of Bessell et al. (1998) and are plotted. Regions affected by telluric absorption have been masked. The most prominent spectral lines are indicated. Two synthetic spectra by Munari et al. (2005) are shown which bracket the best fit model to the observed spectrum. Both model spectra have and has been rotationally broadened by 83 km s and convolved by a truncated Moffat profile representing the observational seeing and slit width. The top light-grey model spectrum has  K and  cgs, while the dark grey model spectrum has  K and  cgs. The models have been reddened using the extinction law by Cardelli et al. (1989) and and then scaled to match the observed spectrum. The model spectra are offset by 0.3 and 0.4 vertical units for clarity.

4.3 Properties of star B

The spectroscopic observation of star B shows strong absorption lines of Na D and Ca II, while H is weak. There is some suggestion of absorption from the TiO bands near 6300 and 7000 Å. A comparison against templates from the library by Le Borgne et al. (2003) favors a K7V spectral type over K2V or M0V templates. This classification agrees with the previous work by Sutaria et al. (2003), who classified star B as a K-type dwarf.

To obtain the radial velocity and spectral parameters of star B, we compare the observed spectrum (Fig. 5) with synthetic spectral templates from the spectral library of Munari et al. (2005). The templates from this library cover a range in effective temperature , surface gravity , rotational velocity and metallicity . To account for the resolution of the observed spectrum, we convolved the synthetic spectra with a Moffat profile that has a width equal to the seeing and is truncated at the width of the slit. For each combination of temperature, surface gravity and metallicity, we used the convolved templates to find the best fitting values for rotational velocity and radial velocity . These convolved, broadened and shifted spectra were then compared with the observed spectrum while fitting the ratio of observed and normalized flux with a second order polynomial.

To exclude possible pollution of the observed spectrum by telluric absorption, the wavelength ranges between 5600 and 6000 Å and 8400 and 8800 Å were selected for the comparison. These ranges contain the strong lines of the Na D doublet and the Ca II triplet, but also the weak H line. The observed spectrum is best represented by the models with , at an effective temperature of  K and surface gravity  cgs, broadened by  km s at a radial velocity with respect to the solar-system barycenter of  km s. The reduced for 875 degrees-of-freedom.

4.4 Mass, distance and reddening

The effective temperature and surface gravity of star B are consistent with those of a low-mass main-sequence star, as higher mass stars will have evolved off the main-sequence and have much lower surface gravities ( cgs).

In Fig.6 we plot the predicted effective temperatures of the PARSEC stellar evolution models of by Bressan et al. (2012), Tang et al. (2014) and Chen et al. (2014, 2015) to estimate the mass and distance of star B. Here we use our Johnson-Cousins photometry as well as the SDSS photometry from the ATLAS survey (Shanks et al., 2015), where star B has , , and . We further more use the model by Green et al. (2014, 2015) to obtain the reddening as a function of distance towards PSR J10240719. This model predicts that the reddening reaches a maximum value of for distances larger than 800 pc. Combined with the extinction coefficients for the and filters from Schlafly & Finkbeiner (2011), the observed magnitudes are corrected for absorption and compared with the predicted absolute magnitudes from the PARSEC models to obtain the plotted distance modulus and hence distance .

Figure 6: Mass and distance predictions from the PARSEC stellar evolution models by Bressan et al. (2012); Tang et al. (2014); Chen et al. (2014, 2015) at the observed temperature and broadband and magnitudes. The observed magnitudes were corrected for reddening using from the model by Green et al. (2014, 2015) and the Schlafly & Finkbeiner (2011) extinction coefficients. The error bar on the left denotes the Lutz & Kelker (1973) bias corrected parallax distance of PSR J10240719.

We find that models with masses between 0.43 and 0.46 M can explain the observed temperature. Taking into account the scatter in the distance moduli between the different filters, our dereddened photometry constrains the distance to 1.1 to 1.4 kpc. These distances are consistent with the parallax distance (corrected for Lutz & Kelker 1973 bias) of  kpc for PSR J10240719 determined from pulsar timing, further confirming the association of star B with PSR J10240719.

4.5 Orbit constraints

The question now remains if star B can account for the apparent radial acceleration and higher order spin frequency derivatives seen in the timing of PSR J10240719.

Since it is unlikely that star B is in an unbound orbit around PSR J10240719, we will consider a bound orbit in the following analysis. In Appendix A, we derive expressions for the line-of-sight acceleration , jerk , snap and crackle acting on the pulsar due to orbital motion and also show they relate to the spin frequency derivatives. This derivation is similar to that presented in Joshi & Rasio (1997), though we also provide expressions for the projected separation of the binary members on the sky, as our optical and timing astrometry provides an additional constraint.

Besides the line-of-sight acceleration , the observed second derivative of the spin frequency translates into a jerk of  cm s along the line-of-sight. The limits on and constrain line-of-sight snap and crackle to  cm s and  cm s, respectively. Furthermore, the angular separation of between the pulsar and star B at the parallax distance of  kpc corresponds to a projected separation  AU.

As shown in Appendix A, six parameters are required to describe the orbit; orbital period , semi-major axis , inclination , eccentricity , argument of perigee and true anomaly . By treating the companion masses as known, i.e. a pulsar mass of  M and a companion mass of  M, we can relate the orbital period to the semi-major axis , and also to the semi-major axis of the pulsar orbit . To investigate which orbits are allowed by the observations, we performed a Monte Carlo simulation where random values for , and were drawn from uniform distributions between , and . Furthermore, we choose random values of from a uniform distribution between and  cm s to account for the unknown intrinsic spin period derivative of PSR J10240719. These values were used to compute and such that the equations for and can then be divided to obtain and . Based on these parameters, values for , and were computed and compared with the observational constraints. We consider only those orbits that yield values within of the observed values or limits.

Figure 7 shows the results of the Monte Carlo simulation. A total of 10000 possible orbits are depicted as a function of orbital period. Two families of solutions are found. The solutions with are eccentric () and have orbital periods up to 2000 yr and favor and , i.e. placing the pulsar at apastron and the apsides pointing towards the observer. The second family of solutions is less constrained and allows for essentially all eccentricities, also circular orbits, but requires the orbital period to be longer than 6000 yr and . Depending on the eccentricity and orbital period, all values for and are possible.

Figure 7: The results of a Monte Carlo simulation to find orbits which account for the observed properties of PSR J10240719. Random values for , and are drawn from uniform distributions and used to compute and . Only those orbits which predict , and consistent with the measurements or limits are plotted.
Figure 8: The present Galactic position of PSR J10240719 is indicated with the black dot, while the line traces the Galactic orbit of the system over the next 5 Gyr. The location of the Sun at a Galactic radius of  kpc is denoted with , while the Galactic center is at the plus sign. It is clear that PSR J10240719 is in an eccentric and inclined orbit in the Galaxy.

4.6 Space velocity

The radial velocity measurement of star B of  km s complements the transverse velocity  km s from the proper motion measurement of PSR J10240719 and allows us to determine the 3D velocity of the system in the Galaxy. Based on the orbital solutions presented in § 4.5, we find that the corrections on the transverse velocity of PSR J10240719 and the radial velocity of star B with respect to the binary center of mass are negligible in comparison with the uncertainties on and . As such, we treat and as representative of the velocity of the binary center of mass, yielding a total velocity of  km s with respect to the solar-system barycenter. Converting these velocities into the Galactic frame using the algorithm by Johnson & Soderblom (1987)222With a rotation matrix for the J2000 input frame. and correcting for solar motion with the values from Hogg et al. (2005), the resulting peculiar motion with respect to the LSR is  km s.

The distance, location and proper motion of PSR J10240719 place it 0.82 kpc above and moving towards the Galactic plane. Using the MWPotential2014 Galactic potential and the orbit integrator implemented in the Galpy software package by Bovy (2015), we numerically integrate the orbit of PSR J10240719 in the Galaxy. The result of the integration is plotted in Fig. 8, which shows that the Galactic orbit of PSR J10240719 is eccentric and inclined. The Galactic orbit varies between Galactic radii of 4 and 9 kpc () and reaches up to 5.8 kpc above the Galactic plane. The vertical velocity of the system as it crosses the Galactic plane varies between between 120 and 220 km s. With these properties, PSR J10240719 can be classified as belonging to the Galactic halo (Bland-Hawthorn & Gerhard, 2016).

5 A triple star formation scenario

The standard formation scenario for millisecond pulsars is not applicable in the case of PSR J10240719. In that scenario, a period of mass and angular momentum transfer from a stellar binary companion spins up an old neutron star (Alpar et al., 1982; Bhattacharya & van den Heuvel, 1991). The outcome of that evolutionary channel is a “recycled” radio millisecond pulsar orbiting the white dwarf remnant of the main-sequence star in a short ( d) and circular orbit (e.g. Tauris 2011; Tauris et al. 2012). This channel has been confirmed observationally by the identification of several white dwarf binary companions to MSPs (e.g. van Kerkwijk et al. 2005; Bassa et al. 2006; Antoniadis et al. 2012; Kaplan et al. 2013).

In the case of PSR J10240719, the current binary companion is a main-sequence star and thus can not have recycled the pulsar. Currently we know of only one other MSP with a main-sequence companion; PSR J1903+0327, which is a 2.15 ms period MSP in a 95 day eccentric () orbit (Champion et al., 2008) around a  M G-dwarf companion (Freire et al., 2011; Khargharia et al., 2012). This system is believed to have evolved from a hierarchical triple consisting of a compact low-mass X-ray binary, accompanied by the G dwarf in a wide orbit. The inner companion then was either evaporated by the pulsar (Freire et al., 2011), or the widening of the orbit of the X-ray binary due to mass transfer led to a dynamical instability which ejected the inner companion from the triple (Freire et al., 2011; Portegies Zwart et al., 2011).

For PSR J10240719, we argue below that the ejection scenario is unlikely based on its orbital constraints. The evaporation scenario proposed in Freire et al. (2011) is only described qualitatively. To obtain a quantitative description, we investigate if PSR J10240719 could have evolved from a triple system in which the inner companion was evaporated. Further considerations of mass transfer, stability, and the effects of supernovae in triple systems can be found in Tauris & van den Heuvel (2014).

5.1 Summary of our model

There are a number of criteria which must be fulfilled in order to produce a stable triple system which would eventually leave behind PSR J10240719: i) the triple system must remain bound after the supernova (SN) explosion which creates the neutron star (NS), ii) to survive the SN, the secondary and the tertiary star must therefore be brought close to the exploding star via significant orbital angular momentum loss, e.g. in a common envelope, iii) the SN explosion must result in a large 3D systemic velocity and cause the tertiary star to be almost (but not quite) ejected, in order to explain the observed properties of PSR J10240719, iv) the post-SN triple system must remain dynamically stable on a long timescale (i.e. avoid chaotic three-body interactions) and v) the inner binary must evolve to evaporate the secondary star after the LMXB recycling phase, leaving behind PSR J10240719.

Figure 9: Illustration of the formation of PSR J10240719 via a triple star scenario: 1) A triple system is formed with a massive star (an O-star, the progenitor of the neutron star, NS), an F-star (needed to recycle the NS) and a K-dwarf (star B as presently observed). 2) When the massive star becomes a giant it expands and engulfs its two orbiting stars. 3) The supernova (SN) explosion. 4) The dynamically stable post-SN system. 5) Mass-transfer and recycling of the NS during the low-mass X-ray binary (LMXB) phase of the inner binary. 6) Evaporation of the remnant of the secondary star from the energetic wind of the recycled millisecond pulsar. 7) The present day system PSR J10240719. See text for a discussion of the various stages.

We find that the following model is able to account for the criteria listed above, albeit with significant fine-tuning during the SN explosion. The model is illustrated in Fig. 9. Starting on the zero-age main-sequence (ZAMS), the system consists of a roughly 16 to 18 M primary O-star and two low-mass companions: a secondary F-star with a mass of to 1.5 M and a tertiary K-dwarf with a mass  M, the presently observed star B (stage 1). After a common envelope phase (stage 2), where the extended envelope of the primary engulfed the other two stars prior to its ejection, the orbital period of the inner binary, was reduced to a few days and the orbital period of the tertiary star was  d. Following a phase of wind-mass loss333Given the low metallicity of star B, we only expect a moderate amount of wind-mass loss from the helium star prior to its collapse. from the exposed (initially)  M helium core, that star reaches a mass of  M (and  d) when it collapses and undergoes a SN explosion (stage 3). Whereas the inner binary remains compact (but now with a fairly high eccentricity), the dynamical effects of the SN causes the tertiary star to be almost (but not quite) unbound, reaching an orbital period of to  yr in an extremely eccentric orbit (), and yet remaining in a long-term dynamically stable three-body system (stage 4, see Section 5.3). A combination of tidal damping and magnetic braking eventually causes the secondary star to fill its Roche-lobe while still on the main sequence and initiate mass transfer to the NS, i.e. the inner binary becomes an LMXB system (stage 5). The LMXB system is converging, i.e. evolving to a small orbital period, and the secondary star may either form a low-mass helium white dwarf companion with a mass of  M (Istrate et al., 2014) or continue evolving through the period minimum and become a semi-degenerate ”black widow”-type companion with a mass of a few  M (Podsiadlowski et al., 2002; Chen et al., 2013). Meanwhile, the accreting NS is recycled and turns on as an energetic millisecond pulsar which will evaporate its companion star completely (stage 6), thereby leaving behind a system with the presently observed properties of PSR J10240719 (stage 7).

5.2 The dynamical effects of the SN explosion

We now describe in more detail the dynamically important properties of our model. The dynamical effects of asymmetric SNe in binaries, leading to either bound or disrupted systems, has been studied in full detail (Flannery & van den Heuvel, 1975; Hills, 1983; Tauris & Takens, 1998). A general discussion of dynamical effects of asymmetric SNe in hierarchical multiple star systems is found in Pijloo et al. (2012), and references therein. Here, we follow their method and apply a two-step process, where we first calculate the effects of an asymmetric SN explosion in the inner binary and then treat the inner binary as an effective point mass with respect to the tertiary star and calculate the dynamical effects of the SN on the outer binary. Hence, we use the obtained recoil velocity of the inner binary as a “kick” velocity added to this effective point-mass star representing the inner binary and then calculate the solution for the outer binary. We assume the SN explosion is instantaneous and for simplicity we neglect in our example any dependence on the orbital phase or inclination of the inner binary with respect to the outer binary.

Assuming a circular pre-SN orbit, the change of the binary semi-major axis as a result of an asymmetric SN, is given by (Hills, 1983):


where is the pre-SN semi-major axis (radius), the post-SN semi-major axis, the effective mass loss during the SN (here in our example, we shall assume 2.9 M to be lost when the exploding star leaves behind a NS with a gravitational mass of  M), the pre-SN total mass of the binary, the pre-SN orbital velocity of the exploding star in a reference fixed on the secondary companion star, the magnitude of the kick velocity, and the angle () between the kick velocity vector, and the pre-SN orbital velocity vector, . When calculating the dynamical effects for the outer binary, we simply use Eq. (1) with and .

In the case of a purely symmetric SN (), a binary will always disrupt if . This follows from the virial theorem and is also seen when the denominator in Eq. (1) becomes negative). Since the tertiary star (star B) is almost ejected from the system, we can estimate a critical mass of the exploding star to be 4.55 M (for which , if we assume  M). In our selected model, we then chose a slightly smaller mass of  M since we also need to apply a kick to obtain a large systemic velocity.

From Eq. (1), we can find a critical value for the kick angle, such that a binary will always dissociate if . Hence, by integrating over all kick angles, assuming an isotropic distribution of kick directions, we can easily calculate the probability ) for the binary to remain bound (Hills, 1983):

Figure 10: Probability for a triple system to survive a given recoil velocity, obtained by the inner binary due to a SN. The three lines are for a pre-SN orbital period of the tertiary star of , 115 and 3000 d. The bullet points represent average values of for the stated kick magnitudes (in km s) imparted onto the newborn 1.4 M NS in the inner binary, calculated from Monte Carlo simulation of SNe and using an isotropic distribution of kick directions for both and . When calculating , only triple systems for which the inner binary avoided merging were considered. However, the long-term stability of the surviving triple systems is not included here. The assumed pre-SN mass of the exploding star is  M, the secondary star has a mass of  M and the tertiary star has a mass of  M. The pre-SN orbital period of the inner binary is  d. The red star at  km s ( d) marks the example investigated further in this scenario.

In Fig. 10, we plot examples of the probability444The plotted probabilities do not take into account the specific requirements on the value of the post-SN orbital period of the tertiary star necessary for forming PSR J10240719. If including this specific constraint, the probabilities shown would be much lower. Similarly, for estimating the resulting values of we did not filter out those surviving triples which would not survive the SN in long-term stable orbits (see Section 5.3). of the outer binary to remain bound as a function of the resulting recoil velocity, of the inner binary, for three different pre-SN values of the orbital period of the tertiary star (, 115 and 3000 days). The dots indicate average values of for a given kick velocity added to the newborn NS (between 50 and 350 km s). To calculate each of these values, we simulated SN explosions using Monte Carlo techniques. As a result of mass loss of the inner binary, the minimum value of for any bound solution of the triple system is 67.8 km s, thus excluding any solutions in the grey-shaded (forbidden) region in Fig. 10. It is clearly seen that the probability of remaining bound is small for large values of (e.g. there are no possible solutions if  d). The star at  km s is for the selected example investigated further in this scenario (Figs. 11 and 12).

The values of the pre-SN orbital periods of the tertiary star () are limited to an interval where it cannot be too large (otherwise the system would dissociate) or too small (in which case the pre-SN triple system would not be stable due to three-body interactions between the stars. Besides reproducing the orbital period of the K-dwarf (star B) in PSR J10240719, we must also make sure to reproduce the large systemic velocity, of PSR J10240719. As derived in Section 4.6, the vertical velocity of PSR J10240719 when it crosses the Galactic plane must be between 120 and 220 km s. Hence, our model must produce at least  km s for an optimal orientation of . A large value of requires a large value of , which depends on the magnitude of the kick, imparted onto the newborn NS during the SN as well as the value of the pre-SN orbital period of the tertiary star .

Figure 11: Post-SN orbital period of the inner binary () as a function of its resulting recoil velocity () and kick velocity () imparted onto the NS. Each colored track is for a given value of (between 50 and 450 km s) and plotted along the tracks are solutions for different values of the kick angle (in steps of ). Grey circles represent systems that merge during the SN. The open star indicates the selected example which we investigate in more detail.

In Fig. 11 we demonstrate how increases with the kick velocity, . For these calculations, we assumed in all cases a pre-SN inner binary with circular orbit and  M,  M,  M,  d and varied between 50 and 450 km s ( is not included as it would lead to dissociation of the system). Each point (colored circle) represents a value of the kick angle, (in steps of 1°) for a given value of . The magnitude of is given at the bottom of each track which ends at . For clarity, we assumed here that the kick was in the pre-SN plane of the inner binary – i.e. the second kick angle or , where is the angle between the projection of the kick velocity vector onto a plane perpendicular to the pre-SN velocity vector of the exploding helium star and the pre-SN orbital plane (e.g. Fig. 1 in Tauris & Takens, 1998). The post-SN orbital period of the inner binary, is seen to decrease monotonically with larger values of , i.e. when the kick is directed backwards compared to the pre-SN orbital motion of the exploding star. Similarly, we have that for (cf. Eq. 1). Grey colors mark systems which immediately merge after the SN, i.e. systems where the periastron separation of the post-SN inner binary is smaller than the radius ( R) of the secondary star (here assuming a  M main-sequence F-star). The secondary star in our model must have a mass between 1.1 and 1.5 M to be massive enough to evolve within a Hubble time and still be light enough to possess a convective envelope, which is necessary for magnetic braking to operate during the LMXB phase.

The open black star symbol in Fig. 11 indicates our selected example model (see below) for which  km s and , leading to post-SN values of  d, and  km s. The surviving systems for which  km s produce too large post-SN eccentricities to ensure a long-term stable triple system.

The post-SN eccentricity is given by:


where the post-SN orbital energy of the system is given by: , and the orbital angular momentum is given by:


where is the reduced mass of the post-SN binary. We neglected any shell impact effects on the companion star since even for SN Ib/c where a significant shell is ejected, the impact effect on the orbital dynamics is small if the pre-SN separation is larger than a few (Liu et al., 2015).

As we have seen already in Fig. 10, a pre-SN orbital period of the tertiary star of  d ( yr) is much too wide to keep the system bound after the SN explosion. Therefore, the value of must have been much smaller and in our chosen example  d. As a result, the eccentricity of the post-SN system must become very large in order for the K-dwarf (star B) to reach post-SN orbital periods of  yr, as derived in Section 4.5.

Figure 12: Post-SN orbital period of the tertiary star, as a function of the eccentricity of its orbit using our selected inner binary model from Fig. 11. The pre-SN model has an orbital period of  d. Solutions are shown along the blue curve for different values of the kick angle, in steps of . The grey lines shown minimum orbital periods for a given eccentricity to ensure a long-term stable three-body configuration. The red line shows the expected orbital period after evaporation of the secondary star – see text for discussions.

In Fig. 12, we plot the post-SN orbital period of the tertiary star, as a function of the eccentricity of its orbit using our selected inner binary model with  km s and , leading to  km s,  d and for the post-SN inner binary, and assuming a pre-SN orbital period of the tertiary star of  d. We only find solutions for a system mimicking PSR J10240719 for a very small interval of given our particular choice of pre-SN parameters. The values of are, as expected, close to the critical value of where the system dissociates (). We can calculate for this system from Eq. (5.2). The red line indicates the final orbital period of the tertiary star as a consequence of orbital widening when the secondary star has been evaporated (assuming here that the NS has accreted 0.1 M during the LMXB phase). The final 3D systemic velocity of our model is  km s, which is sufficient to explain the Galactic motion of PSR J10240719 (Fig. 8). As we discuss below, the post-SN orbit of the tertiary star (star B, the K-dwarf) cannot be too eccentric since this would prevent the post-SN triple system from remaining bound on a long timescale.

5.3 Long-term stability of a triple system

For all our solutions of bound post-SN triples, we must check if their mutual stellar orbits are dynamically stable on a long-term timescale. One must keep in mind that it is expected to last several Gyr from the SN explosion until the secondary star evolves through the LMXB phase and finally becomes evaporated. In all this interval of time, the triple system must be stable. This means that for a given eccentricity of the outer binary, we must require that the ratio of the periastron separation of the tertiary star to the semi-major axis of the inner binary, () be larger than a critical value. Otherwise, three-body interactions will cause the hierarchical triple to be unstable and eject of one of the stars via tidal disruption or chaotic energy exchange (Mardling & Aarseth, 2001). This criterion translates into the requirement that the post-SN orbital period of the tertiary star, , with a given eccentricity, must be larger than a critical value.

A number of stability criteria for triple systems have been proposed over the last four decades (see Mikkola, 2008, for an overview). The grey lines in Fig. 12 show stability criteria for triple systems as suggested by Mardling & Aarseth (2001); Eggleton & Kiseleva (1995); Bailyn (1987); Harrington (1972), indicated by MA, EK, B and H, respectively. The earlier of these studies have been derived empirically for restricted configurations. The most stringent of the stability criteria is by far that of Mardling & Aarseth (2001, Eq. 90);


where , and which can be combined with Kepler’s third law to yield critical orbital periods as plotted in Fig. 12. We note that for our chosen example, the triple system would become unbound for  yr. Nevertheless, for  yr (corresponding to a present binary orbital period of  yr after evaporation of the secondary star) we can reproduce PSR J10240719 for the triple model investigated here.

Finally, it should be noted that the strict stability criteria of Mardling & Aarseth (2001) was derived for coplanar prograde motion in accordance with a recoil velocity of the inner binary in the orbital plane of the outer orbit. Hence, this criteria represents an upper limit for stability of non-coplanar systems since triple systems whose inner and outer orbits are inclined, despite Kozai interactions (Kozai, 1962), will be more stable than coplanar prograde systems with the same mass ratios and eccentricities. Mardling & Aarseth (2001) estimate that for non-coplanar systems, the ratio () could be around 30 per cent smaller and still result in long-term stability against escape in the three-body system. This would slightly enhance the probability for producing a system like PSR J10240719.

We have demonstrated that PSR J10240719 could have formed via a triple system in the Galactic disk. However, estimating the probability for this formation channel would require an investigation combining population synthesis, gravitational dynamics, detailed stellar evolution and hydrodynamical interactions of the common envelope phase, which is much beyond the scope of this paper. While it is difficult to estimate the rate at which binaries like PSR J10240719 could be produced from a triple star scenario, we have demonstrated that this scenario is possible, although severe fine-tuning is required for this to happen. Nevertheless, so far we have only detected one such system in our Galaxy.

5.4 Alternative models

A different formation scenario has been suggested by Kaplan et al. (2016). In the scenario, PSR J10240719 has been recycled in the center of a globular cluster, where an exchange encounter with another cluster star or binary led to the ejection of the pulsar from the core of the cluster. In this scenario, the major part of the high space velocity is simply the original velocity of the globular cluster.

The ejection of the inner companion due to a dynamical instability in a triple system, studied in detail by Portegies Zwart et al. (2011) for PSR J1903+0327, does not seem possible for PSR J10240719. The reason is that the orbit of the outer companion (the K-dwarf/star B) is expected to harden (i.e. become tighter) if the inner companion is ejected. However, given the current extremely wide orbit, any release of orbital binding energy would have been insufficient to energetically eject the inner companion. The possibility that the ejection event of the inner companion also caused the outer companion to be perturbed into its current extremely wide orbit, from an initially much closer orbit, is not allowed by energy considerations.

Detailed simulations of either the dynamically stable triple evolution scenario presented here, or the globular cluster ejection scenario suggested by Kaplan et al. (2016) will be required to ascertain which scenario is most probable. When better constraints on the orbit of PSR J10240719 become available in the future, predictions from these scenarios can be tested.

6 Discussion and conclusions

Since its discovery in 1994, PSR J10240719 was thought to be an isolated millisecond pulsar. Our extended EPTA timing ephemeris of the pulsar, combined with optical astrometry, photometry and spectroscopy of the 2MASS J102438690719190 (star B), show that PSR J10240719 is gravitationally bound to star B, and that the gravitational interaction between the two can account for the apparent discrepancy in the system distance and the steep spectrum timing noise of PSR J10240719. We furthermore find that, even though we cannot fully constrain the binary orbit, it must be extremely wide, with orbital periods in excess of 200 yrs. A rare formation scenario of PSR J10240719 is required to explain the extremely wide binary orbit and the very high space velocity in the Galaxy. We show that PSR J10240719 could have evolved from a triple system, though significant fine-tuning during the supernova explosion is required.

The binary orbit of PSR J10240719 and star B is presently not fully constrained. As a result, the higher order frequency derivatives in the radio timing solution of PSR J10240719 will remain free parameters and will hence be a source of timing noise which may limit using PSR J10240719 as a pulsar timing array source. To fully constrain the orbit by timing PSR J10240719 alone, a measurement of five spin frequency derivatives would be required (Joshi & Rasio, 1997). Considering the component masses as known, and assuming a value for either or , the frequency derivatives would constrain , , , and and provide a description of the timing using the spin frequency and spin frequency derivative as with other pulsars.

Improvements in the timing accuracy of PSR J10240719 are expected in future IPTA data releases which combine EPTA, PPTA and NANOGrav observations. The first IPTA data release (Verbiest et al., 2016) contains 15.9 yrs of EPTA and PPTA data on PSR J10240719. That release contains a subset of the EPTA data used in this paper and hence does not provide stronger constraints than those in Table 1. Combining the timing data presented here with that of NANOGrav (Matthews et al., 2016; Kaplan et al., 2016) and the PPTA (Reardon et al., 2016) should provide an increase in the timing accuracy which, to first order, scales with the square root of the number of TOAs. The uncertainties on the spin frequency derivatives with decrease over time with , where for , for , for and so on. Hence, the uncertainties on will already be halved by extending the timing baseline another 4 years.

Independent astrometry of PSR J10240719 will be provided in the near future through an ongoing VLBA pulsar astrometry project (Deller et al., 2011). The expected uncertainties on position and proper motions will be comparable to that from pulsar timing, while the uncertainty on the pulsar parallax will be significantly better (factor 5 to 10) compared to the current uncertainty of 0.11 mas on the timing parallax (Deller, priv. comm.).

An improvement in the position and proper motion of star B would also help constrain the orbit through Eqns. 6 and 7. Though this would require including the longitude of ascending node as a free parameter, a measurement of differential position and proper motion as well as two spin frequency derivatives would allow a complete description of the orbit. Improving the astrometry presented here through further ground-based observations will be challenging. However, star B is bright enough for its parallax, position and proper motion to be measured by Gaia. The photometric transformations by Jordi et al. (2010) estimate a Gaia -band magnitude of for its -band magnitude and color. The astrometric performance of Gaia555http://www.cosmos.esa.int/web/gaia/science-performance predicts uncertainties on the parallax, position and proper motion of 0.32 mas, 0.23 mas and 0.17 mas yr, respectively. These uncertainties are comparable to those of PSR J10240719 from the timing solution.

Combining the Gaia astrometry of star B with the pulsar timing astrometry of PSR J10240719 will provide much stronger constraints on the orbit. For the nominal distance of  kpc, these uncertainties translate to a distance uncertainty of 500 pc, an uncertainty in the projected position on the sky of 0.3 AU and a projected velocity uncertainty of 0.9 km s. Based on the orbital solutions shown in Fig. 7, the projected velocities on the sky at the nominal distance are in the range of 0.3 to a few km s, where the velocity is larger for less eccentric orbits and with the binary members away from apastron, and hence may be measurable from the astrometry of both binary members. This approach will also provide independent constraints on the distance and the intrinsic spin period derivative of PSR J10240719.

The full Gaia catalog is expected to be released in 2022. Given the prospect of much tighter constraints on the orbital parameters with this future release, and considering the large pulsar timing data set that has already been accumulated for PSR J10240719, we consider it worthwhile to keep observing the pulsar as part of pulsar timing arrays. Besides contributing to the constraints on nano-hertz gravitational waves, these decade long, multi-frequency, pulsar timing data sets like that presented here enable us to measure other effects perturbing the stable rotation of millisecond pulsars and make unexpected discoveries.


We thank David Kaplan for fruitful discussions about the system and the orbital constraints, which helped resolve a mistake in the derivation of the spin frequency derivatives due to orbital motion. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. This research is based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 95.D-0973(A). Part of this work is based on observations obtained the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. Pulsar research at Jodrell Bank and access to the Lovell Telescope is supported by a Consolidated Grant from the UK’s Science and Technology Facilities Council. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). The Westerbork Synthesis Radio Telescope is operated by the ASTRON (Netherlands Institute for Radio Astronomy) with support from the Netherlands Foundation for Scientific Research (NWO). The authors acknowledge the support of the collegues in the European Pulsar Timing Array. CGB acknowledges support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement nr. 337062 (DRAGNET; PI Jason Hessels). SO is supported by the Alexander von Humboldt Foundation. KL acknowledges the financial support by the European Research Council for the ERC Synergy Grant BlackHoleCam under contract no. 610058. RNC and PL acknowledge the support of the International Max Planck Research School Bonn/Cologne and the Bonn-Cologne Graduate School. GT and IC acknowledge financial support from ‘Programme National de Cosmologie and Galaxies’ (PNCG) of CNRS/INSU, France.


  • Abbott et al. (2016) Abbott B. P., et al., 2016, Physical Review Letters, 116, 061102
  • Alpar et al. (1982) Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982, Nature, 300, 728
  • Antoniadis et al. (2012) Antoniadis J., van Kerkwijk M. H., Koester D., Freire P. C. C., Wex N., Tauris T. M., Kramer M., Bassa C. G., 2012, MNRAS, 423, 3316
  • Appenzeller et al. (1998) Appenzeller I., et al., 1998, The Messenger, 94, 1
  • Arzoumanian et al. (2015) Arzoumanian Z., et al., 2015, ApJ, 813, 65
  • Bailes et al. (1997) Bailes M., et al., 1997, ApJ, 481, 386
  • Bailyn (1987) Bailyn C. D., 1987, PhD thesis, Harvard University, Cambridge, MA.
  • Bassa et al. (2006) Bassa C. G., van Kerkwijk M. H., Koester D., Verbunt F., 2006, A&A, 456, 295
  • Bassa et al. (2016) Bassa C. G., et al., 2016, MNRAS, 456, 2196
  • Bessell et al. (1998) Bessell M. S., Castelli F., Plez B., 1998, A&A, 333, 231
  • Bhattacharya & van den Heuvel (1991) Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203, 1
  • Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn J., Gerhard O., 2016, preprint, (arXiv:1602.07702)
  • Bovy (2015) Bovy J., 2015, ApJS, 216, 29
  • Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  • Cardelli et al. (1989) Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
  • Champion et al. (2008) Champion D. J., et al., 2008, Science, 320, 1309
  • Chen et al. (2013) Chen H.-L., Chen X., Tauris T. M., Han Z., 2013, ApJ, 775, 27
  • Chen et al. (2014) Chen Y., Girardi L., Bressan A., Marigo P., Barbieri M., Kong X., 2014, MNRAS, 444, 2525
  • Chen et al. (2015) Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068
  • Cognard & Theureau (2006) Cognard I., Theureau G., 2006, in IAU Joint Discussion.
  • Cognard et al. (2013) Cognard I., Theureau G., Guillemot L., Liu K., Lassus A., Desvignes G., 2013, in Cambresy L., Martins F., Nuss E., Palacios A., eds, SF2A-2013: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. pp 327–330
  • Cutri et al. (2003) Cutri R. M., et al., 2003, VizieR Online Data Catalog, 2246
  • Damour & Taylor (1992) Damour T., Taylor J. H., 1992, Phys. Rev. D, 45, 1840
  • Deller et al. (2011) Deller A. T., et al., 2011, in Alef W., Bernhart S., Nothnagel A., eds, 20th Meeting of the European VLBI Group for Geodesy and Astronomy, held in Bonn, Germany, March 29-30, 2011, Eds: W. Alef, S. Bernhart, and A. Nothnagel, Institut für Geodäsie und Geoinformation, Rheinischen Friedrich-Wilhelms-Universität Bonn, p. 178-182. pp 178–182 (arXiv:1110.1979)
  • Demorest et al. (2013) Demorest P. B., et al., 2013, ApJ, 762, 94
  • Desvignes et al. (2016) Desvignes G., et al., 2016, preprint, (arXiv:1602.08511)
  • Detweiler (1979) Detweiler S., 1979, ApJ, 234, 1100
  • Edwards et al. (2006) Edwards R. T., Hobbs G. B., Manchester R. N., 2006, MNRAS, 372, 1549
  • Eggleton & Kiseleva (1995) Eggleton P., Kiseleva L., 1995, ApJ, 455, 640
  • Espinoza et al. (2013) Espinoza C. M., et al., 2013, MNRAS, 430, 571
  • Flannery & van den Heuvel (1975) Flannery B. P., van den Heuvel E. P. J., 1975, A&A, 39, 61
  • Folkner et al. (2009) Folkner W. M., Williams J. G., Boggs D. H., 2009, Interplanetary Network Progress Report, 178, 1
  • Freire et al. (2001) Freire P. C., Kramer M., Lyne A. G., 2001, MNRAS, 322, 885
  • Freire et al. (2011) Freire P. C. C., et al., 2011, MNRAS, 412, 2763
  • Green et al. (2014) Green G. M., et al., 2014, ApJ, 783, 114
  • Green et al. (2015) Green G. M., et al., 2015, ApJ, 810, 25
  • Guillemot et al. (2016) Guillemot L., et al., 2016, A&A, 587, A109
  • Hamuy et al. (1992) Hamuy M., Walker A. R., Suntzeff N. B., Gigoux P., Heathcote S. R., Phillips M. M., 1992, PASP, 104, 533
  • Hamuy et al. (1994) Hamuy M., Suntzeff N. B., Heathcote S. R., Walker A. R., Gigoux P., Phillips M. M., 1994, PASP, 106, 566
  • Harrington (1972) Harrington R. S., 1972, Celestial Mechanics, 6, 322
  • Hellings & Downs (1983) Hellings R. W., Downs G. S., 1983, ApJ, 265, L39
  • Hills (1983) Hills J. G., 1983, ApJ, 267, 322
  • Hobbs et al. (2006) Hobbs G. B., Edwards R. T., Manchester R. N., 2006, MNRAS, 369, 655
  • Hogg et al. (2005) Hogg D. W., Blanton M. R., Roweis S. T., Johnston K. V., 2005, ApJ, 629, 268
  • Horne (1986) Horne K., 1986, PASP, 98, 609
  • Istrate et al. (2014) Istrate A. G., Tauris T. M., Langer N., 2014, A&A, 571, A45
  • Johnson & Soderblom (1987) Johnson D. R. H., Soderblom D. R., 1987, AJ, 93, 864
  • Jordi et al. (2010) Jordi C., et al., 2010, A&A, 523, A48
  • Joshi & Rasio (1997) Joshi K. J., Rasio F. A., 1997, ApJ, 479, 948
  • Kaplan et al. (2013) Kaplan D. L., Bhalerao V. B., van Kerkwijk M. H., Koester D., Kulkarni S. R., Stovall K., 2013, ApJ, 765, 158
  • Kaplan et al. (2016) Kaplan D., Kupfer T., Nice D., Irrgang A., Heber U., 2016, in prep.
  • Karuppusamy et al. (2008) Karuppusamy R., Stappers B., van Straten W., 2008, PASP, 120, 191
  • Khargharia et al. (2012) Khargharia J., Stocke J. T., Froning C. S., Gopakumar A., Joshi B. C., 2012, ApJ, 744, 183
  • Kopeikin (1996) Kopeikin S. M., 1996, ApJ, 467, L93
  • Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  • Kuijken (2011) Kuijken K., 2011, The Messenger, 146, 8
  • Landolt (1992) Landolt A. U., 1992, AJ, 104, 340
  • Lasker et al. (2008) Lasker B. M., et al., 2008, AJ, 136, 735
  • Lazarus et al. (2016) Lazarus P., Karuppusamy R., Graikou E., Caballero R. N., Champion D. J., Lee K. J., Verbiest J. P. W., Kramer M., 2016, preprint, (arXiv:1601.06194)
  • Le Borgne et al. (2003) Le Borgne J.-F., et al., 2003, A&A, 402, 433
  • Lentati et al. (2014) Lentati L., Alexander P., Hobson M. P., Feroz F., van Haasteren R., Lee K. J., Shannon R. M., 2014, Mon. Not. R. Astron. Soc., 437, 3004
  • Lentati et al. (2015) Lentati L., et al., 2015, MNRAS, 453, 2576
  • Lentati et al. (2016) Lentati L., et al., 2016, MNRAS, 458, 2161
  • Liu et al. (2015) Liu Z.-W., Tauris T. M., Röpke F. K., Moriya T. J., Kruckow M., Stancliffe R. J., Izzard R. G., 2015, A&A, 584, A11
  • Lutz & Kelker (1973) Lutz T. E., Kelker D. H., 1973, PASP, 85, 573
  • Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  • Manchester et al. (2013) Manchester R. N., et al., 2013, Publ. Astron. Soc. Australia, 30, 17
  • Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  • Matthews et al. (2016) Matthews A. M., et al., 2016, ApJ, 818, 92
  • Mikkola (2008) Mikkola S., 2008, in Hubrig S., Petr-Gotzens M., Tokovinin A., eds, Multiple Stars Across the H-R Diagram. p. 11, doi:10.1007/978-3-540-74745-1_2
  • Monet et al. (1998) Monet D., Bird A., Canzian B., et al. 1998, The PMM USNO-A2.0 Catalog. US Naval Observatory
  • Monet et al. (2003) Monet D. G., et al., 2003, AJ, 125, 984
  • Munari et al. (2005) Munari U., Sordo R., Castelli F., Zwitter T., 2005, A&A, 442, 1127
  • Nice & Taylor (1995) Nice D. J., Taylor J. H., 1995, ApJ, 441, 429
  • Petit (2010) Petit G., 2010, Highlights of Astronomy, 15, 220
  • Pijloo et al. (2012) Pijloo J. T., Caputo D. P., Portegies Zwart S. F., 2012, MNRAS, 424, 2914
  • Podsiadlowski et al. (2002) Podsiadlowski P., Rappaport S., Pfahl E. D., 2002, ApJ, 565, 1107
  • Portegies Zwart et al. (2011) Portegies Zwart S., van den Heuvel E. P. J., van Leeuwen J., Nelemans G., 2011, ApJ, 734, 55
  • Qi et al. (2015) Qi Z., et al., 2015, AJ, 150, 137
  • Reardon et al. (2016) Reardon D. J., et al., 2016, MNRAS, 455, 1751
  • Roeser et al. (2010) Roeser S., Demleitner M., Schilbach E., 2010, AJ, 139, 2440
  • Schlafly & Finkbeiner (2011) Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103
  • Shanks et al. (2015) Shanks T., et al., 2015, MNRAS, 451, 4238
  • Shannon et al. (2013) Shannon R. M., et al., 2013, Science, 342, 334
  • Shemar & Lyne (1996) Shemar S. L., Lyne A. G., 1996, MNRAS, 282, 677
  • Shklovskii (1970) Shklovskii I. S., 1970, Soviet Astronomy, 13, 562
  • Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  • Stetson (1987) Stetson P. B., 1987, PASP, 99, 191
  • Stetson (2000) Stetson P. B., 2000, PASP, 112, 925
  • Sutaria et al. (2003) Sutaria F. K., Ray A., Reisenegger A., Hertling G., Quintana H., Minniti D., 2003, A&A, 406, 245
  • Tang et al. (2014) Tang J., Bressan A., Rosenfield P., Slemer A., Marigo P., Girardi L., Bianchi L., 2014, MNRAS, 445, 4287
  • Tauris (2011) Tauris T. M., 2011, in Schmidtobreick L., Schreiber M. R., Tappert C., eds, Astronomical Society of the Pacific Conference Series Vol. 447, Evolution of Compact Binaries. p. 285 (arXiv:1106.0897)
  • Tauris & Takens (1998) Tauris T. M., Takens R. J., 1998, A&A, 330, 1047
  • Tauris & van den Heuvel (2014) Tauris T. M., van den Heuvel E. P. J., 2014, ApJ, 781, L13
  • Tauris et al. (2012) Tauris T. M., Langer N., Kramer M., 2012, MNRAS, 425, 1601
  • Verbiest et al. (2009) Verbiest J. P. W., et al., 2009, MNRAS, 400, 951
  • Verbiest et al. (2010) Verbiest J. P. W., Lorimer D. R., McLaughlin M. A., 2010, MNRAS, 405, 564
  • Verbiest et al. (2016) Verbiest J. P. W., et al., 2016, MNRAS, 458, 1267
  • Zacharias et al. (2013) Zacharias N., Finch C. T., Girard T. M., Henden A., Bartlett J. L., Monet D. G., Zacharias M. I., 2013, AJ, 145, 44
  • van Kerkwijk et al. (2005) van Kerkwijk M. H., Bassa C. G., Jacoby B. A., Jonker P. G., 2005, in Rasio F. A., Stairs I. H., eds, ASP Conference Series Vol. 328, Binary Radio Pulsars. ASP, p. 357

Appendix A Spin frequency derivatives due to orbital motion

We start by expressing the orbit of a binary system on the sky in terms of a projected separation and a position angle :


Here, is the true anomaly, the argument of perigee, the longitude of the ascending node and the inclination. is the location along the line-of-sight, while and represent the position on the sky, but without correcting for . We use the pulsar timing definition by Damour & Taylor (1992) and Kopeikin (1996) which has positive away from the observer. As a result, a right-handed system would have increase in a clockwise direction. This is opposite to the more standard definition of measuring position angles anti-clockwise on the sky, i.e. from North to East.

We follow Freire et al. (2001) to obtain the velocities in the , and axes, which are defined as:


Here, is the semi-major axis and we have used that


Differentiating further for the radial position gives the acceleration , jerk , snap and crackle :