Discovery of extended and variable radio structure from the gamma-ray binary system PSR B1259-63/LS 2883

PSR B1259-63 is a 48 ms pulsar in a highly eccentric 3.4 year orbit around the young massive star LS 2883. During the periastron passage the system displays transient non-thermal unpulsed emission from radio to very high energy gamma rays. It is one of the three galactic binary systems clearly detected at TeV energies, together with LS 5039 and LS I +61 303. We observed PSR B1259-63 after the 2007 periastron passage with the Australian Long Baseline Array at 2.3 GHz to trace the milliarcsecond (mas) structure of the source at three different epochs. We have discovered extended and variable radio structure. The peak of the radio emission is detected outside the binary system near periastron, at projected distances of 10-20 mas (25-45 AU assuming a distance of 2.3 kpc). The total extent of the emission is ~50 mas (~120 AU). This is the first observational evidence that non-accreting pulsars orbiting massive stars can produce variable extended radio emission at AU scales. Similar structures are also seen in LS 5039 and LS I +61 303, in which the nature of the compact object is unknown. The discovery presented here for the young non-accreting pulsar PSR B1259-63 reinforces the link with these two sources and supports the presence of pulsars in these systems as well. A simple kinematical model considering only a spherical stellar wind can approximately trace the extended structures if the binary system orbit has a longitude of the ascending node of omega ~ -40 deg and a magnetization parameter of sigma~0.005.


