Multi-wavelength (radio, X-ray and gamma-ray) observations of the gamma-ray binary LS I +61 303

We present the results of the first multiwavelength observing campaign on the high-mass X-ray binary LS I +61 303 comprising observations at the TeV regime with the MAGIC telescope, along with X-ray observations with Chandra, and radio interferometric observations with the MERLIN, EVN and VLBA arrays, in October and November 2006. From our MERLIN observations, we can exclude the existence of large scale (~100 mas) persistent radio-jets. Our 5.0 GHz VLBA observations display morphological similarities to previous 8.4 GHz VLBA observations carried out at the same orbital phase, suggesting a high level of periodicity and stability of the processes behind the radio emission. This makes it unlikely that variability of the radio emission is due to the interaction of an outflow with variable wind clumps. If the radio emission is produced by a milliarcsecond scale jet, it should also show a stable, periodic behavior. It is then difficult to reconcile the absence of a large scale jet (~100 mas) in our observations with the evidence of a persistent relativistic jet reported previously. We find a possible hint of temporal correlation between the X-ray and TeV emissions and evidence for radio/TeV non-correlation, which points to the existence of one population of particles producing the radio emission and a different one producing the X-ray and TeV emissions. Finally, we present a quasi-simultaneous energy spectrum including radio, X-ray and TeV bands.

1. INTRODUCTION LS I +61 303 is a high-mass X-ray binary consisting of a low-mass [M ∼ (1 − 4) M ⊙ ] compact object orbiting around an early type B0 Ve star along an eccentric (e = 0.7) orbit (Casares et al. 2005, and references therein). The modulation of both, the radio (Gregory & Taylor 1978;Gregory 2002) and Xray (Taylor et al. 1996;Paredes et al. 1997) emissions, display a period of P orb = 26.496 d, attributed to the orbital motion. The position of the maximum of the radio emission along the orbit, as well as its intensity, are modulated with a superorbital period of P sup = 1667 ± 8 days (Gregory 2002). LS I +61 303 is positionally coincident with an EGRET γ-ray source (Kniffen et al. 1997). Moreover, variable emission at TeV energies has been recently detected with the MAGIC telescope (Albert et al. 2006a). These authors found that the peak flux at TeV energies occurs at orbital phase φ orb ≈ 0.65, while no high-energy emission is detected around periastron passage (φ orb ≈ 0.23). From ∼50 mas resolution radio images of LS I +61 303 obtained with MERLIN, Massi et al. (2004) suggested the existence of a precessing relativistic (β ≈ 0.6) jet up to angular scales of ∼0.1 arcsec, which led them to interpret LS I +61 303 within the framework of the microquasar scenario (Bosch-Ramon et al. 2006;Romero et al. 2007). However, recent VLBA imaging obtained by Dhawan et al. (2006) over a full orbit of LS I +61 303 has shown the radio emission to come from angular scales smaller than about 7 mas (projected size 14 AU at an assumed distance of 2 kpc). This radio emission appeared cometary-like, interpreted to be pointing away from the high mass star within a particular scenario, and no relativistic motion, nor halos, nor larger-scale structures were detected at any phase of the orbit. Based on these findings, Dhawan et al. (2006) concluded that the radio and TeV emissions from LS I +61 303 are originated by the interaction of the wind of a young pulsar with that of the stellar companion (Maraschi & Treves 1981;Dubus 2006).
In this article, we present and discuss the results of a multi-wavelength campaign including radio, X-ray, and TeV γ-ray observations on LS I +61 303 , aimed at shedding light on the physical processes going on in the system, as well as yielding useful input for a detailed, timedependent modeling of this relevant system. For this, we study the correlations between the different detected emissions in terms of their morphological, temporal and spectral features, both in the intra-day and day-to-day timescales. Table 1 and Figure 1 summarize our observations of LS I +61 303 , carried out during October and November 2006. In particular, we set up a simultaneous, multiwavelength campaign on LS I +61 303 for the 25 and 26 October 2006 using the MERLIN, EVN, and VLBA interferometers (radio), and the TNG (infra-red), Chandra (X-rays), and MAGIC TeV γ-rays) telescopes. Unfortunately, bad weather conditions did not allow us to get any infra-red data, nor did we get useful data from the MAGIC telescope coincident in time with simultaneously scheduled Chandra observations. The period 16-20 The orbital and super-orbital phases are computed using MJD 0 = 43366.275, P orb = 26.4960 (Gregory 2002). The superorbital phase is φsup = 0.4 for all the observation period. b The average measured fluxes are in mJy for radio observations (MERLIN, EVN, VLBA), 10 −11 ergs cm −2 s −1 for the Chandra X-ray observations (absorbed power-law), and 10 −11 photons cm −2 s −1 for the MAGIC VHE γ-ray observations. Upper limits are given at the 95% confidence level, following the prescription by Rolke et al. (2005). The data were taken at 5.0 GHz for the radio observations (except for the MERLIN data in November 2006, which were taken at 6.0 GHz), at 0.5-10 keV by Chandra, and at 0.3-5.0 TeV by MAGIC.

