Correlated X-ray and Very High Energy emission in the gamma-ray binary LS I +61 303

The discovery of very high energy (VHE) gamma-ray emitting X-ray binaries has triggered an intense effort to better understand the particle acceleration, absorption, and emission mechanisms in compact binary systems, which provide variable conditions along eccentric orbits. Despite this, the nature of some of these systems, and of the accelerated particles producing the VHE emission, is unclear. To answer some of these open questions, we conducted a multiwavelength campaign of the VHE gamma-ray emitting X-ray binary LS I +61 303 including the MAGIC telescope, XMM-Newton, and Swift during 60% of an orbit in 2007 September. We detect a simultaneous outburst at X-ray and VHE bands, with the peak at phase 0.62 and a similar shape at both wavelengths. A linear fit to the simultaneous X-ray/VHE pairs obtained during the outburst yields a correlation coefficient of r=0.97, while a linear fit to all simultaneous pairs provides r=0.81. Since a variable absorption of the VHE emission towards the observer is not expected for the data reported here, the correlation found indicates a simultaneity in the emission processes. Assuming that they are dominated by a single particle population, either hadronic or leptonic, the X-ray/VHE flux ratio favors leptonic models. This fact, together with the detected photon indices, suggests that in LS I +61 303 the X-rays are the result of synchrotron radiation of the same electrons that produce VHE emission as a result of inverse Compton scattering of stellar photons.

