The Binary Nature of PSR J2032+4127
PSR J2032+4127 is a -ray and radio-emitting pulsar which has been regarded as a young luminous isolated neutron star. However, its recent spin-down rate has extraordinarily increased by a factor of two. We present evidence that this is due to its motion as a member of a highly-eccentric binary system with a -M Be star, MT91 213. Timing observations show that, not only are the positions of the two stars coincident within 0.4, but timing models of binary motion of the pulsar fit the data much better than a model of a young isolated pulsar. MT91 213, and hence the pulsar, lie in the Cyg OB2 stellar association, which is at a distance of only 1.4–1.7 kpc. The pulsar is currently on the near side of, and accelerating towards, the Be star, with an orbital period of 20–30 years. The next periastron is well-constrained to occur in early 2018, providing an opportunity to observe enhanced high-energy emission as seen in other Be-star binary systems.
keywords:stars: neutron – pulsars: general
PSR J2032+4127 is a 143-ms pulsar discovered by the Large Area Telescope (LAT) of the Fermi Gamma-ray Space Telescope (Abdo et al., 2009) and subsequently detected at radio frequencies (Camilo et al., 2009). Using -ray data as well as radio data, pulse arrival-time analysis by those authors, and latterly by Ray et al. (2011), showed that the pulsar had a large slow-down rate, indicating that it was a young pulsar, a notion which has been supported by a recent glitch in its rotation rate. The timing analysis showed that the projected position of the pulsar lay close to a V=11.95 Be star, MT91 213 in the Cyg OB2 stellar association. However, Camilo et al. (2009) concluded that it was not a binary companion of PSR J2032+4127, based mainly upon the apparent lack of variations in the pulsar rotation rate that would arise from the Doppler effects of any reasonable (circular) orbital binary motion. Subsequent studies of the properties of the pulsar have all been based upon the assumption that it is a solitary young energetic pulsar. For instance, Aliu et al. (2014) present a VERITAS detection of TeV J2032+4130 and they discuss its likely association with PSR J2032+4127, noting that the extended nature of the source and the prevalence of pulsar wind nebulae (PWNe) in the Galactic TeV source population argues that it is a PWN powered by the pulsar. However they note that the extended X-ray emission that is also spatially coincident with the source (Horns et al., 2007; Butt et al., 2006; Murakami et al., 2011) is quite weak if it is from a PWN. Six years of both -ray and radio timing data are now available, and we revisit the possibility that the pulsar is in orbit with MT91 213.
2 Observations and Timing Noise Model
Fermi LAT data were used from shortly after the start of the mission in August 2008 up to June 2014 (MJD 54682–56824). Following the maximum likelihood method of Ray et al. (2011), we constructed times-of-arrival (TOAs) with a 14-day cadence from “reprocessed” Pass 7 Fermi LAT data. For this analysis, we used SOURCE class photons, excluding events with a zenith angle and when the spacecraft rocking angle exceeds . To improve sensitivity, we employed photon weighting using the spectral models available in The Second Fermi Large Area Telescope Catalog of Gamma-ray Pulsars (Abdo et al., 2013).
Radio timing observations were made with the NRAO Green Bank Telescope (GBT) and the Lovell Telescope (LT) at Jodrell Bank. The GBT observations were primarily in bands centred on 820 MHz and 2000 MHz during the early part of this period (MJD 54836–55589) (Camilo et al., 2009), as well as MJD 56855–56857. Radio observations with the LT were made at approximately weekly intervals at around 1520 MHz for MJD 55222–56830. As described by Camilo et al. (2009), the -ray and radio pulse profiles are very different in form. In order to establish the true alignment of the profiles, the dispersion measure (DM) was initially determined using only radio data. Inspection of the radio pulse profiles shows that there is little change in shape between 800 MHz and 2000 MHz, and a single standard profile was used to obtain TOAs for all the radio observations. The DM was determined at three epochs, using the three 800-MHz TOAs and their neighbouring 1520-MHz and 2000-MHz TOAs. This allowed the time offset of the contemporaneous -ray TOAs to be established using the jump facility in tempo2 (Hobbs et al., 2006). The mean value of these time offsets was consistent with the alignment obtained by Camilo et al. (2009) and was applied to all the -ray TOAs, allowing the delay between these TOAs and the radio TOAs to be used to determine the DM throughout the data set.
The radio TOAs have errors of s, about half the errors of those provided by the Fermi LAT. In total, about 400 good TOAs are available over six years (MJD 54682–56857).
The pulsar shows substantial variations in rotation rate, indicating a large amount of “timing noise”, which is not unexpected in a young pulsar (i.e. having a large intrinsic spin-down rate). However, it has recently become clear that the magnitude of these variations far exceeds any previously encountered timing noise in a pulsar (e.g. Hobbs et al. 2010). Fig. 1a presents the evolution of the barycentric rotation frequency of the pulsar over the 6-year span of the timing data, showing a substantial reduction in frequency during this time. Also apparent is an increase in the magnitude of the slope, corresponding to an increase in the magnitude of the frequency derivative, shown in Fig. 1b, by more than a factor of two, and a decrease in the characteristic age by the same factor, confounding any interpretation in terms of conventional pulsar spin-down. There is also a barely-discernable minor glitch which occurred in Sept 2011 (MJD 55811) and which can be fitted by steps in frequency () and frequency derivative (). The whole set of TOAs can be well-fitted by a model involving the position, three glitch parameters, the rotation frequency and 7 derivatives (Table 1, column 2). The timing residuals (the difference between the observed and model-based TOAs) relative to this model, which we will henceforth refer to as the “noise model”, are shown in Fig. 1c.
The behaviour in the rotation seen in Figs. 1a and 1b has never been seen in any other isolated pulsar, and is not characteristic of either a post-glitch recovery, during which the magnitude of the slow-down rate is usually observed to decrease (e.g. Espinoza et al. 2011), or of timing noise, which is usually characterised by switches, often quasi-periodic, between values of slow-down rate rather than the smooth variation seen here (Lyne et al., 2010; Lyne, 2013). Most anomalous of all is the doubling of the slow-down rate, also not seen in any other isolated pulsar.
|Right Ascension, (J2000.0)|
|Epoch of frequency, (MJD)||55714.0||55700.0||55700.0|
|Freq 1st deriv, (s)||1.16583(12)||0.61(5)||0.61(5)|
|Glitch Epoch, (MJD)||55810.51(2)||55810.76(7)||55810.76(7)|
|Freq, ( Hz)||1.907(2)||1.907(2)||1.907(2)|
|Freq 1st deriv, (s)||4(3)||10(2)||10(2)|
|DM (pc cm)||114.66(3)||114.65(3)||114.65(3)|
|DM deriv, DM1 (pc cmy)||0.00(1)||0.00(1)||0.00(1)|
|Timing noise parameters:||Binary parameters:|
|Freq 2nd deriv, (s)||0.3651(7)||Orbital period, (d)||8578(1200)||8625(1200)|
|Freq 3rd deriv, (s)||0.439(2)||Epoch of Periastron, (MJD)||58161(31)||58153(31)|
|Freq 4th deriv, (s)||0.690(14)||Proj. semi-major axis, (lt-s)||8819(1250)||11151(1250)|
|Freq 5th deriv, (s)||1.26(9)||Eccentricity,||0.93(3)||0.95(3)|
|Freq 6th deriv, (s)||6.2(3)||Long. of Periastron (deg)||12(4)||10(4)|
|Freq 7th deriv, (s)||27(3)|
|Rms timing residual, (ms)||0.53||0.57||0.56|
3 A binary Model
We believe that the only plausible origin of such a large variation in the observed slow-down rate of a pulsar lies in the Doppler effects of binary motion with another star. While we note the remark by Camilo et al. (2009) that there is no evidence of any short-period binary motion and that the binary period (of any circular orbit) must be in excess of 100 years, we have explored the possibility that the pulsar is actually a member of a long-period binary system with a large orbital eccentricity.
First, we sought fits of binary models to the rotation-frequency data of Fig. 1a. While good fits to the data were possible using models for eccentric binary orbits with orbital periods in excess of about 6000 days (16 years), it soon became clear that there were strong covariances between some of the fitted parameters, arising from the small orbital phase range of the available data. In particular, there were large covariances between the intrinsic pulsar slow-down rate, , the orbital period and the projected semi-major axis of the orbit , where is the semi-major axis of the orbit, is the inclination of the plane of the orbit to the plane of the sky and is the speed of light. These covariances allowed many good, but not unique, fits to the data. We therefore explored fits to the data for a range of fixed values of between 6000 and 12000 days at 500-day intervals, and, for each of these, a range of fixed values of , corresponding to fixed values of mass function of 2, 5, 10 and 20 . is a function of the masses of the neutron star () and its companion () and orbital inclination and is determined from and from Kepler’s laws by:
where is Newton’s gravitational constant. If and are in solar masses, is in light-sec and in days, then
Fig. 2 shows an example of one of the best fits to the data, for = 8578 days, = 10 M. Note that the data span occupies only about 20% of the orbit. Remarkably, the root mean square (rms) of the frequency residuals (the differences between the measured and model frequency values) for this simple model was only approximately of the total frequency variation during this time and is consistent with the measurement errors.
For each , pair, the ephemeris resulting from the frequency fit was subsequently used as the basis of a coherent timing analysis of the TOAs using tempo2. Keeping and at their fixed values, the three other Keplerian parameters, the pulsar position, DM and its first derivative, DM1, the rotation frequency and its first derivative and three glitch parameters were all fitted to the TOAs. Fig. 3 summarises the results of the fits. In particular, Fig. 3a shows the rms of the timing residuals relative to the fitted model as a function of and . There is a broad minimum in the rms with a value of about 0.6 ms for 7500 d 9500 d and all values of provide indistinguishably good fits. Best-fit models were obtained for and by including a fit for and the results are given in Table 1, columns 4 and 5. Because of the strong covariance between some of the parameters, we used a Monte-Carlo method to estimate the errors in the fitted parameters that are given in the table. The Monte-Carlo analysis used simulations of the pulsar observations seeded with the best-fit parameter file. The simulated data sets had the same epochs of the combined LAT and radio TOAs, giving the same cadence and orbital coverage as the real data. A noise source was added to the data, defined by the power spectral density:
where the white noise, yr, and the power-law parameters are yr, yr and . The noise parameters were chosen to match the observed noise spectrum. The toasim plugins from tempo2 were used to generate realisations of the noise and the full fitting process was applied to each realisation. The error estimates were determined from the variance of the resultant fit parameters.
Fig. 1e shows the timing residuals for the binary fit for given in Table 1, column 4, showing that the rms is almost entirely limited by measurement errors, with any remaining unmodelled systematic trends at a level that is comparable with those from pulsars of a similar age. We note that these binary models involve the fitting of 14 free parameters to achieve an rms residual of 0.57 ms. This is one parameter fewer than required for the timing-noise model with similar rms residuals that is presented in Fig. 1c and Table 1, and based upon a 7-derivative polynomial model, achieving an rms residual of 0.53 ms. Fig. 1d shows the residuals based upon a 6-derivative timing-noise model that has the same number of fitted parameters as the binary model and is clearly a poorer description of the data.
We conclude that there is a range of simple binary models that describe the TOA data very well. It is also clear that the reason we cannot define the orbit more precisely at present is that, because the data span covers only about 1/5 of an orbit, the differences between the models are buried within the measurement errors of the TOAs and likely timing noise.
However, Figs. 3b–3f show that the values of many of the fitted parameters do not change much near the rms residual minimum. For all the models near this minimum (Fig. 3a), the eccentricity increases from about 0.90 to 0.95 with increasing (Fig. 3c). However, the corresponding lies in a restricted range between MJD 58150 and 58200 (Figs. 3b and 4). Thus the next periastron passage will occur sometime during 2018 February or March. The intrinsic frequency derivative for all these models is about s (Fig. 3d), corresponding to a characteristic age of about 180 kyr, which is substantially greater than the (rapidly decreasing!) characteristic age of the pulsar in the young-isolated-pulsar model, which has changed from kyr to kyr over the six years. The values of right ascension and declination are not strongly model-dependent (Figs. 3e and 3f) and are about and , consistent with the position of the Be star, MT91 213, which is at and 111Owing to a typographical error, Camilo et al. (2009) listed the last digits of declination of MT91 213 as 24.28; the correct value is 24.48 as used here., and which Camilo et al. (2009) considered, but rejected, as a possible companion of the pulsar. The conflicting conclusion of that study with this paper arose from their assumption of a circular orbit and the happenstance that their observations were conducted near apastron, where the gravitational influence of the companion on the pulsar is small, so that Doppler effects result in a small value of the pulsar second rotational frequency derivative.
If MT91 213 is indeed the companion star, what additional constraints might be placed upon the system? Kiminki et al. (2007) estimate the mass of the star to be M, based upon their estimate of its spectral type (B0V) and the study by Martins et al. (2005) who suggest that such spectrally-determined mass estimates may be in error by 35–50%. However, more recently Hohle et al. (2010) have estimated the mass of B0V stars to be M using 2MASS spectrometry and Hipparchos parallax data. We adopt this latter value in the following discussion. Unfortunately, the inclination of the plane of the orbit to the plane of the sky is unknown. For random orientations of a binary orbit, the median of the probability density function of inclination occurs at . For an orbit with this inclination and a pulsar of mass , the mass function (Equation 1) would thus be , which corresponds approximately to the first binary model in Table 1, column 4, and the orange lines in Fig. 3. The greatest possible value of the mass function for these two stars is about 15 M which occurs if . There is no doubt that an orbit with a mass function of M, having an orbital period lying between 7000 and 11000 days, or about 20–30 years, is entirely consistent with the data (Fig. 4). Although larger mass functions also provide satisfactory fits to the data, they are inconsistent with a companion that has the mass of MT91 213.
Two other pulsar/Be-star binary systems are known, PSR B125963/SS 2883 (Johnston et al., 1992) and PSR J16384725 (Lorimer et al., 2006; Lyne, 2008), having orbital periods of 1237 d and 1941 d respectively, and companion masses of and . The PSR J2032+4127/MT91 213 system is more extreme than either of these, with d and M. Fig. 5 shows the most likely orbital configuration. The pulsar is presently on the Earth side of the Be star and is beginning to move rapidly in towards periastron, after which it will move behind the star and any circumstellar disk. Such stars have substantial stellar winds and, in the cases of both PSR B125963 and PSR J16384725, as well as radio eclipses, there are increases in both dispersion measure (DM) and multi-path scattering close to periastron, which arise as the pulsar becomes more embedded in the dense circumstellar environment. No significant change has yet been detected in the scattering or DM of PSR J2032+4127, with the fitted values of its first derivative DM1 consistent with zero (Table 1), but it is still far from periastron and this is not unexpected. Following periastron in early 2018, the pulsar will move behind MT91 213 and may be eclipsed, at least in the radio, by its atmosphere or disk, depending upon the inclination of the orbit to the line-of-sight.
Such a large, long-period system is unique, but its existence and its survival of the velocity kick that the pulsar probably experienced when it was formed is perhaps not surprising, because of the deep gravitational well of the -M companion. The separation of the markers in Fig. 5 gives an indication of the orbital velocity, showing that, although the velocity at apastron is only a few km s, at periastron the velocity is in excess of 100 km s. It is a system that, unlike the majority of pulsars, did not quite become unbound after the neutron star formation.
In all the allowed binary models, periastron will occur in just over three years time. Much of the degeneracy present in the fits will be resolved as periastron is approached. In particular, the rapid increase in slow-down rate during the next 1–2 years will establish the magnitude of the eccentricity and hence the stellar separation at periastron. All evidence indicates that MT91 213 is in the Cyg OB2 stellar association, which is located at a distance of 1.4–1.7 kpc, determined both by spectroscopic parallax and by VLBI determination of the trigonometric parallax of masers associated with the association stars (Massey & Thompson, 1991; Hanson, 2003; Rygl et al., 2012); in the binary hypothesis, this is therefore also the distance of PSR J2032+4127. We note that, at this distance, the projected major axis of the pulsar binary orbit will subtend about 25 mas viewed from the Earth. It should be possible to track the pulsar around the orbit using a VLBI network such as the VLBA or the EVN, hence establishing the size and the projected eccentricity of the orbit. Together with timing measurements, such observations permit the direct determination of the inclination of the orbit to the plane of the sky and the distance of the system from the Earth.
The spin-down luminosity of PSR J2032+4127 in the binary models presented here is erg s, 62% of the Camilo et al. (2009) isolated-pulsar value. The efficiency in converting rotational kinetic energy to -ray luminosity , , is therefore increased by a factor of 1.6. Using the distance of 1.4–1.7 kpc and the updated from Abdo et al. (2013), and assuming as usual an isotropic -ray beam, then 15%–20%, which is unremarkable for a young pulsar (the much higher value of presented in Abdo et al. (2013) uses the DM-based distance of 3.7 kpc).
We also note that this system is still fairly young, the pulsar having a characteristic age of 180 kyr and the Be star having a total lifetime of only a few million years, which needs to encompass the lifetimes of the pulsar progenitor star and of the pulsar itself.
The fact that PSR J2032+4127 is a member of a long period binary system with a massive companion star means that we should revisit the interpretation of the spatially coincident extended TeV -ray and X-ray emission (Aliu et al., 2014). Kargaltsev et al. (2014) have recently shown that the similar system PSR B125963 exhibits variable extended X-ray emission that is located outside the binary system and is moving away. It may be that this is emission that is associated with a circum-binary shock which can be generated between the winds of a pulsar and a massive star (Bosch-Ramon et al., 2012). We note also that if the interpretation of the TeV emission extension is to do with the proper motion of the binary then the increased age will result in a reduced transverse velocity of less than 30 km s compared to that derived by Aliu et al. (2014).
In conclusion, there is a restricted range of binary orbits that describe the observed timing data very well and explain the extraordinary and unique growth in the observed rate of rotational slow-down. This and the coincidence of the precise positions of the pulsar and MT91 213 all point to the two stars comprising a neutron-star/Be-star binary system in the Cyg OB2 association. Near periastron, the pulsar may well become obscured in the radio by free-free absorption and by severe pulse scattering in the Be-star circumstellar wind and disk as it passes behind the star, in the same way as the Be star SS 2883 obscures PSR B125963 at periastron (Johnston et al., 1996). However, studies of PSR J2032+4127’s DM and rotation measure (RM) variations are likely to present a rare opportunity to study the density of the stellar wind and any circumstellar disk and the magnetic field of the Be star. Although it may become obscured in the radio for a short while, it should be possible to track its rotation and orbit around periastron by observation in rays. Close to periastron, the MT91 213 Be star may also display optical and/or X-ray variability as a direct result of its interaction with the pulsar. Even now MT91 213 displays X-ray variability (e.g. Rauw et al. 2014) and optical variability (e.g. Salas et al. 2014; see also Camilo et al. 2009), although these are not atypical of the intrinsic variability observed in other Be stars and almost surely have nothing to do with the neutron star.
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.
Pulsar research at JBCA is supported by a Consolidated Grant from the UK Science and Technology Facilities Council (STFC).
The GBT is operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.
We are grateful to Jules Halpern for illuminating discussions.
- Abdo et al. (2009) Abdo A. A. et al., 2009, Science, 325, 840
- Abdo et al. (2013) Abdo A. A. et al., 2013, ApJS, 208, 17
- Aliu et al. (2014) Aliu E. et al., 2014, ApJ, 783, 16
- Bosch-Ramon et al. (2012) Bosch-Ramon V., Barkov M. V., Khangulyan D., Perucho M., 2012, A&A, 544, A59
- Butt et al. (2006) Butt Y. M., Drake J., Benaglia P., Combi J. A., Dame T., Miniati F., Romero G. E., 2006, ApJ, 643, 238
- Camilo et al. (2009) Camilo F. et al., 2009, ApJ, 705, 1
- Espinoza et al. (2011) Espinoza C. M., Lyne A. G., Stappers B. W., Kramer M., 2011, MNRAS, 414, 1679
- Hanson (2003) Hanson M. M., 2003, ApJ, 597, 957
- Hobbs et al. (2010) Hobbs G., Lyne A. G., Kramer M., 2010, MNRAS, 402, 1027
- Hobbs et al. (2006) Hobbs G. B., Edwards R. T., Manchester R. N., 2006, MNRAS, 369, 655
- Hohle et al. (2010) Hohle M. M., Neuhäuser R., Schutz B. F., 2010, Astronomische Nachrichten, 331, 349
- Horns et al. (2007) Horns D., Hoffmann A. I. D., Santangelo A., Aharonian F. A., Rowell G. P., 2007, A&A, 469, L17
- Johnston et al. (1992) Johnston S., Manchester R. N., Lyne A. G., Bailes M., Kaspi V. M., Qiao G., D’Amico N., 1992, ApJ, 387, L37
- Johnston et al. (1996) Johnston S., Manchester R. N., Lyne A. G., D’Amico N., Bailes M., Gaensler B. M., Nicastro L., 1996, MNRAS, 279, 1026
- Kargaltsev et al. (2014) Kargaltsev O., Pavlov G. G., Durant M., Volkov I., Hare J., 2014, ApJ, 784, 124
- Kiminki et al. (2007) Kiminki D. C. et al., 2007, ApJ, 664, 1102
- Lorimer et al. (2006) Lorimer D. R. et al., 2006, MNRAS, 372, 777
- Lyne (2013) Lyne A., 2013, in IAU Symposium, Vol. 291, IAU Symposium, van Leeuwen J., ed., pp. 183–188
- Lyne et al. (2010) Lyne A., Hobbs G., Kramer M., Stairs I., Stappers B., 2010, Science, 329, 408
- Lyne (2008) Lyne A. G., 2008, in American Institute of Physics Conference Series, Vol. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More, Bassa C., Wang Z., Cumming A., Kaspi V. M., eds., pp. 561–566
- Martins et al. (2005) Martins F., Schaerer D., Hillier D. J., 2005, A&A, 436, 1049
- Massey & Thompson (1991) Massey P., Thompson A. B., 1991, AJ, 101, 1408
- Murakami et al. (2011) Murakami H., Kitamoto S., Kawachi A., Nakamori T., 2011, PASJ, 63, 873
- Rauw et al. (2014) Rauw G. et al., 2014, ApJS, in Press
- Ray et al. (2011) Ray P. S. et al., 2011, ApJS, 194, 17
- Rygl et al. (2012) Rygl K. L. J. et al., 2012, A&A, 539, 79
- Salas et al. (2014) Salas J., Maíz Apellániz J., Barbá R. H., 2014, ArXiv e-prints (1410.6767)