Free Access
Issue
A&A
Volume 548, December 2012
Article Number A103
Number of page(s) 15
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201219553
Published online 30 November 2012

© ESO, 2012

1. Introduction

Gamma-ray binaries are extreme systems that produce non-thermal emission from radio to very high energy (>TeV) gamma rays, with the energy output in the spectral energy distribution (SED) dominated by the MeV–GeV photons (see Paredes 2011, for a review). Their broadband emission is usually modulated by the orbital cycle of the system, which suggests that the physical conditions are also periodic and reproducible. The wide range of different orbital periods and eccentricities of the known gamma-ray binaries (see Table 4 in Casares et al. 2012) provides a diversity of different ambient conditions in which the physical processes take place. The diversity of systems, together with the reproducibility of the conditions within each system, makes gamma-ray binaries excellent physical laboratories in which high-energy particle acceleration, diffusion, absorption, and radiation mechanisms can be explored. Nevertheless, the number of known gamma-ray binaries is still very limited. Five binary systems have been classified as gamma-ray binaries and present orbital modulation of the GeV and/or TeV gamma-ray emission: PSR B1259−63 (Aharonian et al. 2009), LS 5039 (Aharonian et al. 2006; Abdo et al. 2009b), LS I +61 303 (Albert et al. 2009; Abdo et al. 2009a), HESS J0632+057 (Aharonian et al. 2007; Maier et al. 2011), and 1FGL J1018.6−5856 (Fermi LAT Collaboration et al. 2012). The latter has not been clearly detected at TeV energies yet. On the other hand, Cygnus X-3 presents active periods of GeV emission in which orbital periodicity is detected (Fermi LAT Collaboration et al. 2009; Tavani et al. 2009). Cygnus X-1 showed evidence of gamma-ray emission at TeV energies on a single night of observations (Albert et al. 2007). However, in these two cases the peak of the SED is at keV energies and therefore they are classified as X-ray binaries, the high-energy activity of which is clearly powered by accretion. Another phenomenological difference between X-ray and gamma-ray binaries is that variations in the non-thermal broadband emission in X-ray binaries are usually linked to X-ray state changes, whereas in gamma-ray binaries it is coupled to the orbital period.

LS 5039 is a binary system consisting of the homonym bright star LS 5039, of spectral type ON6.5 V((f)), and a compact object (Casares et al. 2005). The orbital period of the system is 3.9 days and the eccentricity of the orbit is e ~ 0.35 (Casares et al. 2005; Aragona et al. 2009; Sarty et al. 2011). Given the mass function of the binary system, the compact object would be a black hole for low inclinations of the orbit, and a neutron star for high inclinations (see the discussion in Casares et al. 2005; and recent constraints in Casares et al., in prep.). No short-period pulsations were found that could demonstrate the presence of a pulsar either in radio (McSwain et al. 2011) or X-rays (Martocchia et al. 2005; Rea et al. 2011). The distance to the system has recently been updated to 2.9 ± 0.8 kpc Casares et al. (in prep.). LS 5039 displays non-thermal, persistent, and periodic gamma-ray emission up to 4 TeV (Paredes et al. 2000; Aharonian et al. 2005; Abdo et al. 2009b). The X-ray and gamma-ray emission show a periodic modulation of 3.9 days (Bosch-Ramon et al. 2005; Aharonian et al. 2006; Kishishita et al. 2009; Abdo et al. 2009b). The total radio emission of LS 5039 is variable, although no periodicity or strong outbursts have been detected so far (Martí et al. 1998; Ribó et al. 1999, 2002). The radio emission above 1 GHz is non-thermal, with a spectral index of  − 0.46 (Martí et al. 1998) that is inverted at lower frequencies (Godambe et al. 2008; Bhattacharyya et al. 2012; but see Pandey et al. 2007). The synchrotron radio emission appears to be extended when observed with VLBI on scales of milliarcseconds (mas). The source shows a main core and extended bipolar emission that has been observed in directions with position angles (PA) between 120 and 150° and projected angular distance between 1 and 180 mas (3–500 AU) from the core (Paredes et al. 2000, 2002; Ribó et al. 2008). In particular, Ribó et al. (2008) measured a change in the source morphology at scales of 2–6 mas from VLBI images obtained five days apart. Finally, for X-rays, Durant et al. (2011) discovered an extended component in X-rays up to 2′ from LS 5039.

Other gamma-ray binaries also present extended and variable radio emission at mas scales. PSR B1259 − 63 is the only gamma-ray binary in which radio pulsations have been detected up to now. During the periastron passage of this eccentric system, which occurs every 3.4 yr, the close interaction between the pulsar and the massive Be star and/or Be disc produces broadband emission during a few months. Moldón et al. (2011a) discovered transient extended emission during and after the periastron passage with a total projected extent of  ~50 mas (120 AU). The peak of the radio emission was detected outside the binary system. VLBI observations at different orbital phases of the 26.5-d period binary LS I +61 303, which also hosts a Be star, were conducted by Dhawan et al. (2006), showing an orbital variability that these authors interpreted as the signature of a cometary tail produced in the colliding winds scenario. Another interpretation of the data suggests that the changes are compatible with a precessing microblazar (Massi et al. 2012). The system HESS J0632+057, also hosting a Be star, has an orbital period of 321 days (Bongiorno et al. 2011), and faint radio extended emission at mas scales was detected  ~130 days after the periastron passage (Moldón et al. 2011b). For 1FGL J1018.6 − 5856, with a period of 16 d and also hosting an O6V((f)) star (Fermi LAT Collaboration et al. 2012), no VLBI observations have been reported so far. In addition, some of these systems also present extended X-ray emission at larger scales of  ~1′, as hinted at for LS I +61 303 (Paredes et al. 2007) and shown for PSR B1259 − 63 (Pavlov et al. 2011).

Several theoretical models were developed to explain the multiwavelength behaviour of LS 5039. The very high energy gamma-ray emission can be interpreted as the result of inverse Compton upscattering of stellar UV photons by relativistic electrons. The acceleration of electrons can be explained by two exclusive scenarios: acceleration in the jet of a microquasar powered by accretion (Paredes et al. 2006; Bosch-Ramon et al. 2006; and the review in Bosch-Ramon & Khangulyan 2009), or shocks between the relativistic wind of a young non-accreting pulsar and the wind of the stellar companion (Maraschi & Treves 1981; Tavani & Arons 1997; Dubus 2006; Khangulyan et al. 2007). A simple and shockless microquasar scenario is disfavoured by previous VLBI radio observations of LS 5039 (Ribó et al. 2008). Also, no signs of an accretion disc have been detected so far (although see Barkov & Khangulyan 2012; Okazaki et al. 2008, for alternative explanations). Sarty et al. (2011) used the stability of the optical photometry of LS 5039 to constrain the orbit inclination to be below 60° and the mass of the compact object above 1.8 M. Models based on the young non-accreting pulsar scenario, first proposed for LS 5039 in Martocchia et al. (2005), had provided descriptions of the HE/VHE light curves and the spectral evolution of the source as a function of the orbital phase (see Sierpowska-Bartosik & Torres 2007; Dubus et al. 2008; Khangulyan et al. 2008; Takahashi et al. 2009). Apart from the acceleration of particles in the main shock between winds, additional gamma-ray emission can be produced by the unshocked pulsar wind (Cerutti et al. 2008) and by secondary cascading (Bosch-Ramon et al. 2008; Cerutti et al. 2010). The X-ray light curve and spectrum show orbital variability, although there are no signatures of variable X-ray absorption or X-ray occultations, as discussed in Reig et al. (2003), Martocchia et al. (2005), Bosch-Ramon et al. (2007) and Szostek & Dubus (2011). Zabalza et al. (2011) constrained the spin-down luminosity of the putative pulsar to be  ~3−6 × 1036 erg s-1.

Very recently, detailed hydrodynamical simulation have been obtained to understand the wind shocks and the flow dynamics in gamma-ray binaries. Two-dimensional hydrodynamic simulations were performed to study the wind-wind collision interaction at scales of the binary system for PSR B1259 − 63 (Bogovalov et al. 2008, 2012) and for LS I +61 303 (Romero et al. 2007). Also for PSR B1259 − 63, Okazaki et al. (2011) and Takata et al. (2012) presented 3D simulations of the tidal and wind interactions. Lamberts et al. (2012) considered the general case of colliding wind binaries and described the outflow structure at larger distances. Finally, the flow structure for a case similar to LS 5039 was described in Bosch-Ramon et al. (2012) for scales up to 100 times the orbit size.