INTRODUCTION
The binary system PSR B1259−63/LS 2883 1 is formed by a young 48 ms radio pulsar in an eccentric orbit of 3.4 years around a massive main-sequence star (Johnston et al. 1992. The parameters of the system are shown in Table 1. The spectral type of the massive star, O9.5 Ve, and some of the binary parameters have been recently updated by Negueruela et al. (2011), who obtained a distance to the system of 2.3 ± 0.4 kpc. Close to the periastron passage the system displays non-thermal unpulsed emission that has been detected in radio (Johnston et al. 2005), X-rays (Cominsky et al. 1994;Uchiyama et al. 2009;Chernyakova et al. 2009), hard X-rays up to 200 keV (Grove et al. 1995), and very high energy (VHE; 0.1-100 TeV) γ-rays above 380 GeV (Aharonian et al. 2005). In the range ∼0.1-100 GeV strict upper limits were obtained by EGRET close to the 1994 periastron passage (Tavani et al. 1996), and the source has not yet been observed by AGILE and Fermi close to periastron. The VHE emission is variable on orbital timescales, and is interpreted as the result of inverse Compton upscattering of stellar UV photons by relativistic electrons, which are accelerated in the shock between the relativistic wind of the young non-accreting pulsar and the wind of the stellar companion (see Maraschi & Treves 1981;Tavani & Arons 1997;Kirk et al. 1999;Dubus 2006;Bogovalov et al. 2008, and references therein).
The high-mass binaries LS 5039 and LS I +61 303 have also 1 Star 2883 in the catalog of Luminous Stars in the Southern Milky Way (Stephenson & Sanduleak 1971). The use of SS 2883 should be avoided, see Negueruela et al. (2011). been detected at VHE, and show a broadband spectral energy distribution (SED) similar to that of PSR B1259−63/LS 2883 (Dubus 2006). In contrast to X-ray binaries, all three sources have the peak of the SED at MeV-GeV energies. For these reasons, they can be considered gamma-ray binaries. However, the nature of the compact objects in LS 5039 and LS I +61 303 is unknown because their masses are not well constrained by the system mass functions (Casares et al. 2005a,b;Aragona et al. 2009), and no pulsations have been found. These systems have been extensively observed at VHE during several orbital cycles, while the observations of PSR B1259−63, the only system with a confirmed pulsar, are scarce due to the long orbital period. Any observational link between the three gamma-ray binaries would shed light in the understanding of this kind of systems.
The radio emission from the PSR B1259−63 system is described most recently in Johnston et al. (2005), which includes multiwavelength data from ATCA observations obtained during the periastron passages of 1994, 1997, 2000, and 2004 (hereafter we will refer to the epoch of each periastron passage as τ ). The emission has two non-thermal components. The first one is pulsed emission with a flux density of ∼2-5 mJy at 2.5 GHz and a nearly flat spectral index, which disappears approximately from 16 days prior to periastron (τ −16) to 15 days after periastron (τ +15). This has been interpreted as an eclipse of the pulsar when it crosses behind the equatorial circumstellar disk present around the massive star (Melatos et al. 1995). The second component is transient unpulsed synchrotron emission that appears at τ −20 days and shows two peaks centered around τ −10 and τ +20 days, with flux densities at 2.5 GHz up to ∼15-20 mJy and ∼30-50 mJy, respectively. After the post-periastron peak, the flux density of the unpulsed emission decreases continuously, and it has been detected up to ∼ τ +100 days. This transient emission remains optically thin during the outbursts.
The broadband transient emission of PSR B1259−63/ LS 2883 is produced around periastron, when strong interaction is produced between the stellar and pulsar winds. The shocked material is contained by the stellar wind behind the pulsar, producing a nebula extending away from the stellar companion. Along this adiabatically expanding flow, the accelerated particles produce synchrotron emission from radio to X-rays (Tavani & Arons 1997;Kirk et al. 1999;Dubus 2006;Takata & Taam 2009). The expected morphology depends on the magnetization parameter of the pulsar wind, σ, defined as the upstream ratio of magnetic to kinetic energy. These models predict that the radio emission extends up to several AU (corresponding to several milliarcseconds, mas, at 2.3 kpc), and that its structure should be variable on orbital timescales. This radio behavior has not been tested in PSR B1259−63, which only emits over a period of a few months every 3.4 year.
Here we present the first VLBI radio images of PSR B1259−63, obtained close to its 2007 periastron passage. The high-resolution images at three different orbital phases provide a direct view of the small-scale morphology of the source, which is comparable to those previously observed in LS 5039 and LS I +61 303.  Table 2. The small number of antennas provides a rather poor uv-coverage that makes the data calibration and imaging more difficult than for arrays with more elements, because we have less baseline redundancy. The data were recorded at a bit rate of 512 Mbps per telescope distributed in eight sub-bands (four for each right-and left-handed polarization) with a bandwidth of 16 MHz, each of them correlated using 64 frequency channels, two-bit sampling, and 2 s of integration time. Hobart and Ceduna recorded at 256 Mbps (only two sub-bands for polarization). The data were correlated at Swinburne University using the DiFX software correlator (Deller et al. 2007) without applying pulsar gating or binning.
The observations were performed using phase referencing on the calibrator J1337−6509 (B1334−649), which has an angular separation of 4. • 0 from PSR B1259−63 and was correlated at α J2000.0 = 13 h 37 m 52. s 4443 and δ J2000.0 = −65 • 09 ′ 24. ′′ 900. This reference position from the global VLBI solution 2006d_astro 2 has an uncertainty of 13 mas. The cycle time was 6 minutes, spending half of the time on the phase calibrator and the target source alternatively. The total flux density of the phase calibrator was 457 ± 11, 367 ± 20, and 491 ± 6 mJy for runs A, B, and C, respectively. The source J0538−4405 (B0537−441) was used as a fringe finder for runs A and B, and J1337−6509 (B1349−439) was used for run C. No astrometric check source was observed during the runs.   (1), from the fourth orbital solution obtained in Wang et al. (2004).
The data reduction was performed in AIPS 3 . Total electron content models based on GPS data obtained from the CD-DIS data archive 4 were used to correct phase variations due to the ionosphere. Standard instrumental corrections were applied (parallactic angle, instrumental offsets, and slopes between and within bands and bandpasses). Fringe fitting on the phase calibrator was performed with the AIPS task FRING, and the solutions were applied to the target source. A selfcalibrated model of the calibrator was used as an input model for CALIB, which was used to reduce the amplitude instabilities on timescales greater than 10 minutes. The data were averaged in frequency and time, and clean images were produced. For run B only, a single round of phase self-calibration was applied, to mitigate the residual atmospheric phase instabilities, which were more noticeable at this epoch. The final images were produced using a natural weighting scheme (robust 5 within the AIPS task IMAGR). For run B, robust 2 and tapering were applied to avoid the presence of possible unreliable high-resolution features due to sidelobes of the synthesized beam. Self-calibration slightly affects measured properties such as extent and position angle (P.A.), while it preserves the morphology.  Table 3 should be taken as uncertain at the ∼ 10% level. The flux densities of runs A and B are compatible with previous ATCA observations at the corresponding orbital phases (the ATCA data from the current observations were not correlated as an independent array). The flux density in run C is compatible with the flux density of the pulsar. This is expected considering the lack of unpulsed emission at τ + 150 and τ + 180 in previous ATCA observations (Johnston et al. 2005).
The phase-referenced observations allow us to obtain relative astrometry between runs. Since no astrometric check source was observed, we do not have real measurements of the astrometric uncertainties, and we will use the formal errors of a Gaussian fit obtained with JMFIT within AIPS. Given the relatively large calibrator-target separation of ∼ 4 • , we expect an additional systematic component to the position error due to the unmodeled ionosphere of 1-5 mas (Brisken et al. 2000). As a reference position for the plots we use the peak position of run C, α J2000.0 = 13 h 02 m 47. s 6435(1) and δ J2000.0 = −63 • 50 ′ 08. ′′ 636(1), which we consider to represent the pulsar position at MJD 54623.48. Assuming a mass of the neutron star of 1.4 M ⊙ and a stellar mass of 31±5 M ⊙ (see Table 1 for the system parameters), the mass function provides  0.0 ± 0.6 0.0 ± 1.1 · · · · · · NOTE. -The columns are run label, total and peak flux density at 2.3 GHz, synthesized beam (HPBW) parameters, size and P.A. of the extended emission, the position offset of the peak of the emission from the reference position (position of run C), and the range of possible separations between the peak and the pulsar position. i = 22. • 2 and the semimajor axis of the pulsar orbit is 7.2 AU, or 3.1 mas for a distance to the system of 2.3 kpc. The red crosses in Figure 1 mark the region where the pulsar should be located in each epoch. Their centers are placed at the position of run C corrected for proper motion at the corresponding epochs (MJD 54309.25 for run A and MJD 54329.18 for run B). Since we do not know the orientation of the orbit in the sky, we plot as error bars the projected orbital separation of the pulsar with respect to run C. The error bars also include the 1σ uncertainties on the mass of the star (and hence the distance uncertainty), on the astrometry of run C, and on the offset due to the proper motion. Finally, we have included the astrometry in runs A and B to compute the range of possible separations between the peak of the emission and the pulsar, as shown in Table 3. We consider all possible values of the longitude of the ascending node (Ω), which determines the orientation of the orbit in the plane of the sky.

KINEMATICAL INTERPRETATION
The radio morphology at a given epoch depends on the spatial distribution of synchrotron emitting particles and their emission processes. Given the limitations of our data (only two images, and without accurate astrometry), we have used a simple kinematical model to check if it can trace the extended structures detected.
We have considered the shock between the relativistic pulsar wind and a spherical stellar wind. The shock is produced at the standoff distance, the region where the pulsar and stellar wind pressures balance, as described in Dubus (2006). The evolution of the nebular flow after the shock is described in Kennel & Coroniti (1984). This should be considered as a first approximation to the much more complex hydrodynamic behavior of a shocked flow in a binary system, as shown in Bogovalov et al. (2008). In the Kennel & Coroniti approximation of a non-turbulent adiabatically expanding flow, the flow speed depends only on the magnetization parameter σ when assuming σ ≪ 1. This allows us to compute the past trajectory of the flow produced behind the standoff distance, which depends on the components separation along the orbit, the mass loss rate of the star (0.6 × 10 −7 M ⊙ yr −1 for a typical O9 star; Vink et al. 2000), the terminal wind velocity, v ∞ , and the spin-down luminosity of the pulsar,Ė sp . With these restrictions, the only free parameters are the longitude of the ascending node, Ω, which describes the orientation of the orbit, and the magnetization parameter, σ. Projecting the past trajectory of the flow on our VLBI images of runs A and B, we found that the best match with the detected morphologies is obtained for Ω ≃ −40 • and a magnetization parameter of σ ≃ 0.005, assuming an orbital inclination of i = 22. • 2. The obtained trajectories are shown in Figure 2 for three different values of σ. We note that the uncertainty in the distance to the system scales the size of the contours, but not the orbit and the trajectories, which are computed in AU. The range of magnetizations in the plots has been chosen to approximately show the effect of keeping σ = 0.005 and changing either the distance from 1.9 to 2.7 kpc, either Ω from −35 to −45 • , or a variation of the productṀ v ∞ of two orders of magnitude.
The parameters of the circumstellar equatorial disk of the star are uncertain, but consideringṀ e = 5 × 10 −8Ṁ yr −1 , (Johnston et al. 1996), and a wind velocity of ∼10 km s −1 , the productṀ e v ∞,e is ∼100 times smaller than the spherical wind contribution. If this equatorial component dominates during certain orbital phases, the flow trajectory would be closer to the orbit than the dashed models in Figure 2 for some specific regions of the more recent part of the trajectory. However, we do not attempt to accurately model this circumstellar disk, as its density, velocity, and crossing time are unknown.
We note that the astrometric errors can be of the order of ∼10 AU (see above). We also emphasize that this simple model cannot account for the complex magnetohydrodynamical turbulences of the real flow and should be considered as a first approximation to constrain σ. Previous studies have assumed a magnetization parameter σ of around 0.01-0.02 for this pulsar (see Tavani & Arons 1997;Dubus 2006), considerably higher than our best-fit value, although comparable lower values for σ have been suggested for other systems (e.g., the Crab pulsar; Kennel & Coroniti 1984).

DISCUSSION AND CONCLUSIONS
The results presented in this Letter show that the particle accelerator within PSR B1259−63/LS 2883 can produce a flow of particles emitting synchrotron radiation that can travel several AU. The total projected extent of the nebula is ∼ 50 mas, or 120 ± 20 AU, and the peak of the emission is clearly displaced from the binary system orbit (see Figure 1 and Table 3). Similar morphologies and displacements have been found in the other two known gamma-ray binaries, LS 5039 and LS I +61 303, although for smaller sizes and on shorter timescales Moldón et al. 2010;Massi et al. 2004;Dhawan et al. 2006;Albert et al. 2008).
There is an ongoing debate on the nature of the compact object and particle acceleration mechanisms in LS 5039 and LS I +61 303. Although initially suggested to be accreting/ejecting microquasar systems (Paredes et al. 2000;Massi et al. 2004), they are now thought to contain young non-accreting pulsars (Maraschi & Treves 1981;Dubus 2006;Dhawan et al. 2006;Ribó et al. 2008;Moldón et al. 2008), which can explain their multiwavelength emission (Sierpowska- Bartosik & Torres 2008;Cerutti et al. 2008;Bogovalov et al. 2008;Zdziarski et al. 2010). However, there are three basic issues that are not well understood for LS 5039 and LS I +61 303, and comparisons with PSR B1259−63/LS 2883 can help to clarify the situation. First, the putative pulsar properties of these two sources are unknown, as no pulsations have been detected for these systems. The lack of pulsations can be explained by the intense stellar wind that produces an extremely high absorption for these small orbits, a 3.9 day orbit with separations between 0.1 and 0.2 AU for LS 5039, and 26.5 days with separations of 0.1-0.7 AU for LS I +61 303 (Casares et al. 2005a,b;Aragona et al. 2009). As a reference, the separation for PSR B1259−63 is in the range 0.9-13.4 AU, and the pulsations disappear for distances below ∼ 1.6 AU. Second, it is not clear if the massive stellar wind can confine the pulsar wind (Romero et al. 2007). VLBI images, like the ones presented here, can shed light on the shock conditions and geometry. Third, the observed SED and variability at GeV energies is not well understood (Abdo et al. 2009a,b;Torres 2010). Fermi and AGILE are observing PSR B1259−63 for the first time during the 2010 periastron passage, and are providing GeV data that can be compared with those already obtained for LS 5039 and LS I +61 303. In this context, the high-resolution VLBI radio observations presented here establish a common link to test the similarities between the three systems.
In conclusion, our results provide the first observational evidence that non-accreting pulsars orbiting massive stars can produce variable extended radio emission at AU scales. Similar structures are also seen in LS 5039 and LS I +61 303, in which the nature of the compact object is unknown because the detection of pulsations is challenging. The discovery presented here for the young non-accreting pulsar PSR B1259−63 reinforces the link with these two sources and supports the presence of pulsars in these systems as well. Planned LBA observations of PSR B1259−63 during the 2010 periastron passage will allow us to provide a full comparison with the behavior observed in these sources, which have been extensively monitored during several orbital cycles. We have also shown that the orientation of the orbit and the magnetization of the pulsar can be inferred from VLBI observations of the source. Several images at different orbital phases covering a wider range of true anomalies will allow for a complete modeling of the orbital changes of the extended emission. Finally, accurate VLBI observations of the pulsed emission during several orbits can provide the pulsar trajectory, from which we can directly obtain the proper motion of the binary system, the inclination and Ω of the orbit, and the distance to the system.