a simultaneous outburst at X-ray and VHE bands, with the peak at phase 0.62 and a similar shape at both wavelengths. A linear fit to the simultaneous X-ray/VHE pairs obtained during the outburst yields a correlation coefficient of r = 0.97, while a linear fit to all simultaneous pairs provides r = 0.81. Since a variable absorption of the VHE emission towards the observer is not expected for the data reported here, the correlation found indicates a simultaneity in the emission processes. Assuming that they are dominated by a single particle population, either hadronic or leptonic, the X-ray/VHE flux ratio favors leptonic models. This fact, together with the detected photon indices, suggests that in LS I +61 303 the X-rays are the result of synchrotron radiation of the same electrons that produce VHE emission as a result of inverse Compton scattering of stellar photons. Subject headings: binaries: general -gamma rays: observations -stars: emission-line, Be -stars: individual (LS I +61 303) -X-rays: binaries -X-rays: individual (LS I +61 303) 1. INTRODUCTION LS I +61 303 is one of the few X-ray binaries that have been detected in very high energy (VHE) gamma rays (see, e.g., Paredes 2008 for a recent review). It is a high-mass X-ray binary system located at a distance of 2.0±0.2 kpc (Frail & Hjellming 1991). The system contains a rapidly rotating early type B0 Ve star with a persistent equatorial decretion disk and mass loss, and a compact object with a mass between 1 and 4 M ⊙ orbiting it every ∼26.5 d in an eccentric orbit (see Casares et al. 2005;Grundstrom et al. 2007;Aragona et al. 2009, and references therein). LS I +61 303 was classified as a microquasar based on structures detected from five to several tens of milliarcseconds with the European VLBI Network (EVN) and MERLIN (Massi et al. 2004, and references therein), although analysis of later EVN and MERLIN data sets does not reveal the presence of such structures (Dhawan et al. 2006;Albert et al. 2008a). Very Long Baseline Array (VLBA) images with a resolution of 1 mas (2 AU at the source distance) obtained during a full orbital cycle show an elongated structure that rotates as a function of the orbital phase (Dhawan et al. 2006). Later VLBA images show repeating structures at the same orbital phases, reinforcing the idea that the milliarcsecond structure depends on the orbital phase (Albert et al. 2008a). This may be consistent with a model based on the interaction between the relativistic wind of a young non-accreting pulsar and the wind/decretion disk of the stellar companion (Dubus 2006a), as occurs in PSR B1259−63 (but confining the pulsar wind may be problematic; see Bogovalov et al. 2008).
LS I +61 303 shows periodic non-thermal radio outbursts on average every P orb =26.4960±0.0028 d, with the peak of the radio emission shifting progressively from phase 0.45 to 0.95, using T 0 =JD 2,443,366.775, with a modulation period of 1667±8 d (Gregory 2002). According to the most precise orbital parameters, the periastron takes place at phase 0.275 and the eccentricity of the orbit is 0.537 ± 0.034 (Aragona et al. 2009).
LS I +61 303 has been observed several times in Xrays (see Smith et al. 2009 and references therein). It generally displays quasi-periodic X-ray outbursts, with the maximum occurring between orbital phases 0.4 and 0.8 (Goldoni & Mereghetti 1995;Taylor et al. 1996;Paredes et al. 1997;Harrison et al. 2000;Esposito et al. 2007), although the lack of a sensitive long-term monitoring has prevented to search for a super-orbital modulation. The source also shows short-term variability on timescales of hours (Sidoli et al. 2006).
At VHE gamma rays, LS I +61 303 has been clearly detected several times both by MAGIC (Albert et al. 2006(Albert et al. , 2008a(Albert et al. , 2009) and VERITAS (Acciari et al. 2008(Acciari et al. , 2009). The source also displays VHE gamma-ray periodicity, with min-ima taking place near periastron, where only upper limits were found in MAGIC observations, and maxima occurring on average at phase 0.6-0.7, although the source has also shown a second peak at phase 0.84 in one of the cycles (Albert et al. 2009). There are indications of correlated X-ray/VHE emission, based on non-simultaneous data taken more than 6 hr apart (Albert et al. 2009), and 1 d apart (Albert et al. 2008a).
The lack of a systematic behavior from cycle to cycle at X-ray and VHE bands, and the occurrence of short-term variability in the X-ray flux, has hampered the definitive detection of an X-ray/VHE correlation from the comparison of non-simultaneous data. We therefore conducted a multiwavelength campaign in 2007 September, covering the epoch of maximum VHE emission. Here, we report the first simultaneous VHE and X-ray observations of LS I +61 303 obtained with the MAGIC Cherenkov telescope and the XMM-Newton and Swift X-ray satellites that show correlated emission in both energy bands. We also briefly comment on radio and optical spectroscopic observations obtained during the campaign.

VHE Gamma Rays
The VHE gamma-ray observations were performed from 2007 September 4 to 21 using the MAGIC telescope on the Canary Island of La Palma (28.75 • N, 17.86 • W, 2225 m a.s.l.), from where LS I +61 303 is observable at zenith angles above 32 • . The essential parameters of MAGIC are a 17 m diameter segmented mirror of parabolic shape (currently the largest dish for an air Cherenkov telescope), an f /D of 1.05, and an hexagonally shaped camera of 576 hemispherical photo multiplier tubes with a field of view of 3.5 • diameter. MAGIC can detect gamma rays from 60 GeV to several TeV, and with a special setup the trigger threshold can be reduced to 25 GeV (Aliu et al. 2008). Its energy resolution is ∆E/E = 20% above energies of 200 GeV. The current sensitivity is 1.6% of the Crab Nebula flux for a 5σ detection in 50 hr of on-source time. The improvement compared to the previous sensitivity was achieved by installing new 2 GHz FADCs (Albert et al. 2008b).
The total observation time was 58.8 hr, including 39.6 hr in dark conditions and 19.2 hr under moonlight or twilight (Britzger et al. 2009). The range of zenith angles for these observations was [32 • , 50 • ], with 96% of the data having zenith angles below 43 • . The final effective total observation time is 54.2 hr. Table 1 gives the MJD, the effective observation time, and the orbital phase for each of the MAGIC observations. The observations were carried out in wobble mode (Fomin et al. 1994), i.e., by alternately tracking two positions at 0.4 • offset from the source position. This observing mode provides a reliable background estimate for point-like sources such as LS I +61 303. The data analysis was performed using the standard MAGIC analysis and reconstruction software (Albert et al. 2008c). The quality of the data was checked and data taken with anomalous event rates (very low atmospheric transmission, car light flashes, etc.) were rejected following the standard procedure, as previously done in Albert et al. (2009). From the remaining events, image parameters were calculated (Hillas 1985). In addition, the time parameters described in Aliu et al. (2009) were calculated as well. For the γ/hadron separation, a multidimensional classification procedure based on the image and timing parameters with the Random Forest method was used (Albert et al. 2008d). We optimized the signal selection cuts with contemporaneous Crab Nebula data taken at the same zenith angle range. The energy of the primary γ-ray was reconstructed from the image parameters using also a Random Forest method. The differential energy spectrum was unfolded taking into account the full instrumental energy resolution (Albert et al. 2007). We estimate the systematic uncertainty to be about 30% for the derived integral flux values and ±0.2 for the obtained photon index. For more details on the systematic uncertainties present in the MAGIC data, see Albert et al. (2008c).

X-rays
We observed LS I +61 303 with XMM-Newton during seven runs from 2007 September 4 to 11, lasting from 12 to 18 ks (see Table 2), amounting to a total observation time of 104.3 ks. The EPIC pn detector used the Large Window Mode, while the EPIC MOS detectors used the Small Window Mode. All detectors used the medium thickness optical blocking filter. Data were processed using version 8.0.0 of the XMM-Newton Science Analysis Software (SAS). Known hot or flickering pixels were removed using the standard SAS tasks. Further cleaning was necessary to remove from the data set periods of high background corresponding to soft proton flares, reducing the net good exposure durations to 67.0 and 92.6 ks for the pn and MOS detectors, respectively.
Cleaned pn and MOS event files were extracted for spectral analysis. Source spectra were extracted from a ∼70 ′′ radius circle centered on the source (point-spread function of 15 ′′ ) while background spectra were taken from a number of source-free circles with ∼150 ′′ radius (three for the pn de-tector, four for MOS1, and five for MOS2). The extracted spectra were analyzed with XSpec v12.3.1 (Arnaud 1996). The spectra were binned so that a minimum of 20 counts per bin were present and energy resolution was not oversampled by more than a factor 3. An absorbed power-law function (wabs * powerlaw) yielded satisfactory fits for all observations. De-absorbed fluxes in the 0.3-10 keV range were computed from the spectral fits.
Additional observations of 2-5 ks each, consisting of several pointings of 0.2-1.0 ks, were obtained with the Swift/XRT from 2007 September 11 to 22. The total observation time was 28.5 ks. The Swift data were processed using the FTOOLS task xrtpipeline (ver. 0.12.1 under HEASoft 6.6). The spectral analysis procedures were the same as those used for the XMM-Newton data. Since the hydrogen column density was poorly constrained, with typical values of (0.5 ± 0.2) × 10 22 cm −2 , we fixed it to 0.5 × 10 22 cm −2 , a typical value for LS I +61 303 also found in the XMM-Newton fits and close to the average value for the Swift fits.
To search for short-term X-ray variability, we also extracted 0.3-10 keV background-subtracted light curves for each observation, binned to 100 s for XMM-Newton and at half-time of the sparse 0.2-1.0 ks Swift pointings within each observation. For XMM-Newton, we considered the sum of the count rate in the three detectors. From these light curves, we computed the degree of variability as the standard deviation with respect to the mean of the count rate divided by this mean, and the significance of this variability computed from the χ 2 probability of the count rate being constant. We also computed hardness ratios as the fraction between the count rates above and below 2 keV.

VHE Gamma Rays
The measured fluxes of LS I +61 303 at energies above 300 GeV are listed in Table 1, and the corresponding light curve is shown in Figure 1 (top). The periodically modulated peak emission (Albert et al. 2009) is prominently seen as the highest flux value of these observations at phase 0.62 (followed by a smaller flux at phase 0.66). In addition to this established feature we measured, as in 2006 December, significant flux on a nightly basis at phase 0.85. We note that there is significant flux in the phase range 0.8-1.0 at the level of (5.2 ± 1.0) × 10 −12 cm −2 s −1 , which is compatible with the 2σ upper limit we obtained for the shorter 2006 observations (Albert et al. 2009).
We also determined the differential energy spectrum from the ∼10 hr of observations conducted in the orbital phase range 0.6-0.7. A power-law fit yields: , compatible within errors with our previous measurements (Albert et al. 2006(Albert et al. , 2009).

X-rays
We summarize in Table 2 the parameters of the spectral fits and variability obtained for both the XMM-Newton and Swift/XRT data sets. All X-ray fluxes quoted hereafter are de-absorbed. We note that there is no significant hardness  ratio change within each of the observations, in contrast to what was found in the XMM-Newton observations reported by Sidoli et al. (2006). In principle this implies that, for each observation, the flux and corresponding uncertainty obtained from the spectral fit is a good estimate of the flux during the whole observation. However, moderate count-rate variability on timescales of ks is present in most observations, ranging from 6% to 17% for the XMM-Newton data and from 9% to 25% for the Swift data (see Table 2). The rms of this variability should be considered in addition to the statistical uncer-tainty when providing a flux measurement spanning several ks. We converted this count-rate variability into flux variability by multiplying the degree of variability defined in Section 2.2 with the fluxes coming from the spectral fits. Since no spectral change is detected within each observation, we added this flux variability in quadrature to the spectral fits flux errors. This procedure provides the more realistic total flux uncertainties quoted in parentheses in Table 2, and used hereafter.
We show in Figure 1 (bottom) the 0.3-10 keV light curve of LS I +61 303. The source displays a steady flux during the first four observations and shows a steep increase at phase 0.62, which is followed by a slower decay up to phase 0.69 (XMM-Newton) and probably up to phase 0.72 (Swift). The behavior is very similar to the one seen at VHE gamma rays, although at X-ray energies the baseline has a significant flux. Later on there is a significant increase of the X-ray flux up to phase 0.8. This high flux is detected with a sparse sampling up to phase 0.9, and the source goes back to its baseline flux at phase 1.0. This high X-ray flux between phases 0.8 and 1.0 occurs when the source is also detected at VHE gamma rays.

X-Ray/VHE Gamma-Ray Correlation
A clear correlation between the X-ray and VHE gammaray emissions is seen during the outburst, with a simultaneous peak at phase 0.62 (see Figure 1). To study the significance of this correlation, we selected the X-ray data sets that overlap with MAGIC observations. There are six overlapping MAGIC/XMM-Newton data sets, for which strictly simultaneous observations range from 3.3 to 3.9 hr. For the four overlapping MAGIC/Swift data sets, the strictly simultaneous observations range from 2.2 to 4.1 hr (although the Swift runs have gaps). We plot in Figure 2 the X-ray fluxes against the VHE fluxes (from Tables 1 and 2) for all 10 simultaneous pairs, which are marked with arrows in Figure 1. The linear correlation coefficient for the six simultaneous MAGIC/XMM-Newton pairs that trace the outburst is r = 0.97. For the ten simultaneous pairs, we find r = 0.81 (a |r| larger than that has a probability of about 5 × 10 −3 to be produced from independent X-ray and VHE fluxes). Minimizing χ 2 we obtain χ 2 = 7.68 for 8 degrees of freedom and the fol-  Figure 1, have been considered. Error bars correspond to a 1σ confidence level in all cases. The solid line represents a χ 2 linear fit to all data points.
However, as can be seen in Figure 2, the flux uncertainties are relatively large. This calls for a test of the reliability of the correlation strength considering the errors of individual data points. To this end, we use the z-Transformed Discrete Correlation Function (ZDCF), which determines 68% confidence level intervals for the correlation coefficient from unevenly sampled data (see, e.g., Edelson & Krolik 1988;Alexander 1997). The Fisher z-transform of the linear correlation coefficient is used to estimate the 68% confidence level interval. Applying the ZDCF to the X-ray and VHE light curves reported here we obtain the following uncertainties for the linear correlation coefficient: r = 0.81 +0.06 −0.21 . The sensitivity of the MAGIC telescope requires observation times of several hours to be able to get few-sigma detections for the VHE fluxes from LS I +61 303. During such long periods, the X-ray fluxes have quite large intrinsic variability (see Table 2). This results in quite large uncertainties on both the VHE and X-ray fluxes that restrict the ability to determine a possible correlation out of the outburst.

DISCUSSION
We have discovered an X-ray/VHE gamma-ray correlation in the gamma-ray binary LS I +61 303 based on simultaneous multiwavelength data obtained with MAGIC, XMM-Newton, and Swift during a single orbital cycle. The correlation is mainly due to the very similar trends of the detected flux during the outburst around orbital phase 0.62, while the uncertainties prevent to be sure about the existence of the correlation outside the outburst. Given the variability in the X-ray flux (up to 25% on hour scales), it is necessary to program simultaneous observations to perform correlation studies in LS I +61 303. This conclusion has also been reached recently by the VERITAS Collaboration, when reporting the lack of Xray/VHE correlation based on contemporaneous data with a VHE sampling that is not dense enough (Acciari et al. 2009). Although a similar X-ray/VHE correlation has been obtained for the VHE gamma-ray emitting X-ray binary LS 5039, this result was based on non-simultaneous data acquired years apart (Aharonian et al. 2005Takahashi et al. 2009).
The VHE emission within a binary system can suffer photon-photon absorption via pair creation, mainly with the stellar optical/ultraviolet photons. In LS I +61 303, this absorption is only expected to be significant toward the observer for E > 300 GeV just before periastron (Dubus 2006b;Bednarek 2006;Sierpowska-Bartosik & Torres 2009), a phase range not explored here. In addition, the quoted Xray fluxes are already de-absorbed (and the hydrogen column densities and the associated errors are low). Overall, there are no absorption effects to be considered. Therefore, the Xray/VHE correlation we have found for LS I +61 303 cannot be an artifact due to variable absorption toward the source. This indicates that the emission processes at both wavelengths occur at the same time and are probably the result of a single physical mechanism. In this context, it is reasonable to assume that the X-ray and VHE emissions are produced by a single particle population.
It is interesting to note that the MAGIC spectrum in the 0.6-0.7 phase range yields ∼ 11 × 10 −12 erg cm −2 s −1 for E > 300 GeV, while the X-ray flux is ∼ 19 × 10 −12 erg cm −2 s −1 . Therefore, the total X-ray flux is approximately twice the VHE flux. However, if we subtract an apparent baseline Xray flux of 10 × 10 −12 erg cm −2 s −1 , the resulting X-ray flux is similar to the total VHE flux in the phase range 0.6-0.7. If the radiation mechanisms are dominated by a single particle population, the X-ray/VHE correlation and the smaller/similar VHE fluxes favor leptonic models. In hadronic models, the X-ray emitting e ± and the VHE photons would come from the same protons (for reasonable values of the magnetic field), and the luminosity of the e ± radiation should be 1/2 that of VHE gamma-rays (see Figure 5 of Kelner et al. 2006 for reasonable proton energy distributions) unlike it is observed. In addition, the inverse Compton (IC) cooling channel is less efficient than the synchrotron channel to produce the detected X-ray emission for reasonable values of the magnetic field (see Takahashi et al. 2009 for a similar discussion for LS 5039). This clearly suggests that the X rays are the result of synchrotron radiation of the same electrons that produce VHE emission as a result of IC scattering of optical/ultraviolet stellar photons.
The observed photon indices of the simultaneous X-ray and VHE spectra are consistent with one population of electrons following a power-law energy distribution with index ∼2.1. These electrons would produce X-rays via synchrotron and VHE photons via IC with an interaction angle π/2. We note that an electron index of ∼2.1 is too hard if synchrotron cooling dominates in the X-ray range, since it implies an injected electron index of 1.1. On the other hand, dominant adiabatic cooling implies an injection index of 2.1, a more reasonable value (see, for example, the discussion in Takahashi et al. 2009). Although small changes in the VHE spectrum would be expected due to variations in the IC interaction angle or the electron index (as seen in X-rays), at present they are not detectable due to the large uncertainties of the VHE photon index.
Finally, we note that contemporaneous radio light curves obtained with RATAN, VLBA images, and Hα spectroscopy are consistent with previous results (Gregory 2002;Dhawan et al. 2006;Zamanov et al. 1999). Details on these observations will be reported elsewhere. Therefore, the Xray/VHE correlation occurred when the source was showing a standard behavior in both its outflow (radio) and decretion disk (Hα line).