From a theoretical point of view, a general discussion on the properties of the radio emission of LS 5039 can be found in Bosch-Ramon et al. (2008). However, only few predictions on the expected radio morphology have been discussed. Dubus (2006) presented emission maps at different orbital phases for LS 5039 and predicted orbital morphological changes at mas scales, as well as displacements of the peak of the emission of a few mas. Bosch-Ramon & Khangulyan (2011) obtained synthetic maps at radio wavelengths from a model based on secondary particles created from gamma-ray absorption, which may account for a significant fraction of the total radio flux density of the source and also predicts extended radio emission at mas scales.

In this paper we present the results from a dedicated VLBA campaign that covers a full orbital cycle and provides astrometric and morphological information during five consecutive days. The observations and data reduction are presented in Sects. 2 and 3, including a discussion on the calibration caveats. In Sect. 4 we present the obtained astrometry and morphology of the source at different orbital phases and analyse the components characterising the extended emission of the source. In Sect. 5 we show a compilation of VLBI data on LS 5039 and we describe and compare the morphology at different orbital phases for observations at frequencies between 1.7 and 8.5 GHz taken from 1999 to 2009. In Sect. 6 we present a model to check if the measured changes are compatible with a wind-wind collision scenario. The model allows us to estimate the orientation of the orbit on the sky, in particular its inclination, and consequently constrain the mass of the compact object. Finally, we present a summary of the results and the main conclusions of this work in Sect. 7.

2. Observations

We conducted VLBI observations on LS 5039 during five consecutive days to cover an orbit of 3.9 d and an extra day to disentangle between orbital or secular variability. We observed at 5 GHz (6 cm wavelength) with the Very Long Baseline Array (VLBA) of the National Radio Astronomy Observatory (NRAO). The stations forming the array are Br, Fd, Hn, Kp, La, Mk, Nl, Ov, Pt, and Sc. The VLBA project code is BR127, and the five observations were conducted from July 5 to 9, 2007 (MJD 54 286 to 54 290). The five runs, hereafter run A–E, spanned from 03:30 to 09:30 UTC and were scheduled in the same way to reduce differences between runs. The corresponding orbital phases were computed using the ephemerides in Casares et al. (in prep.), T0 = HJD 2 453 478.09(6) and Porb = 3.90603(8), where the values in parentheses refer to the uncertainty in the last digit. The phase of the periastron passage is 0.0. The orbital phase at the centre of each observation was 0.03, 0.29, 0.55, 0.80, and 0.06, respectively, with an uncertainty of 0.02, or approximately 1.9 h. A log of the observations is shown in Table 1.

Table 1

Log of VLBA observations of project BR127, conducted at 5 GHz.

The data were obtained with single circular left-handed polarisation, with eight sub-bands of 8 MHz, and were correlated with two bits per sample, obtaining a total aggregate bit-rate of 256 Mbps. The data were processed at the VLBA hardware correlator in Socorro, which produced the final visibilities with an integration time of 2 s.

We conducted the observations using the phase-referencing technique on the phase reference calibrator J1825−1718, located at 2.5° from LS 5039 in a PA of  −176.4° (see Fig. 1). We observed in 4.2-min cycles, spending 2.5 min on the target source and 1.7 min on J1825−1718. The total scheduled time on the target source was 2.7 h. As an astrometric check source we used J1818−1108, including one scan of 2.5 min every 30 min. Finally, we included three 5-min scans on the fringe finder J1733−1304. The reference position for the observations is the correlation position of the phase reference source J1825−1718 (or  ± 1.3 mas) and (or  ± 2 mas), taken from Fomalont et al. (2003) via the SCHED database (NASA catalogue 2005f_astro). LS 5039 is close to the galactic plane, at a galactic latitude of  −1.29°. The phase calibrator and the astrometric check source suffer from galactic scatter broadening (a general discussion on this effect can be found in Fey et al. 1991). The phase calibrator J1825−1718 has a total flux density of 350 mJy at 5 GHz, but the amplitude of the visibilities decreases with the uv distance and the phase calibration fails for baselines longer than  ~90–100 Mλ. The compact core has a size of 3 mas. The resolution of the phase-referenced images is therefore limited by the scatter broadening of the phase calibrator. On the other hand, the astrometric check source J1818−1108 is not compact. The self-calibrated images show a main core of 680 mJy and a secondary blob of 180 mJy located at 50 mas westwards from the core. Both components have a size of  ~8 mas. The visibility amplitudes quickly drop with the uv distance and there is no signal beyond  ~30 Mλ. The distance between J1818 − 1108 and the phase calibrator is 6.4° in PA of 16° (see Fig. 1), which makes it difficult to transfer the phase solutions for this low-elevation source.

thumbnail Fig. 1

Distribution on the sky of LS 5039, the phase calibrator J1825 − 1718, and the astrometric check source J1818 − 1108. The dashed line indicates the Galactic Plane. The arrows indicate the two phase calibrations obtained: direct, using J1825 − 1718 as phase reference; inverse, using LS 5039 as phase reference.

3. Data reduction

Data reduction was principally performed in AIPS1. The Difmap package (Shepherd 1997) was used for imaging and self-calibration. The raw visibilities were loaded into AIPS, and a priori flagging on telescope off-source times because of antenna slewing was applied. We searched for data with instrumental problems and flagged them, as well as all visibilities with antenna elevation below 5°. Initially, we updated the geometrical model of the correlator using the Earth orientation parameters (EOPs). However, the correction implied a change on the final measured positions one order of magnitude below our final uncertainties. Considering the small effect and the problems with this correction in the past (wrong parameters during the correlations between 2003 and 20052 and a bug in the AIPS task CLCOR from 2009 to 2011), we did not include them in the final data processing. For these observations, the EOPs correction produced a  ~10° phase offset in the visibilities, which was in any case removed by the phase calibration (see below).

3.1. Ionospheric correction

The different ionospheric conditions above VLBI antennas are one of the main contributions to systematic errors for astrometry at low frequencies. The unmodelled phase and delay contribution of the ionosphere randomly modifies the observed phases on timescales of minutes, producing two effects. On one hand, it spreads the signal and broadens the source, producing an overall decrease of the detected signal-to-noise ratio (S/N). On the other hand, it introduces an unknown position offset. This correction is specially important in this case because the sources have low declination, which implies low elevations during the observations, and because the sources are separated mainly in the north-south direction (see Fig. 1), where the atmospheric conditions change more significantly.

We used ionospheric total electron content (TEC) models based on GPS data obtained from the CDDIS data archive3 to correct the phase variations caused by the ionosphere. These maps provide one TEC measure every two hours on a grid of 5° × 2.5° in terrestrial longitude and latitude. This is a coarse correction because the ionosphere changes occur on a shorter timescale, although it accounts for significant astrometric offsets. The ionospheric models are produced by different institutes, which estimate the TEC above different locations and provide electron distribution maps (IONEX files) that can be loaded into AIPS to apply the phase corrections. We compared the effects of the phase correction on the visibilities using the models provided by four institutes: the Jet Propulsion Laboratory (JPL), the Center for Orbit Determination in Europe (CODE), the ESOC Ionosphere Monitoring Facility (ESA), and the Universitat Politècnica de Catalunya (UPC). We tested the models by comparing the S/N and the stability of the position of the peak of the emission of LS 5039 in different short intervals of 1 and 3 h along the observations. We did not measure appreciable differences in the peak flux density when using the different models. The ionospheric models provide an average position offset of 0.02 and  −1.7 mas in right ascension and declination, respectively. The models are very similar, and the dispersion of these corrections are 0.11 and 0.2 mas, respectively. For each run we applied the model that provides offsets closer to the averaged value. The JPL model was used to correct the runs A and C, and the CODE model to correct runs B, D, and E.

We checked the stability of the correction for the observed times by computing the variability of the TEC provided by the GPS maps. The most variable stations were Sc (at the end of the observation the Sun was already rising) and Mk (the observations started before sunset). Comparing the mean variability of the TEC content for every run during the central hours of the project (from 4:00 to 8:00 UTC), we found that runs B, D, and E have a relative average variability of 18%, for run C it was 15%, and for run A 22%. This shows that during run A the ionosphere was, in general, more variable. This can in part justify the astrometric problems encountered for this run (see Sect. 3.3), while the rest of the weather conditions were similar for all runs.

3.2. Amplitude and phase calibration