November included only MAGIC and MERLIN observations.
The range of observed orbital phase in the October campaign was φ orb = 0.6-0.7, for which the TeV maximum was detected in the past (Albert et al. 2006a), and φ orb = 0.44-0.57 in the November campaign. The super-orbital phase for both campaigns was φ sup = 0.4.

Radio Observations
We observed LS I +61 303 with several radio interferometric arrays during 25-26 October 2006, including the Multi-Element Radio Linked Interferometer (MERLIN) in the UK, the European VLBI Network (EVN), and the Very Long Baseline Array (VLBA) in the USA, all of which observed at 5.0 GHz. Additionally, we also monitored the flux and large angular scale structural and flux density variations of LS I +61 303 using MERLIN at 5.0 GHz on 27 and 28 October 2006, and at 6.0 GHz for the period from 16 to 20 November 2006.
The MERLIN array included six antennas (Defford, Cambridge, Knockin, Darnhall, Jodrell Bank (MkII), and Tabley) for most of our observations, yielding synthesized beams of about 50 to 70 mas (corresponding to projected linear sizes of 100 to 140 AU). The EVN observations on 26 October 2006, which were the first ever carried out for LS I +61 303 using the e-VLBI technique (Szomoru 2006), included five antennas spread over Europe (Cambridge, Jodrell Bank MkII, Medicina, Torun, and the Westerbork phased array), yielding a synthesized beam of about 7 mas (or 14 projected AU). The data were directly streamed to the correlator at JIVE (Joint Institute for VLBI in Europe, Dwingeloo, The Netherlands) through the Internet, and this was one of the first science observations that achieved 256 Mbps sustained data rate with the e-VLBI technique. The VLBA observations on 25 and 26 October 2006 included ten 25-m antennas spread over the USA, yielding a beam size of 4.6×2.1 mas, corresponding to a projected linear resolution of ∼9.2 AU and ∼4.2 AU in right ascension and declination, respectively.
All three interferometric arrays observed LS I +61 303 in phase-referenced mode.
LS I +61 303 and the bright (S 5GHz ≈ 600mJy), nearby International Celestial Reference Frame (ICRF) source J0244+6228 were alternately observed through each observing run, the phases of J0244+6228 being transfered to the position of LS I +61 303 in the posterior data analysis. J0244+6228 also served as amplitude calibrator for our observations. We performed standard calibration and data reduction within the NRAO Astronomical Imaging Package System (AIPS; Diamond 1995). We also used standard hybrid mapping techniques within AIPS to obtain the flux densities shown in Figure 1 and Table 1, and the radio images shown in Figures 2 and 3