ACCOR was used to fix the amplitudes in the cross-correlation spectrum from the VLBA correlator. The amplitude calibration was performed using the antenna gains and the system temperatures measured at each station in real time during the observations, using the APCAL procedure. The parallactic angle correction was automatically applied with VLBAPANG. To correct the instrumental offsets and slopes between and within the different bands, we ran FRING on 30 s of a scan of the fringe finder 3C 345 (manual phase-cal). To correct the dependence of the visibility amplitudes with the frequency we used the auto-correlation values from 3C 345 to smooth the bandpass shape of the amplitudes.

An initial fringe fitting was performed using the FRING routine in AIPS, using a point-source model on the phase calibrator J1825−1718. These data were exported to Difmap and were self-calibrated. We fitted the visibilities for each run with a circular Gaussian. The mean flux density fitted for the core component was 345 mJy, with a standard deviation between runs of 2.3 mJy (0.67%), while the calibrator structure did not change significantly. Therefore, we combined the five data sets of the calibrator and performed several iterations of imaging and self-calibration to create a master calibrator image. We ran FRING again, but now we removed the source structure contribution to the phase calibration by using the combined image of J1825−1718 as an input model. One solution for each scan for all IFs was searched within a delay and rate search windows of 80 nanoseconds and 20 mHz, respectively. The minimum allowed S/N for solutions was 8, providing 95.8% of good solutions. The 4.2% of the data without solutions mostly correspond to long baselines (above 90 Mλ). The solutions for every frequency band were corrected with MBDLY, which fits the phase delay between bands. The final phases were explored and flagged when necessary. Finally, we applied the calibration and flag tables, and we averaged all frequency channels within each IF to obtain three single-source files, one for LS 5039, one for the phase calibrator J1825−1718 and one for the astrometric check source J1818−1108.

We also obtained inverse phase referencing based on LS 5039 because the source can be self-calibrated at 5 GHz. We followed an analogue procedure to fringe-fit the LS 5039 data, although we did not produce a combined data set. The S/N cut was set to 5 because of the faintness of the source. We flagged the unstable phases and transferred the solutions from LS 5039 to J1825−1718 and J1818 − 1108, for which inverse astrometry was obtained. A sketch of the phase referencing schemes is shown in Fig. 1.

3.3. Phase referencing imaging

To reduce the systematic errors of the phase referencing produced by data taken when the sources were at low elevations, we only used the central hours of the runs. We produced images with time intervals of 30 and 60 min and inspected the image quality and the astrometric stability. We also produced images for the individual scans. The data obtained at the beginning and at the end of the observations are not consistent because we measured artificial displacements of the peak of the emission between  ~0.5 and 3 mas, while the source became elongated in random directions. This behaviour can be explained by the ionosphere affecting the data with low elevation. We note that the effect is more significant for run A, which suffered severe displacements of more than 6 mas in the time blocks for the first hour and the last two hours. After inspecting these systematic errors for different time blocks we only used the data in the time intervals when the source was reasonably stable, from 04:40 to 07:40 UT for runs B, C, and D, from 04:30 to 07:00 UT for run A, and from 06:00 to 08:30 UT for run E. Taking into account that the declination of the phase-reference calibrator and the target source are  −17.3° and  −14.8°, respectively, that the observations were centred at the culmination of the sources, and that the observations were conducted under normal weather conditions (without rain, strong winds, or snow), the time ranges used correspond to elevations of most of the antennas above  ~25°. Some antennas contributing to the longest baselines (Mk, Sc, Hn, and Br) were below this value during part of the observations, and the phase calibration was not succesful for them during these time intervals, altough this can also be caused by the intrinsic calibrator structure.

The phase-referenced data were imaged using task IMAGR with a pixel size of 0.2 mas. We used a weighting scheme with robust 0 as a compromise between angular resolution and sensitivity. Pradel et al. (2006) showed that removing the shorter baselines of the VLBA can decrease the astrometric systematic errors by a 15%. For LS 5039, the source was detected with S/N above 10 when imaging the baselines with uv distances between 15–20 and 100 Mλ. The range of minimum baselines was lowered in those cases in which secondary lobes were more important or when the S/N of the peak emission was too low. The same strategy was used to image the astrometric check source J1818−1108. In this highly resolved source we only used baselines with uv distances between 12 and 40 Mλ.

For the inverse phase referencing, where phases from LS 5039 were transferred to J1825−1718 and J1818−1108, a uv range of 20 to 60 Mλ was used for J1825−1718 and 12 to 60 Mλ for the check source. The maximum baseline of the latter is higher than in the direct phase referencing because of the smaller angular distance between J1818−1108 and the phase reference source in this case (see Fig. 1). In both cases we used a pixel size of 0.5 mas.

3.4. Self-calibration

To perform the self-calibration of the LS 5039 data, we loaded the individual data files into Difmap. We averaged the data with a binning size of 50 s. The visibilities with long baselines were down-weighted, using a Gaussian taper at radius of 80 Mλ. We performed several iterations of cleaning and phase self-calibration. For each of them, phase solutions were found, first for the most compact part of the core by imaging the data with a uniform weighting scheme and then solving for the more extended emission using a natural weighting cleaning. After each cycle, an amplitude self-calibration was obtained, each time with a shorter integration time, up to a minimum of 30 min. Finally, two hybrid images were produced for each run, one with natural weight (lower resolution and better sensitivity) and one with uniform weight (better resolution and lower sensitivity). The final images were exported to AIPS. Their rms noise and the total flux density were obtained by fitting the pixel flux density distribution of the image with IMEAN. The extended emission was characterised by fitting two or three Gaussian components to the images using the task JMFIT.

We also fitted Gaussian components to the calibrated uv data sets of LS 5039 (modelfit). This is an iterative process, where the initial values for the fit have to be fixed manually, and the final solution can depend on the starting parameters. Therefore, the model fitting results can be slightly subjective, in particular for this case, where the source structure is not clearly defined by individual components. Despite this caveat, the modelfit provides morphological information directly from the uv data and therefore is not limited by a synthesized beam because it is independent of the deconvolution procedure used. We obtained the modelfit of the uv data with the task UVFIT in AIPS. As a starting point, we used the solutions from JMFIT obtained from the naturally weighted images. Then we inspected several combinations of number of components and shape restrictions until a robust solution was found. As a parallel check, we performed another modelfit in Difmap to verify the consistency of the results. As both approaches provide similar results, we only present here the solutions obtained with AIPS.

4. Analysis and results

4.1. Astrometry

The average position of LS 5039 from the five runs of project BR127 is (or  ± 0.4 mas) and (or  ± 0.6 mas). This position is measured with respect to the correlation position of J1825−1718, which is listed in Sect. 2. The absolute astrometry in this project was used to obtain an accurate proper motion of LS 5039 (Moldón et al. 2012). The astrometric uncertainties cannot be obtained directly from the position fit because of the unknown ionospheric effects, and therefore we used the astrometric check source J1818−1108 to estimate them. The standard deviation of the peak positions of the check source was 2.3 mas in right ascension and 2.4 mas in declination. These deviations were converted into uncertainties in the determination of the position of LS 5039 by multiplying them by two correction factors that account for the different observing conditions for LS 5039 and J1818−1108, following the general theoretical astrometric precision for an interferometer (Thompson et al. 1986) (1)where the wavelength (λ) and the maximum baseline (B) can be represented by the synthesized beam size. First, the uncertainty was scaled taking into account the different resolutions of the images. Basically, the data of J1818−1108 were imaged with poorer resolution, only using short baselines (see Sect. 3.3). This factor was obtained as a ratio of the synthesized beam size of the LS 5039 observations to the mean size of the J1818−1108 beams. Second, the astrometric precision depends on the S/N of the detection. We scaled the uncertainty by the ratio of individual S/N of LS 5039, which was around 15, to the mean S/N of J1818−1108, which was 7. These two corrections yield a scaling factor applied to the dispersion of J1818−1108 of  ~0.19 and  ~0.25 in right ascension and declination, respectively. Finally, the mean uncertainties in the determination of the right ascension and declination of LS 5039 are 0.4 and 0.6 mas, respectively. The same procedure was used for the inverted astrometry.

Another approach to estimate the astrometric uncertainties is to use the expected theoretical systematic uncertainties computed using Eq. (2) in Pradel et al. (2006), although some assumptions considered in the general discussion in that paper are not met in these observations. The minimum theoretical uncertainties for LS 5039 using J1825−1718 as a reference source and assuming mean tropospheric conditions are 0.17 and 0.56 mas in right ascension and declination, respectively, which are compatible with the estimation quoted above. This indicates that the uncertainties in declination seem to be dominated by systematic uncertainties. On the other hand, the minimum theoretical uncertainties for J1818−1108 are 0.43 and 1.42 mas in right ascension and declination, respectively. The real dispersion measured, 2.3 mas in right ascension and 2.4 mas in declination, is significantly higher possibly because the extended nature of the source, the poor uv coverage, and the low S/N. Given these results, we have adopted the more conservative approach described above, although we note that the real uncertainties of the measurement are possibly a non-trivial combination of the systematic and nominal uncertainties.

In Table 2 we present the measured astrometry. We show the relative astrometry (with respect to the average position) for the check source J1818−1108 and the synthesized beam of the phase-referenced images and the relative astrometry for LS 5039. We also show the inverted astrometry obtained using LS 5039 as a phase reference source. The mean synthesized beam for J1818−1108 in the direct phase referencing is 6.5 × 3.3 mas2 at PA of 6°, and for the inverse phase referencing it is 6.9 × 2.7 mas2 at PA of 7°. The measured relative offsets of the peak of the emission of LS 5039 using the two methods are plotted in Fig. 2. Both approaches provide similar relative astrometry. The peak position of LS 5039 during run A is significantly displaced from the mean value. However, in run A, the check source shows a even stronger displacement in the same direction (see first row in Table 2). This fact, combined with the suspicious behaviour described in Sect. 3.3, makes clear that the astrometry of run A suffers from uncorrected ionospheric effects, and therefore we do not consider the peak displacement in run A as reliable. The only significant displacement is found between run B and run C, which corresponds to displacements of 1.0 ± 0.5 mas and  −2.0 ± 0.7 mas in right ascension and declination, respectively. Therefore, the total displacement between phase 0.29 and phase 0.55 is 2.3 ± 0.6 mas in PA of 154°, which, considering the limitations of the data, cannot be considered a robust sign of peak displacement.

Table 2

Relative astrometry of the five runs of project BR127, also plotted in Fig. 2.

thumbnail Fig. 2

Relative offsets of the peak of the emission of LS 5039 during five consecutive days with respect to the average position. For straightforward comparison, we plot the opposite value of the inverse astrometry. The run label is displayed close to each direct measure. The size of the binary system is indicated by the small orbit in the corner, which is plotted face-on and with an arbitrary orientation. The displacements and uncertainties are shown in Table 2.

4.2. Morphology

The self-calibration of the data improves the S/N of the phase-referenced images while preserving the main features of the extended emission, partially avoiding the problems with the phase calibration caused by the ionosphere, at the expense of losing the astrometric information. In the first two rows of Fig. 3 we show the final self-calibrated images for each run, produced with natural and uniform weight, respectively. The five images with natural weight have synthesized beams with sizes about 6.1 × 2.3 mas2 and PA of  −2° and total flux densities of 27.9, 26.1, 24.4, 28.5, and 26.4 mJy with an rms noise of the images of  ~0.06 mJy beam-1. On the other hand, the images with uniform weight have synthesized beams with sizes around 4.1 × 1.3 mas2 and PA of  − 4°. The total flux densities in the five images are 27.8, 26.7, 25.4, 28.5, and 27.0, with an rms noise of 0.08 mJy beam-1. The individual components fitted to the images with JMFIT are summarised in Table 3, where we list the flux density, relative position, and size of each component. The last two columns show the deconvolved size of the components, which describes to a certain degree the intrinsic size of the component, deconvolved from the synthesized beam size. We also list the Gaussian components fitted to the uv data (modelfit), which are not convolved with any beam. Most of the modelfit components are circular instead of elliptical. We plot these components in the last row of Fig. 3 convolved with an artificial circular beam with an area equal to the average synthesized beam of the rest of the images. As a guide, we plot a line towards the direction of the main extended component, which corresponds to PA of  −75, 123,  −61,  −52, and  − 74°, respectively.

thumbnail Fig. 3

VLBA images of LS 5039 at 5 GHz from project BR127, obtained during five consecutive days in July 2007. Each column, labelled with the run name and orbital phase, corresponds to one epoch (see Table 1). For each epoch, the first and second rows correspond to the self-calibrated images obtained with a natural and uniform weighting scheme, respectively. The third row shows the uv components fitted to the data, convolved with a Gaussian circular beam with an area equal to the average synthesized beam of all images. The lines show the direction of the main extended component. The restoring beams are plotted in the bottom-left corner of each panel. Dashed contours are plotted at  − 3 times the rms noise of each image and solid contours start at 3 times the rms and scale with 21/2. The Gaussian components fitted to the images are listed in Table 3.

All images show a main core component and extended emission with changing PA During run A, which was obtained soon after periastron (orbital phase  ~0), the source is extended westwards, although the core is better described with two components. This double core has an important contribution eastwards, at a distance of  ~0.5 mas. For the image obtained 24 h later, at orbital phase 0.29, a fast morphological change occurs and the main extended component is detected towards the south-east. Run C, obtained soon after apastron, displays bipolar extended emission. The emission towards the north-west consists of two components, although only one component could be fitted. The faint component towards the south-east could not be fitted to the data and therefore it is not listed in Table 3. Run D is well described with a main core and a bright component towards the north-west. Finally, the images from run E, also obtained shortly after the periastron passage, recover the same morphology as run A. This shows that the morphology of LS 5039 is periodic on consecutive orbital cycles.

Table 3

Gaussian components fitted to the images shown in Fig. 3.

We also divided each data set into two blocks to check for possible intraday variations. We obtained modelfit solutions and produced images for the divided uv data sets. However, it was not possible to obtain reliable morphological differences because of the very different uv coverage at the beginning and the end of each run.

5. Compilation of VLBI observations

With the aim to confirm the repeatability of the changing morphological structure of LS 5039 found in our five consecutive day campaign (project BR127), we compiled the available archival VLBI data of the source that provide an angular resolution similar to the one explored here. We compiled images from another eleven observations conducted at frequencies between 1.7 and 8.5 GHz between 1999 and 2009. A log of these observations is shown in Table 4. The orbital phases were computed with the ephemerides from Casares et al. (in prep.) and have uncertainties below 0.022.

There are three observations at 5 GHz (6 cm wavelength) from 1999 and 2000 that correspond to the projects BP051 and GR021A–B, which included the VLBA and the phased VLA array as an additional VLBI station. These results were already presented in Paredes et al. (2000) and Ribó et al. (2008), where the interested reader can find detailed information on the data and the data reduction process.

Between 2004 and 2006, three runs were observed with the VLBA at 8.5 GHz (3.6 cm wavelength) as part of a long-term astrometric project (PI: V. Dhawan). The observational codes are BD087G, BD105A, and BD105G. The data were correlated with a lower sensitivity provided by 128 Mbps, and the final resolution is slightly better than for the other projects thanks to the higher frequency of the observations. We produced three self-calibrated images with natural weight. The obtained synthesized beams and noises are shown in Table 4. More details on the data reduction can be found in Moldón et al. (2012). The total flux densities recovered for BD087G, BD105A, and BD105G were 18.6, 17.6, 20.2 mJy, respectively.

We also present three self-calibrated images from project EF018 (PI: R. Fender), which consists of three runs (A–C) at 5.0 GHz obtained with the European VLBI Network (EVN) in 2007. The runs were separated by two days. The antennas used were Effelsberg, Westerbork, Jodrell Bank (Lovell), Onsala, Medicina, Noto, Torun, Nanshan (Urumqi), Sheshan (Shanghai), and Hartebeesthoek. This array provides good sensitivity, although each run lasted only four hours and the elevation of the source was lower than in the previous projects. The longest baselines in these observations were from the stations in Nanshan (Urumqi) and Hartebeesthoek, although the uv coverage for those baselines is poor and the final resolution is lower. The total flux densities measured were 26.2, 19.1, and 27.4 mJy, respectively. The lower flux density and noise of epoch EF018B is probably an artefact of the amplitude self-calibration and not a real change in the source flux density.

In Fig. 4 we show the self-calibrated images of these projects. The images in the panels are ordered in increasing phase, from left to right and from top to bottom. We distributed them in four groups, which correspond to images with similar orbital phases and that show approximately a similar morphology. For a quick reference, the inset table included in Fig. 4 shows the dates of the images in the panels. Although they were obtained with different instruments, frequencies, and at very different epochs, the images in each group share common morphological trends. Each group can be approximately characterised with the morphology described in Sect. 4.2 for the runs BR127A–D.