X-ray Observations
We obtained Chandra X-ray observations of LS I +61 303 on 25 October 2006 through the Director Discretionary Time program (Chandra ObsId 8273). The observations were carried out using ACIS-I for a total exposure time of 20.0 ks. At the time of the observation, LS I +61 303 was expected to be in a high-state and its high X-ray brightness could have resulted in an excessively high count rate producing appreciable pile-up in the observations. In order to minimize pile-up effects, LS I +61 303 was offset by 8 ′ from the ACIS-I aim-point, thus smearing its image and reducing the count rate at the source peak. We further used 1/4 sub-array, reducing the exposure time of individual frames from the nominal exposure time of 3.2 s down to 0.8 s.
Data reduction of the X-ray observations of LS I +61 303 have been performed using the Chandra X-ray Center software CIAO V3.4. The data reduction included the application of standard filters and rejection of events with bad grades and those originating from bad pixels. The background count rate is consistent with the quiescent background (Markevitch 2001), and no time intervals of enhanced background needed to be removed. The processed observations have a useful exposure time of 19.1 ks. Due to the high count rate of LS I +61 303 out-of-time events were not negligible and produced a noticeable streak along the columns. These have been corrected using the CIAO task "acisreadcorr". Finally, notice that as of CALDB v3.1.0 (June 2005), the CIAO task "acis process events" is routinely used by the level 1 processing pipeline to mitigate the charge transfer inefficiency (CTI) that affects the Chandra ACIS-I front illuminated chips. Therefore, no further CTI correction was applied. Lightcurves and spectra of LS I +61 303 were obtained using standard CIAO tasks and analyzed using HEASARC FTOOLS and XSPEC V11.2.0 routines (Arnaud 1996).

VHE γ-ray Observations
Observations in the very high energy (VHE) γ-ray band (E γ > 100 GeV) were scheduled with the MAGIC telescope on 26-28 October and 16-19 November 2006 (see Table 1). These observations are part of an extensive, standalone observing campaign carried out between September and December 2006 (Albert et al. 2008, in preparation). Bad weather conditions at the Observatorio del Roque de los Muchachos, however, prevented us from obtaining useful data on 26 and 28 October 2006.
The observations were carried out in the false-source track (wobble) mode (Fomin et al. 1994), with two directions at 24 ′ distance and opposite sides of the source direction, which allows for a reliable estimation of the background with no need of extra observation time.
The data were analyzed using the standard MAGIC calibration and analysis software (Albert et al. 2006b;Gaug et al. 2005). Data runs with anomalous event rates were discarded for further analysis. Hillas variables (Hillas 1985) were combined into an adimensional γ/hadron discriminator (hadronness) and an energy estimator by means of the Random Forest classification algorithm, which takes into account the correlation between the different Hillas variables (Breiman 2001; Bock et al. 2005). The incoming direction of the primary γ-ray events was estimated using the DISP method, suited for observations with a single IACT (Fomin et al. 1994;Domingo-Santamaría et al. 2005).

RESULTS
In this section we present the results obtained from the different multi-wavelength observations we have performed, and put them into the context of the past measurements in the different bands. They are summarized in Table 1 and Figure 1.

Radio Results
The total radio flux density obtained with MERLIN shows a decline between October 26 (φ orb = 0.66) and 28 (φ orb = 0.70), from ∼ 35 to ∼ 30 mJy, and a peak on 18 November (φ orb = 0.50) at ∼ 80 mJy. At the superorbital phase of the observations (φ sup = 0.4) the radio source is in the weak state. The predicted flux of the radio flux is between 50-100 mJy at the orbital phases φ orb ∼ 0.9. We therefore measure a flux compatible with the predicted one, although at a much earlier phase value. However, for the weak state, secondary peaks of comparable flux show up at other orbital phase values (see, e.g., Dhawan et al. 2006).
In Figure 2 we show the radio images corresponding to 25-26 of October, where all three arrays observed simultaneously LS I +61 303 . In Figure 3 we show the MERLIN images obtained on 27 and 28 October 2006 at 5.0 GHz, and on 16 through 20 November at 6.0 GHz. Our MERLIN observations show no evidence for large angular scale structures, contrary to what had been previously claimed (Massi et al. 2004), hence excluding the existence of persistent jets at these scales. The structure of LS I +61 303 genuinely follows the beam shape of the MERLIN array at all epochs. Model fits to the u−v plane of the MERLIN data show that the projected intrinsic size of the radio emitting region of LS I +61 303 is, at any observing epoch, no larger than ∼6 mas (∼12 AU), as confirmed by our higher resolution images obtained with the EVN and the VLBA.
Our 5 Fig. 1.-Radio, X-ray and VHE γ-ray light-curves obtained with VLBA, EVN, MERLIN (left-hand scale), Chandra (right-hand scale in 10 −11 ergs cm −2 s −1 ) and MAGIC (right-hand scale in 10 −11 photons cm −2 s −1 ) during the two observing periods (October and November 2006). The horizontal error bars show the time spanned by the different observations. The shaded area marks the range of X-ray flux values previously reported (see Table 3). The upper axis shows the orbital phase using the ephemeris from Gregory (2002).  Table  1).  Table 1). Note that the shape of the radio brightness distribution of LS I +61 303 genuinely follows the synthesized beam for each epoch (upper left inset in each panel). radio emission to the east-southeast (see Figure 2). The general aspect of the radio emission recalls the structure found by Dhawan et al. (2006). In fact, panel (a) of Figure 2 shows both our 5 GHz VLBA image on 25 October 2006 (φ orb = 0.62, contours) overlaid on top of the 8.4 GHz VLBA data on 2 February 2006 (φ orb = 0.62, grey scale), kindly provided by V. Dhawan in advance of publication.

MilliARC
The projected distance from the peak of the radio brightness distribution to the edge of the extended region is, as imaged with the VLBA at 5.0 GHz on 26 October 2006, of about ∼8.2 mas (16.4 AU). If the maximum time to fill in this region with radio (synchrotron) emitting particles is less that the time spanned by our two consecutive VLBA observations (1 day), the implied outflow velocity is at least of 21500 km/s, or v > ∼ 0.09 c. The coordinates of the peak of radio brightness distribution of LS I +61 303 , as obtained from our VLBA observations (using task IMFIT), were RA = 02 h 40 m 31.6638849 s , DEC = 61 • 13 ′ 45.592235 ′′ on October 25 and RA = 02 h 40 m 31.6638756 s , DEC = 61 • 13 ′ 45.592496 ′′ on October 26, with and estimated accuracy of 12 µas and 7 µas in RA and DEC, respectively. This shift of the brightness peak corresponds to a dayto-day projected speed of 904±60 km s −1 , in excellent agreement with the typical value of ∼1000 km s −1 found by Dhawan et al. (2006) along a complete orbital cycle.
The peak flux density varied significantly between our two VLBA observing runs, decreasing from 24.3 mJy/b on October 25 to 15.4 mJy/b on October 26. The total 5 GHz flux density varied from 39.5±0.2 mJy (radio luminosity L R = (9.5 ± 0.1) × 10 29 ergs s −1 ) to 32.7±0.2 mJy (L R = (7.8±0.1)×10 29 ergs s −1 ). We therefore suggest that the drop seen in the total radio flux density of LS I +61 303 between our two consecutive VLBA observations is directly related to the change in the peak flux density, which hints towards a physical link with structural variations in the innermost ( < ∼ 3 mas projected radius, or 6 AU at a distance of 2 kpc) region of the source, which our observations cannot resolve.