5.1. Morphology at larger scales

We observed LS 5039 at 1.7 GHz (18 cm wavelength) with the EVN in 2009, project code EM074. This is a deep pointed observation on LS 5039, conducted without phase referencing. The total time on source was 5.8 h, and nine stations participated: Effelsberg, Westerbork, Cambridge (32 m), Jodrell Bank (Lovell), Onsala, Medicina, Noto, Torun, and Urumqi. At this frequency the resolution is lower, although the long east-west baselines provide good resolution in right ascension. We produced a uniform weighting scheme image with a synthesized beam of 26 × 5.8 mas2 at PA 7° and an rms noise level of 0.04 mJy beam-1. The total flux density of the source was 33.5 mJy at 1.7 GHz. The self-calibrated image is shown in Fig. 5. The source shows a morphology more similar to the images in group 2, although it was observed at orbital phase 0.46, which is closer to those of group 3. This can be explained by the longer lifetime of electrons emitting synchrotron radiation at lower frequencies.

Table 4

Log of the VLBI observations of LS 5039.

On the other hand, the self-calibrated image from project ER011 provides lower angular resolution than the images in Fig. 4 and it is not reproduced here because it is already published in Paredes et al. (2002). The image was obtained at orbital phase 0.39. The innermost region of the core shows extended emission towards the south-east, similar to the group 2 images in Fig. 4, although additional bipolar diffuse components extend up to  ~20 mas in opposite directions. Paredes et al. (2002) also presented a lower resolution image obtained with MERLIN showing extended emission with PA of 150° up to 170 mas. Finally, we note that lower resolution VLA images show some hints of extended emission in the same north-west/south-east direction at angular distances of  ≳ 0.1′′ (Martí et al. 1998). This confirms that LS 5039 contains an additional extended and diffuse component that cannot be traced with the high-resolution images described in Sects. 4.2 and 5.

6. Model

In this section we try to reproduce the changing radio morphology of LS 5039 by modelling the emission produced by an outflow of electrons accelerated as a consequence of the interaction of the stellar wind and the wind of the putative young pulsar. We are interested in describing the orientation of the PA of the emission up to scales of  ~5–10 mas. The formation and evolution of the outflow is complex, and hydrodynamical and magnetohydrodynamical simulations are needed to describe in detail the source radio emission (see the discussion in Sect. 1). Our aim here is to check if the general scenario of colliding winds is compatible with the main features and the continuous morphological changes found at scales up to 10 mas.

thumbnail Fig. 4

VLBI self-calibrated images of LS 5039 at 5.0 and 8.5 GHz (see Table 4 for details on the projects). The plotting parameters are the same as in Fig. 3. The panel order is in increasing orbital phase from left to right and from top to bottom, and the images have been grouped according to similar morphological features. The project code and orbital phase are quoted on the top part of each panel. The images from project BR127, obtained during the same orbital cycle, have bold axes. Dashed contours are plotted at  −3 times the rms noise of each image and solid contours start at 3 times the rms and scale with 21/2, except for GR021A-B and BP051, which start at 5σ as in the original publications. The rms noise of each image can be found in Table 4.

We assume a mass for the companion star of 33 M and an eccentricity of the system of 0.35 Casares et al. (in prep.). The semimajor axis of the orbit for a 1.4 M neutron star is approximately a = 2.4 × 1012 cm, or 0.05 mas for a distance to the source of 2.9 kpc Casares et al. (in prep.). The semimajor axis does not change significantly for compact object masses below 5 M. Bosch-Ramon et al. (2012) identified regions where strong shocks are produced at distances between 5 and 10a from the massive star, towards the direction of the pulsar. Here we will assume that the extended radio emission is produced by electrons accelerated at a distance of 10 times the semimajor axis of the orbit,  ~0.5 mas at 2.9 kpc, and neglect the high-energy processes taking place at shorter distances. We show a sketch of the scenario in the first row of panels of Fig. 6, where the grey area indicates the flow of particles accelerated at a distance of 10a from the star and are travelling away from it. We consider that the injected electrons follow an energy power law-distribution, . The global radio spectral index of the source is  −0.46 (Martí et al. 1998), and therefore we use an electron index p = 2. For simplicity, we assume a constant injection rate of Qinj = 1035 erg s-1 (see below) along the orbit, for electron energies of 1 < γe < 104. Higher electron energies do not significantly contribute to the radio emission at the considered frequencies. On the other hand, a range for the lower electron energy of  ~1–100 slightly modifies the innermost part of the emission, but does not significantly affect the larger source structure at 2–5 mas. The model is not detailed enough to constrain this value and therefore we assume the minimum energy possible. We note that much higher minimum energies are expected in the primary shock between the stellar and the pulsar wind (Kennel & Coroniti 1984) but we are here considering a secondary shock at larger distances. Zabalza et al. (2011) inferred a spin-down luminosity for the putative pulsar in LS 5039 of Lsd ~ (3−6) × 1036 erg s-1, and therefore our injected energy accounts for  ~2% of this pulsar spin-down luminosity. The energy spectrum was divided into 200 logarithmic intervals for the computations.

For distances from the massive star larger than 10a, the evolution of the electron distribution is determined by the radiative losses (synchrotron and inverse Compton in the Thomson regime) and adiabatic cooling. We assume that the particles travel away from the massive star forming an outflow with conical shape (see Fig. 6). This is conical shape is a coarse estimation, as the ambient pressure and the hydrodynamical instabilities severely modify the flow structure (Bosch-Ramon et al. 2012). We use an opening angle for the cone of 20°, of the same order as the post-shock flow structure seen in Bosch-Ramon et al. (2012, see also Bogovalov et al. 2008), although the final flux density distribution is not very sensitive to this value. We note that this angle is much lower than the asymptotic opening angle of the contact discontinuity, which for this system is about  ~75° assuming a pulsar with spin-down luminosity of 5 × 1036 erg s-1, with a stellar wind with a mass loss rate of 2.5 × 10-7   M yr-1 and wind velocity at infinity of 2.4 × 103 km s-1 (McSwain et al. 2011). We compute the electron energy distribution density n(r,γ), where r is the linear distance from the massive star, taking into account the conservation of the number of particles along the outflow through the continuity equation N(r,γ)dγ = N(r0,γ0)dγ0.

The outflow is bent as the particles travel away from the system because of the orbital motion of the pulsar (see second row of panels of Fig. 6). The emitting region is bent in the opposite direction of the orbital motion of the pulsar and it forms a spiral (see Dubus 2006, for a general description). The flow, similar to a cometary tail, is a projection of the orbit assuming a ballistic motion. However, the outflow, and in general the whole interaction structure, in fact should not follow a ballistic motion because of the effect of the Coriolis force and consequently the resulting spiral should be less open than for simple ballistic motion (Bosch-Ramon & Barkov 2011; Bosch-Ramon et al. 2012). Nevertheless, we approximate the outflow shape assuming a constant advection velocity of 0.04c, which traces a structure similar but more open than the one found in the hydrodynamical simulations by Bosch-Ramon et al. (2012), because we find a better match with the images in this case. We note that those simulations were computed for a slightly different system and using a circular orbit.

We compute the electron density along the flow for r between 10a and 140a, although the farthest part ( ≳ 100a) does not contribute significantly to the emission. We divide the flow into 200 steps in r. The sizes of the steps is determined by the displacement of the pulsar along the orbit and are computed at constant increments of the true anomaly of the pulsar. Therefore, steps close to periastron are shorter, corresponding to time intervals of about 500 s each, whereas steps close to apastron are longer, corresponding to time intervals of about 3000 s each. These intervals are always at least one order of magnitude shorter than the fastest cooling time at each r and γe. For the lower energy electrons, adiabatic cooling dominates along the flow, whereas for more energetic electrons, inverse Compton losses dominate up to distances of  ≲ 4 mas. We compute the emissivity of each section of the flow assuming a magnetic field of  G, which corresponds to B ~ 0.1 G in the acceleration region. We note that this magnetic field is similar to the equipartition magnetic field inferred for non-thermal particles from VLBI observations, which is  ~0.2 G (Paredes et al. 2000).

thumbnail Fig. 5

Self-calibrated image of LS 5039 at 1.7 GHz obtained with the EVN in 2009. The plot parameters are the same as in Fig. 4, although a larger angular scale is shown here.

thumbnail Fig. 6

Sketch of the orbit orientation and the emission region of the model. From left to right the panels show three different scales: the orbit (a), the acceleration region (10a), and the whole flow (100a). The particle acceleration occurs at a distance of 10a from the massive star (blue circle) behind the pulsar (empty green circle) and the flow of emitting particles (grey area) expands, forming a cone while travelling away from the binary system. The star plane is the plane of the sky intersecting the star. The colours in the orbit and the flow indicate regions situated in front of the star plane (red, closer to the observer than the star) and behind the plane (yellow, farther away than the star). The curved arrows indicate the orbital motion of the pulsar. The first row shows the general scenario without orbital motion, whereas the rest of the rows include the orbital motion and different projections of the orbit, indicated in the labels. The last row of panels shows the orbit projected with our best estimate of i and Ω.

The observed outflow geometry strongly depends on the orientation of the orbit on the sky, which is basically determined by two angles: the longitude of the ascending node, Ω, measured from north to east, and the inclination of the orbit, i, which is a rotation with respect to the node axis (the orbital elements are defined in Smart 1930). In Fig. 6 we show different projections of these two angles at orbital phase 0.5. The images in Fig. 4 show a predominant north-east/south-west direction, in a PA of  ~130°. This constant direction for most part of the orbital phases suggests a high inclination of the orbit, as already discussed in Ribó et al. (2008). The orbital elements determine the flow orientation and the final flux density distribution. We searched for the combination of Ω and i that better describes the source structure seen in Fig. 4. We used the electron density and the system geometry to compute the synchrotron emissivity of the electrons and project it on the sky. For every orbital phase in Fig. 4, we computed the sky flux density distribution at the corresponding frequency and convolved it with a beam equal to the corresponding synthesized beam of the image. The model accounts for part of the core component and the extended emission. We produced synthetic images in steps of 5° in i and Ω independently and found that the images in Fig. 4 are better reproduced by the model for an inclination of i ~ 70° and Ω ~130°. Based on visual inspection, a reasonable range for these parameters is i = 60 − 75° and Ω = 120 − 140°. The rotation of the pulsar around the star is counter-clockwise, and therefore i < 90° (see Smart 1930). Soon after the periastron passage the pulsar transits behind the star plane and after apastron it is in front of the star plane. For the inclination obtained the pulsar position is eclipsed during its superior conjunction.

The resulting synthetic images and the projected orbit are shown in Fig. 7. The first contour level for all plots is set to 0.06 mJy beam-1. The model produces total flux densities of around 10–16 mJy, which accounts for the extended flux density recovered with VLBI. However, this model does not account for the emission produced at distances below  ≲ 0.5 mas, so we artificially included an additional 10 mJy point-like source at the position of the binary system to visually help to compare the results with the images in Fig. 4. This component approximately accounts for the peak flux density of the core component of the images. With this ad-hoc component, we recover a total flux density of about 15–18 mJy at 5 GHz and 13–16 mJy at 8.5 GHz, for the considered orbital phases, which are slightly below the measured values.

thumbnail Fig. 7

Synthetic images produced with our model, computed for the corresponding orbital phases of the different VLBI projects and convolved with the corresponding beam shown in Fig. 4. The black line traces the projected position of the flow axis of cooling particles, which is computed for linear distances from the star between 1.6 and 22 AU. The crosses mark the position of the main star, set at the origin. The panel in the top-right corner shows the system orbit, projected with an inclination of 70° and a longitude of the ascending node of 130°.

This model correctly describes the main features in most of the observed images, except for the observations shortly after the periastron passage. For that case (group 1), our model only accounts for the secondary component eastwards, but not for the main extended component towards the west.The rest of the images are well described by this simple model. It produces a bright extended component for orbital phases 0.15–0.40 (group 2) and continuously develops extended emission towards the north-west between phases 0.5 and 0.6 (group 3), while preserving the bipolar structure. The extended emission is then dominated by the north-west component until the next periastron passage. We note that this model predicts bipolar extended emission at certain orbital phases, although the outflow is not bipolar itself. In particular, the projected flow at phase 0.5, shown in the bottom-right panel of Fig. 6, accounts for the bipolar flux density distribution seen for instance in the image of project GR021A (see Figs. 4 and 7). At the same time, this approach justifies that the bipolar emission is not completely symmetric and that the two components show different PA (Ribó et al. 2008).

The differences between the model (Fig. 7) and the real data (Fig. 4), in particular the missing component in the images from group 1, can be explained by the oversimplification of the flow structure and the constant physical conditions at all orbital phases. For instance, we assumed that the acceleration region was at a constant distance from the massive star and that the energy injection was constant, although it is expected that these parameters depend on the orbital phase in an eccentric binary system. Therefore, additional electron populations, as well as more realistic time-dependent physical conditions should be included to produce this additional component. The model should also include considerations on the relativistic Doppler boosting possibly affecting the synchrotron emission (Dubus et al. 2010). However, including additional particle populations and assuming time-dependent physical conditions would require the addition of several free parameters, which is beyond the scope of this paper. Therefore, we do not try to include a complete flux density analysis or to predict flux density distributions. Also, predictions of displacements of the peak of the emission, such as the hint found in Sect. 4.1, cannot be obtained by this model, as these are produced at smaller distances. Despite these caveats, this model traces, except for the phases close to the periastron passage, all the main features of the VLBI images and the morphological variability measured in multifrequency and multiepoch observations, even when assuming simple and constant physical conditions.

7. Discussion and conclusions

VLBA observations at 5 GHz during five consecutive days show that LS 5039 displays orbital morphological variability, showing one-sided and bipolar structures, but recovering the same morphology when observing at the same orbital phase. The PA of the extended emission with respect to the core changes significantly on timescales of one day. Although the total flux density remains approximately constant during the orbital cycle, the flux density of the different components varies. Remarkably, the PA of the extended emission changes by  ~160° between phases 0.03 and 0.29 (24 h). The sensitivity of the observations is not sufficient to study changes on shorter timescales, within a few hours. We measured a displacement of the peak of the emission between orbital phases 0.29 and 0.55 of 2.3 ± 0.6 mas in PA 154°, nearly in the same direction of the extended emission. However, we consider this measurement to be only a hint of peak displacement because, given the limitations of the astrometry of these observations, we would require a 5-σ significance to consider it a significant displacement. Dubus (2006) predicted a continuous peak displacement forming an ellipse of a few mas along the orbit. The peak of the emission of the gamma-ray binary LS I +61 303 at 8.4 GHz traces an ellipse of semimajor axis 1.5 AU, about four times the binary semimajor axis (Dhawan et al. 2006). Unfortunately, our astrometric accuracy is not sufficient to measure these displacements. The source morphology during two consecutive periastron passages (run A and E) is very similar, showing that the changes within the same orbital cycle are periodic. An interesting future project would be to obtain accurate astrometry to explore the displacements of the peak of the emission, in particular at different frequencies, and also to measure variations of the source structure on timescales of a few hours.

We analysed the available VLBI observations of LS 5039 in a consistent way. Images from data obtained between 1999 and 2009 at frequencies between 1.7 and 8.5 GHz show that the morphology at similar orbital phases is similar. When ordered in increasing phase, the images show an approximately continuous change. The observational conclusion of our VLBI analysis is that the morphological changes of LS 5039 are periodic and show orbital modulation stable over several years. The gamma-ray binary LS I +61 303 also shows a similar behaviour (Dhawan et al. 2006). Therefore, all gamma-ray binaries are expected to display periodic orbital modulation of their VLBI structure.

A simple model of an expanding outflow of relativistic electrons accelerated at a distance of  ~10a and travelling away from the system can account for the main changes seen in the extended emission of LS 5039. The changes are explained for an inclination of the orbit i ~ 70° and a longitude of the ascending node of Ω ~ 130°. In the last row of panels in Fig. 6 we show the orientation of the orbit on the sky projected with these parameters. The inclination of the orbit has deep implications for the physical properties of the binary system, because the orbital parameters restrict the mass of the compact object within the system. In particular, assuming a mass function of the binary system of 0.0045(6) M Casares et al. (in prep.), a mass of the star of 33 M implies that the mass of the compact object is 1.8–2.0 M, for an inclination between 75 and 60°. However, the mass of the main star is barely constrained between 20 and 50 M, and consequently the inclinations inferred from the model are compatible with masses between 1.3–1.5 and 2.4–2.7 M, for the lower and the higher stellar mass limits, respectively. The uncertainty in the mass function adds an additional uncertainty of  ~0.1 M. The values obtained are in any case below 3 M, and therefore are compatible with the presence of a neutron star in the system. However, we note that this was one hypothesis of the model and therefore cannot be a proof. After applying a simple model, we conclude that the presence of a young non-accreting pulsar in LS 5039 agrees with the observed source structure at different orbital phases and produces a coherent scenario. However, this model does not account for all emission produced at shorter distances and also fails to explain the main component of the extended emission seen shortly after the periastron passage. This is possibly due to unaccounted variability in the injection rate, geometry and localisation of the acceleration zone, and flow velocity of the electrons.