X-ray Results
LS I +61 303 is clearly detected in our Chandra observation, with a background-subtracted averaged count rate of 1.067±0.008 cnts s −1 in the energy band 0.5-10.0 keV. This count rate would have indeed resulted in considerable pile-up for a standard observational setup. The X-ray spectrum (Figure 4) shows a broad peak at ∼1.6 keV and a hard energy tail that extends up to 9 keV. Assuming a foreground interstellar absorption with solar abundances and the absorption cross-sections of Balucinska-Church & McCammon (1992), the X-ray spectrum of LS I +61 303 can be reasonably well fitted by either an absorbed power law, or an absorbed Bremsstrahlung model (Fig. 4). The best fit parameters, goodness of the fit (reduced χ 2 ), and implied source flux and luminosity in the energy range of 0.5-10.0 keV for these two models are listed in Table 2. In both cases, the unabsorbed X-ray luminosity of LS I +61 303 in the energy range of 0.5-10.0 keV is ∼ 10 34 ergs s −1 .
We have also investigated the short-term variability of LS I +61 303 X-ray flux and hardness ratio (defined as the ratio of the count rate between 2.0-10.0 keV over that  Table 3 for details). The lower panel shows the relative residuals of the fits (∆I) in terms of the bin standard deviation (σ) for both the absorbed power law (red dots) and Bremsstrahlung (blue dots) models. The errors bars in both the spectra and residual plots are 1-σ. The insets show the χ 2 plots as a function of the power law index, Γ, and N H (upper inset), and of the plasma temperature, kT , and N H (lower inset) of the spectral fits to the Chandra ACIS-I spectrum of LS I +61 303 using absorbed power law and Bremsstrahlung models, respectively. The contours represent 68%, 90%, and 99% confidence levels. 1.84±0.05 L 0.5−10.0 keV 1.18±0.11 1.11±0.04 a Best fit parameters to the X-ray spectrum for LS I +61 303 obtained using absorbed power-law and Bremsstrahlung models. The flux between 0.5 and 10.0 keV is expressed in 10 −11 ergs cm −2 s −1 and the luminosity in 10 34 ergs s −1 and for an assumed distance of 2 kpc.  1.87 ± 0.02 a The first column shows an ID composed of a letter standing for the X-ray satellite (A = ASCA (Leahy et al. 1997), R = RXTE (Harrison et al. 2000), S = BeppoSAX (Sidoli et al. 2006), X = XMM-Newton (Sidoli et al. 2006;Chernyakova et al. 2006), C = Chandra (Paredes et al. 2007, and this work)) and an ordering index. Observation C1 corresponds to the one performed during our multiwavelength campaign. The exposure is in ks. φ orb and φsup are the orbital and super-orbital phases of the beginning of the observation, computed using T 0 = 43366.275, P orb = 26.4960 days and Psup = 1667 days (Gregory 2002). The photon index (Γ) is obtained from a fit to an absorbed power-law model. The flux (f ) is expressed in 10 −11 ergs cm −2 s −1 , and is integrated between 2-10 keV for RXTE, BeppoSAX and XMM-Newton, between 0.3-10 keV for observation C0 and between 0.5-10 keV for C1. between 0.5-2.0 keV), shown in Figure 5. The Chandra ACIS-I count rate of LS I +61 303 in the energy band of 0.5-10.0 keV varied by a 25%, from ∼ 1.25±0.03 cnts s −1 down to ∼ 1.00 ± 0.03 cnts s −1 , within a timescale of 1-2 hours, showing a clear decline in the final third of the observation. The probability for this variation to be a statistical fluctuation of a constant flux is less than 10 −15 , as obtained from a χ 2 fit of a constant function to the data set. Fast flux variations of up to a 70% within few hours time-scale have been reported by Sidoli et al. (2006). In their case, these were accompanied by variations in the hardness ratio (HR = 2-12 keV/0.3-2 keV) of about 30%, which our data do not confirm. It must be noted, however, that assuming a linear correlation between flux and HR, we should observe less than ∼ 10% variations in the HR during our observation, which is below our statistical error. Table 3 summarizes the results of the observations of LS I +61 303 in the ∼1-10 keV energy range existing in the literature, including that obtained in our Chandra observations. We measure a flux of 1.9×10 −11 ergs cm −2 s −1 and a photon index 1.53. The historical values obtained over several years range between [0.4−2.0]×10 −11 ergs cm −2 s −1 for the flux and [1.25−1.83] for the photon index. Therefore, we observed the source in a particularly high and hard state. This is shown also in Figure 6, where we plot the historical values of the photon index against the flux, in the soft X-ray range. The Pearson's correlation coefficient between these two quantities is r = −0.46, but becomes r = −0.91 when the outlier point obtained by Paredes et al. (2007) is removed from the computation. This is, to our knowledge, the first time that this correlation is shown, and confirms the claim by BeppoSax (Sidoli et al 2006) XMM-Newton (Sidoli et al 2006) Chandra (Paredes et al 2007) Chandra (This work) Fig. 6.-Correlation between the X-ray flux and photon index for the existing observations of LS I +61 303 in the ∼1-10 keV energy band, obtained by using public data from the literature and our own data. (See main text for details). Sidoli et al. (2006) that the source is harder when it is brighter.
3.3. VHE γ-ray Results LS I +61 303 was detected with MAGIC only on October 27 (φ orb = 0.66), with a significance of 4.5σ. For the rest of the nights, no detection above 2σ was found, and we derived the corresponding upper limits to the integral flux (see Table 1).
On October 27, the measured average flux above 300 GeV corresponds to 15% of the Crab Nebula flux at these energies. The VHE γ-ray source is point-like for the MAGIC angular resolution (0.1 • ) and the location is compatible with that of LS I +61 303 . The energy spectrum is well fitted by an unbroken power-law with index α = −2.7 ± 0.5 ± 0.2, where the quoted errors correspond to the statistical and systematic uncertainties, respectively. No significant variations of the absolute flux were detected within that night.
Previous observations of this source with the MAGIC (Albert et al. 2006a and VERI-TAS (Acciari et al. 2008) telescopes have shown a pointlike source, with a peak of ∼ 15% Crab flux intensity at orbital phase φ orb ∼ 0.65, and a spectral index α ∼ −2.6, hence in agreement with the results derived from our observations.

DISCUSSION
The comparison between VHE γ-ray and radio data during the October campaign (see Figure 1) shows a detection at TeV energies for φ orb = 0.66 with a flux level of ∼ 15% of the Crab Nebula flux, during a period when the radio emission is constant at 35 mJy. During the November campaign, however, there is no detection at TeV energies, while the radio data show a peak flux twice as high as in October. Albert et al. (2006a) reported on radio and TeV peaks detected almost simultaneously, while for our campaing we see the TeV peak for a flat and low radio flux and a radio peak for no significant TeV emission. Therefore, we exclude a general TeVradio correlation. A plausible explanation is that the emissions are produced by different particle populations. On the other hand, the detections at X-ray and TeV energies, both at particularly high flux values and one day apart, might point to a correlation between X-ray and TeV fluxes, and hence to the fact that both radiations are produced by the same population of particles. Yet, our data are too scarce to make any solid conclusion at present.
The results obtained by radio imaging at different angular scales show that the size of the radio emitting region of LS I +61 303 is constrained below ∼6 mas (∼12 projected AU), and the presence of persistent jets above this scale is therefore excluded.
We resolve a radio-emitting region at 5.0 GHz with VLBA, extending east-southeast from the brighter, unresolved emitting core. The outflow velocity implied by our observations is at least of ∼0.1 c, in agreement to what previously suggested by Chernyakova et al. (2006), for an over-pressured pulsar wind, although recent hydrodynamical simulations (Bogovalov et al. 2008) show that the shocked pulsar wind could be relativistic.
The comparison of our VLBA image from 25 October with an image at 8.4 GHz obtained at a similar orbital phase but 10 orbital cycles apart shows a high level of similarity between them. It is worth noting that Dhawan et al. (2006) obtained high-resolution radio images of the source along a complete orbital cycle, and found an extended feature whose orientation with respect to the central core varied by 360 • along the orbit. Therefore, the similarity in the morphology and orientation between both images suggests a significant level of periodicity and stability of the physical processes involved in the radio emission. This result points to the fact that the extended radio emission is produced by the interaction of steady flows (from a relativistic pulsar wind, jet and/or stellar wind) rather than by the interaction of such an outflow with wind clumps. We note that if the radio emission is produced by a milliarcsecond scale jet, the required stability and periodic behavior of such a jet is difficult to reconcile with the non-persistent nature of a large scale (∼100 mas), putative relativistic jet, as deduced from our MERLIN observations in combination with those obtained by Massi et al. (2004).
Finally, we combine our X-ray data from 25 October, the TeV data from 26 October and the average VLBA in both nights to produce a quasi-simultaneous multiwavelength spectrum, including radio, X-ray and VHE γray observations, which we display in Figure 7. However, the simultaneous data cover only a small range of the orbital phase of LS I +61 303 . Given the high variability of the physical conditions of the system along the orbit, more simultaneous, multiwavelength data, and particularly involving longer exposure times, orbital phase coverage and redundancy, will shed further light in our understanding of this peculiar object. National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooper-ative agreement by Associated Universities, Inc.