For gamma-ray binaries, the nature of the compact object, either accreting or not, is a fundamental ingredient for any model aimed at understanding these peculiar systems. The most straightforward indication would be the detection of pulsations. However, direct detection may not be feasible in binary systems with close orbits and powerful massive companions because of free-free absorption (Dubus 2006). In addition, the mass function of the known gamma-ray binaries does not allow to distinguish between black holes and neutron stars. Therefore, the debate is still open regarding the nature of the compact object of gamma-ray binaries. We have presented high-resolution radio images of LS 5039 with which the scenarios can be tested. A first approach has shown that the source radio morphology can be explained with the presence of a young non-accreting pulsar. PSR B1259 − 63, the only gamma-ray binary with a confirmed pulsar, shows variable extended radio structures (Moldón et al. 2011a), and LS I +61 303 shows variable and periodic morphological changes (Dhawan et al. 2006). The common features among these three systems suggest that the known gamma-ray binaries contain young non-accreting pulsars.


1

The NRAO Astronomical Image Processing System. http://www.aips.nrao.edu/

3

The Crustal Dynamics Data Information System http://cddis.nasa.gov/

Acknowledgments

We thank V. Bosch-Ramon for very useful discussions on theoretical modelling. We are grateful to J. Casares for allowing us to mention several new parameters of LS 5039 prior to publication. We also thank V. Zabalza and K. Iwasawa for helpful comments and discussions. We thank the anonymous referee for very useful comments. The Very Long Baseline Array is operated by the USA National Radio Astronomy Observatory, which is a facility of the USA National Science Foundation operated under co-operative agreement by Associated Universities, Inc. The European VLBI Network (http://www.evlbi.org/) is a joint facility of European, Chinese, South African, and other radio astronomy institutes funded by their national research councils. This research has made use of SAO/NASA’s Astrophysics Data System. We acknowledge support by the Spanish Ministerio de Ciencia e Innovación (MICINN) under grants AYA2010-21782-C03-01 and FPA2010-22056-C06-02. J.M. acknowledges support by MICINN under grant BES-2008-004564. M.R. acknowledges financial support from MICINN and European Social Funds through a Ramón y Cajal fellowship. J.M.P. acknowledges financial support from ICREA Academia.

References

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009a, ApJ, 701, L123 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009b, ApJ, 706, L56 [NASA ADS] [CrossRef] [Google Scholar]
  3. Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005, Science, 309, 746 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  4. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 460, 743 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Aharonian, F. A., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, A&A, 469, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Aharonian, F., Akhperjanian, A. G., Anton, G., et al. 2009, A&A, 507, 389 [Google Scholar]
  7. Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 665, L51 [NASA ADS] [CrossRef] [Google Scholar]
  8. Albert, J., Aliu, E., Anderhub, H., et al. 2009, ApJ, 693, 303 [NASA ADS] [CrossRef] [Google Scholar]
  9. Aragona, C., McSwain, M. V., Grundstrom, E. D., et al. 2009, ApJ, 698, 514 [NASA ADS] [CrossRef] [Google Scholar]
  10. Barkov, M. V., & Khangulyan, D. V. 2012, MNRAS, 421, 1351 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bhattacharyya, S., Godambe, S., Bhatt, N., Mitra, A., & Choudhury, M. 2012, MNRAS, 421, L1 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bogovalov, S. V., Khangulyan, D. V., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2008, MNRAS, 387, 63 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bogovalov, S. V., Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2012, MNRAS, 419, 3426 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bongiorno, S. D., Falcone, A. D., Stroh, M., et al. 2011, ApJ, 737, L11 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bosch-Ramon, V., & Barkov, M. V. 2011, A&A, 535, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bosch-Ramon, V., & Khangulyan, D. 2009, Int. J. Mod. Phys. D, 18, 347 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bosch-Ramon, V., & Khangulyan, D. 2011, PASJ, 63, 1023 [NASA ADS] [Google Scholar]
  18. Bosch-Ramon, V., Paredes, J. M., Ribó, M., et al. 2005, ApJ, 628, 388 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bosch-Ramon, V., Romero, G. E., & Paredes, J. M. 2006, A&A, 447, 263 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Bosch-Ramon, V., Motch, C., Ribó, M., et al. 2007, A&A, 473, 545 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Bosch-Ramon, V., Khangulyan, D., & Aharonian, F. A. 2008, A&A, 482, 397 [Google Scholar]
  22. Bosch-Ramon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A, 544, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Casares, J., Ribó, M., Ribas, I., et al. 2005, MNRAS, 364, 899 [NASA ADS] [CrossRef] [Google Scholar]
  24. Casares, J., Ribó, M., Ribas, I., et al. 2012, MNRAS, 421, 1103 [NASA ADS] [CrossRef] [Google Scholar]
  25. Cerutti, B., Dubus, G., & Henri, G. 2008, A&A, 488, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Cerutti, B., Malzac, J., Dubus, G., & Henri, G. 2010, A&A, 519, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Dhawan, V., Mioduszewski, A., & Rupen, M. 2006, in VI Microquasar Workshop: Microquasars and Beyond [Google Scholar]
  28. Dubus, G. 2006, A&A, 456, 801 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Dubus, G., Cerutti, B., & Henri, G. 2008, A&A, 477, 691 [Google Scholar]
  30. Dubus, G., Cerutti, B., & Henri, G. 2010, A&A, 516, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Durant, M., Kargaltsev, O., Pavlov, G. G., Chang, C., & Garmire, G. P. 2011, ApJ, 735, 58 [NASA ADS] [CrossRef] [Google Scholar]
  32. Fermi LAT Collaboration, Abdo, A. A., Ackermann, M., et al. 2009, Science, 326, 1512 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  33. Fermi LAT Collaboration, Ackermann, M., Ajello, M., et al. 2012, Science, 335, 189 [NASA ADS] [CrossRef] [Google Scholar]
  34. Fey, A. L., Spangler, S. R., & Cordes, J. M. 1991, ApJ, 372, 132 [NASA ADS] [CrossRef] [Google Scholar]
  35. Fomalont, E. B., Petrov, L., MacMillan, D. S., Gordon, D., & Ma, C. 2003, AJ, 126, 2562 [NASA ADS] [CrossRef] [Google Scholar]
  36. Godambe, S., Bhattacharyya, S., Bhatt, N., & Choudhury, M. 2008, MNRAS, 390, L43 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kennel, C. F., & Coroniti, F. V. 1984, ApJ, 283, 710 [NASA ADS] [CrossRef] [Google Scholar]
  38. Khangulyan, D., Hnatic, S., Aharonian, F., & Bogovalov, S. 2007, MNRAS, 380, 320 [NASA ADS] [CrossRef] [Google Scholar]
  39. Khangulyan, D., Aharonian, F., & Bosch-Ramon, V. 2008, MNRAS, 383, 467 [NASA ADS] [CrossRef] [Google Scholar]
  40. Kishishita, T., Tanaka, T., Uchiyama, Y., & Takahashi, T. 2009, ApJ, 697, L1 [NASA ADS] [CrossRef] [Google Scholar]
  41. Lamberts, A., Dubus, G., Lesur, G., & Fromang, S. 2012, A&A, 546, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Maier, G., for the VERITAS Collaboration, Skilton, J., & for the HESS Collaboration 2011, Proc. 32nd International Cosmic Ray Conference (ICRC2011), Vol. 7, OG2.1-2.2, 78 [Google Scholar]
  43. Maraschi, L., & Treves, A. 1981, MNRAS, 194, 1 [Google Scholar]
  44. Martí, J., Paredes, J. M., & Ribó, M. 1998, A&A, 338, L71 [NASA ADS] [Google Scholar]
  45. Martocchia, A., Motch, C., & Negueruela, I. 2005, A&A, 430, 245 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Massi, M., Ros, E., & Zimmermann, L. 2012, A&A, 540, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. McSwain, M. V., Ray, P. S., Ransom, S. M., et al. 2011, ApJ, 738, 105 [NASA ADS] [CrossRef] [Google Scholar]
  48. Moldón, J., Johnston, S., Ribó, M., Paredes, J. M., & Deller, A. T. 2011a, ApJ, 732, L10 [NASA ADS] [CrossRef] [Google Scholar]
  49. Moldón, J., Ribó, M., & Paredes, J. M. 2011b, A&A, 533, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Moldón, J., Ribó, M., Paredes, J. M., et al. 2012, A&A, 543, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Okazaki, A. T., Romero, G. E., & Owocki, S. P. 2008, in Proc. 7th INTEGRAL Workshop [Google Scholar]
  52. Okazaki, A. T., Nagataki, S., Naito, T., et al. 2011, PASJ, 63, 893 [NASA ADS] [Google Scholar]
  53. Pandey, M., Rao, A. P., Ishwara-Chandra, C. H., Durouchoux, P., & Manchanda, R. K. 2007, A&A, 463, 567 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Paredes, J. M. 2011 [arXiv:1101.4843] [Google Scholar]
  55. Paredes, J. M., Martí, J., Ribó, M., & Massi, M. 2000, Science, 288, 2340 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  56. Paredes, J. M., Ribó, M., Ros, E., Martí, J., & Massi, M. 2002, A&A, 393, L99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Paredes, J. M., Bosch-Ramon, V., & Romero, G. E. 2006, A&A, 451, 259 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Paredes, J. M., Ribó, M., Bosch-Ramon, V., et al. 2007, ApJ, 664, L39 [NASA ADS] [CrossRef] [Google Scholar]
  59. Pavlov, G. G., Chang, C., & Kargaltsev, O. 2011, ApJ, 730, 2 [NASA ADS] [CrossRef] [Google Scholar]
  60. Pradel, N., Charlot, P., & Lestrade, J.-F. 2006, A&A, 452, 1099 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Rea, N., Torres, D. F., Caliandro, G. A., et al. 2011, MNRAS, 416, 1514 [NASA ADS] [CrossRef] [Google Scholar]
  62. Reig, P., Ribó, M., Paredes, J. M., & Martí, J. 2003, A&A, 405, 285 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Ribó, M., Reig, P., Martí, J., & Paredes, J. M. 1999, A&A, 347, 518 [NASA ADS] [Google Scholar]
  64. Ribó, M., Paredes, J. M., Romero, G. E., et al. 2002, A&A, 384, 954 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Ribó, M., Paredes, J. M., Moldón, J., Martí, J., & Massi, M. 2008, A&A, 481, 17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Romero, G. E., Okazaki, A. T., Orellana, M., & Owocki, S. P. 2007, A&A, 474, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Sarty, G. E., Szalai, T., Kiss, L. L., et al. 2011, MNRAS, 411, 1293 [NASA ADS] [CrossRef] [Google Scholar]
  68. Shepherd, M. C. 1997, in Astron. Data Anal. Software Syst. VI, eds. G. Hunt, & H. Payne, 77, ASP Conf. Ser., 125 [Google Scholar]
  69. Sierpowska-Bartosik, A., & Torres, D. F. 2007, ApJ, 671, L145 [NASA ADS] [CrossRef] [Google Scholar]
  70. Smart, W. M. 1930, MNRAS, 90, 534 [NASA ADS] [Google Scholar]
  71. Szostek, A., & Dubus, G. 2011, MNRAS, 411, 193 [NASA ADS] [CrossRef] [Google Scholar]
  72. Takahashi, T., Kishishita, T., Uchiyama, Y., et al. 2009, ApJ, 697, 592 [NASA ADS] [CrossRef] [Google Scholar]
  73. Takata, J., Okazaki, A. T., Nagataki, S., et al. 2012, ApJ, 750, 70 [NASA ADS] [CrossRef] [Google Scholar]
  74. Tavani, M., & Arons, J. 1997, ApJ, 477, 439 [NASA ADS] [CrossRef] [Google Scholar]
  75. Tavani, M., Bulgarelli, A., Piano, G., et al. 2009, Nature, 462, 620 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  76. Thompson, A. R., Moran, J. M., & Swenson, G. W. 1986, Interferometry and synthesis in radio astronomy (New York: Wiley-Interscience) [Google Scholar]
  77. Zabalza, V., Bosch-Ramon, V., & Paredes, J. M. 2011, ApJ, 743, 7 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Log of VLBA observations of project BR127, conducted at 5 GHz.

Table 2

Relative astrometry of the five runs of project BR127, also plotted in Fig. 2.

Table 3

Gaussian components fitted to the images shown in Fig. 3.

Table 4

Log of the VLBI observations of LS 5039.

All Figures

thumbnail Fig. 1

Distribution on the sky of LS 5039, the phase calibrator J1825 − 1718, and the astrometric check source J1818 − 1108. The dashed line indicates the Galactic Plane. The arrows indicate the two phase calibrations obtained: direct, using J1825 − 1718 as phase reference; inverse, using LS 5039 as phase reference.

In the text
thumbnail Fig. 2

Relative offsets of the peak of the emission of LS 5039 during five consecutive days with respect to the average position. For straightforward comparison, we plot the opposite value of the inverse astrometry. The run label is displayed close to each direct measure. The size of the binary system is indicated by the small orbit in the corner, which is plotted face-on and with an arbitrary orientation. The displacements and uncertainties are shown in Table 2.

In the text
thumbnail Fig. 3

VLBA images of LS 5039 at 5 GHz from project BR127, obtained during five consecutive days in July 2007. Each column, labelled with the run name and orbital phase, corresponds to one epoch (see Table 1). For each epoch, the first and second rows correspond to the self-calibrated images obtained with a natural and uniform weighting scheme, respectively. The third row shows the uv components fitted to the data, convolved with a Gaussian circular beam with an area equal to the average synthesized beam of all images. The lines show the direction of the main extended component. The restoring beams are plotted in the bottom-left corner of each panel. Dashed contours are plotted at  − 3 times the rms noise of each image and solid contours start at 3 times the rms and scale with 21/2. The Gaussian components fitted to the images are listed in Table 3.

In the text
thumbnail Fig. 4

VLBI self-calibrated images of LS 5039 at 5.0 and 8.5 GHz (see Table 4 for details on the projects). The plotting parameters are the same as in Fig. 3. The panel order is in increasing orbital phase from left to right and from top to bottom, and the images have been grouped according to similar morphological features. The project code and orbital phase are quoted on the top part of each panel. The images from project BR127, obtained during the same orbital cycle, have bold axes. Dashed contours are plotted at  −3 times the rms noise of each image and solid contours start at 3 times the rms and scale with 21/2, except for GR021A-B and BP051, which start at 5σ as in the original publications. The rms noise of each image can be found in Table 4.

In the text
thumbnail Fig. 5

Self-calibrated image of LS 5039 at 1.7 GHz obtained with the EVN in 2009. The plot parameters are the same as in Fig. 4, although a larger angular scale is shown here.

In the text
thumbnail Fig. 6

Sketch of the orbit orientation and the emission region of the model. From left to right the panels show three different scales: the orbit (a), the acceleration region (10a), and the whole flow (100a). The particle acceleration occurs at a distance of 10a from the massive star (blue circle) behind the pulsar (empty green circle) and the flow of emitting particles (grey area) expands, forming a cone while travelling away from the binary system. The star plane is the plane of the sky intersecting the star. The colours in the orbit and the flow indicate regions situated in front of the star plane (red, closer to the observer than the star) and behind the plane (yellow, farther away than the star). The curved arrows indicate the orbital motion of the pulsar. The first row shows the general scenario without orbital motion, whereas the rest of the rows include the orbital motion and different projections of the orbit, indicated in the labels. The last row of panels shows the orbit projected with our best estimate of i and Ω.

In the text
thumbnail Fig. 7

Synthetic images produced with our model, computed for the corresponding orbital phases of the different VLBI projects and convolved with the corresponding beam shown in Fig. 4. The black line traces the projected position of the flow axis of cooling particles, which is computed for linear distances from the star between 1.6 and 22 AU. The crosses mark the position of the main star, set at the origin. The panel in the top-right corner shows the system orbit, projected with an inclination of 70° and a longitude of the ascending node of 130°.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.