Periodic very high energy gamma-ray emission from LS I +61 303 observed with the MAGIC telescope

The MAGIC collaboration has recently reported the discovery of gamma-ray emission from the binary system LS I +61 303 in the TeV energy region. Here we present new observational results on this source in the energy range between 300 GeV and 3 TeV. In total 112 hours of data were taken between September and December 2006 covering 4 orbital cycles of this object. This large amount of data allowed us to produce an integral flux light curve covering for the first time all orbital phases of LS I +61 303. In addition, we also obtained a differential energy spectrum for two orbital phase bins covering the phase range 0.5<phi<0.6 and 0.6<phi<0.7. The photon index in the two phase bins is consistent within the errors with an average index Gamma=2.6+-0.2_{stat}+-0.2_{sys}. LS I +61 303 was found to be variable at TeV energies on timescales of days. These new MAGIC measurements allowed us to search for intra-night variability of the VHE emission; however, no evidence for flux variability on timescales down to 30 minutes was found. To test for possible periodic structures in the light curve, we apply the formalism developed by Lomb and Scargle to the LS I +61 303 data taken in 2005 and 2006. We found the LS I +61 303 data set to be periodic with a period of (26.8+-0.2) days (with a post-trial chance probability of 10^{-7}), close to the orbital period.


INTRODUCTION
The γ -ray binary system LS I +61 • 303 is located at a distance of ∼2 kpc and is composed of a compact object of unknown nature (neutron star or black hole) orbiting a Be star in a highly eccentric orbit (e = 0.72 ± 0.15 or e = 0.55 ± 0.05 following Casares et al. (2005) and Grundstrom et al. (2007), respectively).
The orbital period of the system is 26.4960 days long (Gregory 2002). The periastron passage, derived from the optical spectra, is found to be at phase φ = 0.23 ± 0.02 in Casares et al. (2005) and φ = 0.301 ± 0.011 in Grundstrom et al. (2007), adopting a zero phase at T 0 = JD2443366.775.
Radio outbursts are observed every orbital cycle at phases varying between 0.45 and 0.95 with a 4.6 years modulation (Gregory 2002). Radio-imaging techniques have shown extended, radio-emitting structures with angular extensions of ∼0.01-0.1 arcsec, where the radio emission originates in a twosided, possibly precessing, relativistic jet (β/c = 0.6) (Massi et al. 2004). These extended radio structures have led some authors to adopt the microquasar scenario to explain the nonthermal emission in LS I +61 • 303 (e.g., Bosch-Ramon et al. 2006;Bednarek 2006). Recent high-resolution Very Long Baseline Array (VLBA) measurements show a complex and changing morphology different from what is expected for a typical microquasar jet (see, e.g., radio images in Dhawan et al. 2006;Albert et al. 2008b). Furthermore no solid evidence for the presence of an accretion disk (i.e., a thermal X-ray component) has been observed (Chernyakova et al. 2006). This seems to favor a scenario in which the nonthermal emission in LS I +61 • 303 is powered by the interaction between a pulsar and the primary star winds (Maraschi & Treves 1981). Romero et al. (2007) presented smoothed particle hydrodynamics (SPH) simulations of the system to discuss the possible compositions, although again results are not compelling evidence of one or the other.
At higher energies, LS I +61 • 303 was found to be spatially coincident with the EGRET γ -ray source 3EG J0241 + 6103 (Kniffen et al. 1997). Variable emission at TeV energies was observed with the MAGIC telescope (Albert et al. 2006) and was recently confirmed by VERITAS (Acciari et al. 2008). The system showed the peak TeV γ -ray flux at phase φ ∼ 0.65, while no very high-energy emission was detected around the periastron passage.
Here we present new MAGIC telescope observations of LS I +61 • 303. We briefly discuss the observational technique and the data analysis procedure, investigate the very high energy (VHE) γ -ray spectrum during the high emission phase of the source, and put the results into perspective with previous VHE γ -ray observations of this system. Finally, we analyzed the temporal characteristics of the TeV emission and find a periodic modulation of the signal with the orbital period.

OBSERVATIONS
The observations were performed from 2006 September to December using the MAGIC telescope on the Canary Island of La Palma (28. • 75N, 17. • 86W, 2225 m a.s.l.), from where LS I +61 • 303 is observable at zenith distances above 32 • . The telescope operates in the energy band from 50-60 GeV (trigger 27 Deceased. threshold at zenith angles less than 30 • ) up to tens of TeV, with a typical energy resolution of 20%-30%. The accuracy in reconstructing the direction of the incoming γ -rays is about 0. • 1, depending on the energy. A detailed description of the telescope performance can be found in Albert et al. (2008c).
The data on LS I +61 • 303 were taken between 2006 September 15 and 2006 December 28 covering four orbital periods of the system. In total, 120 hr of data were taken at zenith angles between 32 • and 55 • , with ∼ 97% of the data below 44 • . After preselection of good quality data, a total of 112 hr of data remained for the analysis. About 17% of these were recorded under moderate moonlight conditions. Due to the different observation conditions such as bad weather, too bright moon or too large zenith angle, the data set was not uniform with the orbital phase. In Table 1, the observation times of the analyzed data are summarized.
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 actual source position. This observation mode allows for a reliable background estimate for point-like objects such as LS I +61 • 303.

DATA ANALYSIS
The data analysis was carried out using the standard MAGIC analysis and reconstruction software (Albert et al. 2008c). The images were cleaned by requiring a minimum number of ten photoelectrons (core pixels) and five photoelectrons (boundary pixels), see, e.g., Fegan (1997). The quality of the data was checked and bad data such as accidental noise triggers or data taken during adverse conditions (very low atmospheric transmission, car light flashes, etc.) were rejected. From the remaining events, image parameters were calculated (Hillas 1985).
For the γ /hadron separation, a multidimensional classification procedure based on the Random Forest method (Albert et al. 2008a;Bock et al. 2004) was used. For every event a parameter called hadronness (h) is derived, based on the values of the event's image parameters. The hadronness denotes the probability that an event is a hadronic-induced (background) event. The final separation was achieved by a cut in h which was determined by requiring 80% of the simulated Monte Carlo (MC) γ -ray events to be kept. In addition to the cut in h, a geometrical cut in the squared angular distance of the assumed source position to the shower direction axis (θ 2 cut) was performed so that 70% of all simulated MC γ -ray events from a point-like source are left after the cut. The cut efficiencies were chosen so that the same cut efficiencies yield the highest significance of a Crab nebula data sample taken under similar conditions as the LS I +61 303 data. The same cut procedure was applied to the final LS I +61 • 303 sample. The energy of the primary γ -ray was reconstructed from the image parameters, also using a Random Forest method, leading to an assigned estimated energy for each reconstructed γ -ray event. The differential energy spectrum is unfolded taking into account the full instrumental energy resolution (Albert et al. 2007). For the integral flux calculation of the light curves we used fixed cuts (for all energies) in hadronness and θ 2 . In the case of the energy spectrum determination, we derived fixed hadronness and θ 2 cuts for each energy bin.
The main contributions to the systematic error of our analysis are the uncertainties in the atmospheric transmission, the reflectivity (including stray-light losses) of the mirror and the light catchers, the photon to photoelectron conversion calibration, and the photoelectron collection efficiency in the photomultiplier front end (Albert et al. 2008c). Also MC uncertainties in the detector simulation and systematic uncertainties from the analysis methods contribute significantly to the overall error.

Note.
a Quoted for cases where the flux significance is 2σ (following Rolke et al. 2005).
All errors in this paper are statistical errors, unless otherwise stated explicitly. In addition, there is a 30% systematic uncertainty on flux levels and 0.2 on the spectral photon index.
MAGIC has the capability to operate under moderate moonlight. This permits us to increase the duty cycle by up to 28%, thus considerably improving the sampling of transient sources.
In particular, 17% of the data used in this analysis were recorded under moonlight. The nights which were partly taken under moonlight conditions are labeled with an asterisk in Table 1. For these days we estimate an increased systematic error of ∼ 40% instead of the ∼ 30% in the case of the dark night observations. All spectra are derived from data which were only taken under dark night conditions and thus no additional error is present in the obtained parameters. Figure 1 presents the gamma-ray flux above 400 GeV measured from the direction of LS I +61 • 303 as a function of the orbital phase of the system for the four observed orbital cycles. The probability for the distribution of measured fluxes to be a statistical fluctuation of a constant flux (obtained from a χ 2 fit to the entire data sample) is 4.4×10 −6 (χ 2 /dof = 108.9/51). In all orbital cycles, significant detections (S > 2σ ) occurred during the orbital phase bin 0.6-0.7. The highest measured fluxes are dominantly found in this phase bin. Among those nights around phase 0.65, the night MJD 54035.11 shows the maximum flux, with a statistical significance of 4.5 σ .

Light Curve
At the periastron passage (phase 0.23, according to Casares et al. 2005) the flux level is always below the MAGIC sensitivity and we derive an upper limit with 95% confidence level of 4 × 10 −12 cm −2 s −1 (MJD 53997). If we take for the periastron passage the phase value 0.3 as obtained by Grundstrom et al. (2007), we obtain on MJD 53999 a flux of F (E > 400 GeV) = (5.3 ± 2.4) × 10 −12 cm −2 s −1 . This is not enough to provide certainty for a confident detection at this early phase though a 2.2σ hint is obtained in the phase range 0.2-0.3, possibly a fluctuation (see Table 2). In addition, since the real value of the periastron passage is still the subject of debate the averaged flux value between the phase bin 0.2-0.3 can be used as an upper limit to the emission at periastron. Thus the flux must be less than F (E > 400 GeV) = 2.2 × 10 −12 cm −2 s −1 at the 95% confidence level.
The flux quoted above corresponds to 7% of the integral Crab nebula flux in the same energy range. The mean flux for all other phase bins can be found in Table 2. This is well in agreement with the flux measured by MAGIC in the first campaign (Albert et al. 2006). The data we present here have been reanalyzed with an improved energy estimation.
One very interesting peculiarity found in the light curve is that a second peak in the flux level is seen in the last observed period at phase 0.84. The flux is (16 ± 4) × 10 −12 cm −2 s −1 which is at a similar level compared to the maximum flux detected in phase ∼0.65 (MJD 54035.11). This high flux was not seen in similar phases in any previous cycle, where only upper limits could be set (see Table 1). Six hours after our measurement, the data of the Swift X-ray satellite showed a high flux (0.25±0.01 counts s −1 ) at phase 0.85 (Esposito et al. 2007). These Swift observations did not cover the same orbital phase in any other orbit. Beside this second peak the main X-ray emission peak is found between the phases 0.5-0.8 (∼ 0.24 counts s −1 , Esposito et al. 2007) in exactly the same phase bin where LS I +61 • 303 is detected by MAGIC. This is a hint for a correlated X-ray/TeV emission.
Our measurement is in agreement with the published VERITAS measurements (Acciari et al. 2008), i. e., that LS I +61 • 303 is detected at TeV energies in the phase range 0.5-0.8. The MAGIC and VERITAS data are not strictly simultaneously taken and VERITAS did not observe LS I +61 • 303 in 2006 December where the second peak occurred.
Due to our long observation time and dense sampling of orbital phases we obtained the currently most detailed light curve of LS I +61 • 303.

Spectral Studies
As seen from the light curve ( Figure 1) LS I +61 • 303 is a very variable source which shows high flux levels only at some orbital phases. For the phases 0.5 < φ < 0.6 and 0.6 < φ < 0.7, where we measured significantly high flux levels, we were able to determine differential energy spectra. In both cases, the energy spectra obtained are compatible with pure power laws. In the case of the phase bin 0.6 < φ < 0.7, a power-law fit gives The spectrum of phase bin 0.5 < φ < 0.6 (green), phase bin 0.6 < φ < 0.7 (red), and the spectrum obtained from previous MAGIC measurements (blue dashed; Albert et al. 2006) are shown. The spectra from phase bin 0.6 < φ < 0.7 and the previous MAGIC measurements can be well described by a simple power law with photon index 2.6±0.2. The spectral slope of phase bin 0.5 < φ < 0.6 is compatible with these results within the errors. with a reduced χ 2 /dof = 5.22/5. The spectral fit parameters are in excellent agreement with those previously reported ones using MAGIC (Albert et al. 2006). In addition, we derived the differential energy spectra for the two nights with a signal of > 4.5σ significance which are part of the same phase bin 0.6 < φ < 0.7. Both spectra are also well described by a pure power law (see Table 3). No evidence for spectral variations has been found.
In the case of phase bin 0.5 < φ < 0.6, we obtained with a reduced χ 2 /dof = 1.42/4, showing that the spectral shape is well compatible with a simple power law.
The energy spectra of the two phase bins together with the power-law fits are shown in Figure 2. The corresponding fit parameters and their errors are also shown in Table 3.
The spectral indices of the fitted power laws, for both phase bins and the single night spectra, are compatible within their errors indicating that no significant spectral changes happened between the different phase bins and between the different orbital cycles. In the phase bins 0.0 < φ < 0.5 and 0.7 < φ < 1.0 the γ -ray flux is too low to derive meaningful differential energy spectra.
Another possibility to search for spectral variation is by means of the hardness ratio HR which we define as the ratio of the integral flux between 400 GeV and 900 GeV and above 900 GeV. The HR plotted against the total integral flux above 400 GeV for each night with a signal above 2σ significance is shown in Figure 3. The requirement of the 2σ significance of the signal is to minimize systematic effects on the calculation of the correlation coefficient. We do not find any clear correlation between the HR and the flux level. Thus we do not find any change in the spectral behavior in nights where LS I +61 • 303 is detected at modest significance. The spectral studies on LS I +61 • 303 exhibit that the spectrum is soft (compared to other galactic sources) during all phases in which LS I +61 • 303 is detected at TeV energies. While the flux level changes on timescales of days and reaches a maximum flux (detected above 3σ ) of 15.6 × 10 −12 cm −2 s −1 the source shows a spectral photon index of 2.6 ± 0.2, compatible with being constant.

Search for Intranight Variability
LS I +61 • 303 was found to be variable on timescales of days in the previous observational campaign by MAGIC (Albert et al. 2006). A still open question is whether LS I +61 • 303 also shows variability on shorter timescales. We investigated the data for all nights with a significant flux level (F (E > 400 GeV) > 4 × 10 −12 cm −2 s −1 ) with respect to intranight variability on timescales ranging from 30 to 75 min, in steps of 15 min. This yields 16 suitable nights. Among these, the longest one was MJD 54086.95 (phase 0.61). Its intranight light curve is shown in Figure 4, where each bin has ∼ 1 hr width. The light curve is fitted with a constant flux level with probability of 90% (χ 2 /dof = 1.08/4). For the rest of the nights and tested timescales, the post-trial probabilities of being chance fluctuations of a constant flux are all above 32%. We conclude that the VHE fluxes are compatible with being constant on timescales of 30-75 min within the MAGIC sensitivity.

Search for Periodicity
The emission of the binary system changes periodically in radio, optical, and X-rays, and the modulation is associated with the orbital period of the binary system. At higher energies, EGRET measurements showed hints for variable γ -ray emission (Tavani et al. 1998), although no periodicity could be established with these data. One of the aims of the MAGIC long observational campaign on LS I +61 • 303 was to search for periodic VHE γ -ray emission.
In order to maximize the sensitivity and accuracy of the timing analysis, we used the data presented in this work together with the data taken in the first campaign (Albert et al. 2006), with an observation time of 54 hr and covering six orbital cycles.
The periodicity analysis was carried out using the formalism developed by Lomb and Scargle (Lomb 1976;Scargle 1982). This formalism allows one to analyze unevenly sampled data while still keeping the simple exponential probability density function (PDF, P (z > z 0 ) = e −z 0 ) for Gaussian white noise (GWN) as valid for the classical Fourier analysis of evenly sampled data. A remaining problem caused by the uneven data sampling is that the independent Fourier spacing is broken, i.e., even a single-frequency component can result in a complex power spectrum with a large number of aliasing peaks. Another problem is that the mean and variance values that enter the Lomb-Scargle periodogram have to be estimated from the data themselves.
A practical method to determine the chance probabilities is the following (see, e.g., Frescura et al. 2007).

A large number of random data series is constructed with
a MC simulation of random fluxes while keeping the sampling times fixed. In the lower panel, we show the periodograms after subtraction of a sinusoidal signal (see Figure 7) at the orbital period (yellow line) and a sinusoidal plus a Gaussian wave form (blue line). Vertical dashed line corresponds to the orbital frequency. Inset: zoom around the highest peak, which corresponds to the orbital frequency (0.0377d −1 ). Its post-trial probability is nearly 10 −7 (see Figure 5). The IFS is also shown.
2. For each random series, we construct a periodogram, sampling it for a preselected group of frequencies. 3. For each frequency, we compare the periodogram derived from the real data set with the PDF obtained from the simulated random series, in order to empirically determine the (pre-trial) chance probability. 4. The overall (post-trial) chance probability is computed according to the following generalization: for each simulated data series we inspect the corresponding periodogram, identify the highest Fourier power that occurs at any of the preselected frequencies, and use this value to construct the post-trial PDF. It should be noted that this constructed PDF is based on the null hypothesis of GWN.
Integration of the post-trial PDF gives the complementary Cumulative Distribution Function (cCDF) which is used to determine the (post-trial) chance probability for a given Fourier power value.
In Figure 5, we show the empirical post-trial cCDF of the Lomb-Scargle power, estimated via Monte Carlo simulation of random fluxes. The expected cCDF above a spectral peak z 0 is F (z > z 0 ) = 1 − (1 − e −z 0 ) M , where M is the number of independent frequencies. By fitting the PDF for LS I +61 • 303, we obtain a probability of 75% (χ 2 /dof= 263.9/279) and a number of independent frequencies of M = 550.8 ± 0.6. This result is used to estimate the chance probability of the Lomb-Scargle powers.
In Figure 6 (middle panel), we show the Lomb-Scargle periodogram for an almost (up to detector related effects) independent background sample, obtained simultaneously with the LS I +61 • 303 data (see below for the time intervals). The highest obtained power is 7.5, which yields a probability of 0.3. Thus we obtain no significant probability peaks for The black curve is a fit to a sinusoidal signal. We also fitted a sinusoidal signal plus a Gaussian component (blue dotted line), which adjusts better to the data (fit parameters are given in the inset). The vertical dashed lines mark the two measurements of the periastron passage (according to Casares et al. 2005 andGrundstrom et al. 2007).
any of the scanned frequencies. We construct the background periodograms using an OFF sample obtained simultaneously with the LS I +61 • 303 data. This is treated in different ways for campaigns I and II: for campaign I the data are taken in ON/ OFF mode, i.e., ON and OFF data are taken at different times.
To obtain a simultaneous OFF data sample we use an anticut in alpha. (As the alpha cut for ON data is α < 6 • , we consider the events with 20 • < α < 30 • .) With these numbers of γ -ray candidates and the collection area (given by MC simulations) and observation time, we estimate the corresponding flux. For campaign II the data are taken in wobble mode and thus we use the signal candidates recorded by the antisource position, and estimate the background with the other two wobble regions to evaluate the flux. The background rate in both campaigns is very stable, for campaign I (7.28 ± 0.08) × 10 −3 and for campaign II (5.34 ± 0.07) × 10 −3 s −1 , and thus presents a reliable background sample. We apply the Lomb-Scargle test to the LS I +61 • 303 data and obtain the periodogram shown in Figure 6 (upper panel). The periodogram is performed with the LS I +61 • 303 integral flux above 400 GeV, measured in a time interval [t i − Δt 2 , t i + Δt 2 ], for Δt = 15 min and i = 0, . . . 717 data points. The overall time range is 442 days, which yields an independent Fourier spacing (IFS) of ν IFS = 1/T = 0.0023 d −1 . We scanned the frequency range 0.0023-0.25 d −1 with an oversampling factor of 5.
A maximum peak in the Lomb-Scargle periodogram is clearly seen at frequency ν = 0.0373 d −1 , for which we obtain a power of ∼ 22, corresponding to a post-trial chance probability of 2 × 10 −7 .
Several less prominent but significant peaks are also detected for other frequencies (e.g., 0.041 d −1 with probability 10 −5 ). Those peaks are related to the signal, since they are not present in the contemporaneous background sample (Figure 6, middle panel). These are aliasing peaks of the orbital period of LS I +61 • 303 caused by the various gaps in the data set.
The observational bias due to the moon cycle cannot be responsible for the observed peak since this period should otherwise be also present in the background periodogram. The data folded with the peak frequency (ν = 0.0377 d −1 ) are presented in Figure 7, where a sinusoidal fit is performed (χ 2 /dof= 29.4/16). Subtracting the obtained sinusoid from the data, we produce the periodogram shown in Figure 6 (lower panel, yellow line). The peak associated with the orbital frequency has been removed as expected. Also the satellite peaks are reduced, but the fact that some of the other peaks do not achieve a level consistent with the background test indicates that the signal in the LS I +61 • 303 data is not purely sinusoidal.
To reduce these remaining powers, we fitted the data set with a more complex signal. Motivated by the data shape, we fitted the data set with a sinusoidal function plus a Gaussian signal contribution (χ 2 /dof= 14.3/13), as shown in Figure 7 (blue dotted line). The corresponding periodogram subtracting this function to the LS I +61 • 303 data set is given in Figure 6 (lower panel, blue line). The orbital frequency peak has been removed and some of the periodogram peaks are much more reduced than in the purely sinusoidal subtraction, giving a better agreement with the background periodogram level.
We performed a MC simulation to evaluate the error in the frequency estimation without any signal shape assumption: we simulate light curves where the number of γ -ray and background candidates are selected randomly from Poisson distributions with a mean equal to the actually measured distributions of events, arriving in every given time interval. The periodogram is calculated for 10 3 of those randomly generated series, and the distribution of the resulting peak power frequencies is fitted with a Gaussian function, yielding an error of 0.0003 d −1 . An accurate peak frequency determination is done by scanning more frequencies (increasing the oversampling factor) around the frequency which has maximum probability in the periodogram, and is found to be (0.0373 ± 0.0003) d −1 , corresponding to a period of (26.8 ± 0.2) days.
In Figure 8, we show the period obtained with MAGIC data compared with the measurements in other wavelengths. The most accurate measure of the orbital period is (26.4960 ± 0.0028) days, reported in radio by Gregory (2002). We also show period measurements reported in near-IR and optical Vband (Paredes et al. 1994), optical wavelengths (Mendelson & Mazeh 1989), photometry in the I−B and I bands (Mendelson & Mazeh 1994), Hα measurements (Zamanov et al. 1999), and soft X-ray measurements from Wen et al. (2006) and Paredes et al. (1997). The period obtained with MAGIC data is compatible with the orbital radio measurement within 1.5σ .

CONCLUSION
We find that LS I +61 • 303 is a periodic γ -ray binary with an orbital period of 26.8 ± 0.2 days (and chance probability ∼ 10 −7 ), compatible with the optical, radio, and X-ray period. This result implies that the flux modulation is tied to the orbital period. The high state in VHE γ -rays occurs in the same phases as the X-ray high state. This is especially interesting since we found an additional hint for X-ray/γ -ray variability correlation in the orbital phase 0.85. A strictly simultaneous multiwavelength campaign is needed to investigate this correlation in more detail.
We looked for possible intranight variability and found the flux consistent with being constant within errors in 30-75 min timescales.
We produce energy spectra for two phase bins, 0.5 < φ < 0.6 and 0.6 < φ < 0.7, and averaged flux values for several phase bins. There is clear evidence for a significant change in the VHE γ -ray flux level between different phase bins of LS I +61 • 303. The spectral photon index does not show this dependence on the phase. All derived spectral photon indices are compatible with 2.6 ± 0.2, obtained from the most significant phase bin of LS I +61 • 303.
We can put constraints to the emission at the periastron passage and conclude that the system is detected in γ -rays only in the phases 0.5-0.9. Since significant emission is only detected in an orbital sector of the phases at which the maximum γ -ray flux should occur under photon-photon absorption (see Figure 5 in Dubus 2006), the latter can hardly be the only source of variability in the emission.
Thorough multiwavelength observations will allow us to probe the intrinsic variability of the nonthermal emission from LS I +61 • 303 along the orbit and can prove possible correlations between the X-ray and TeV energy bands. This is a necessary step for understanding the source nature, and the physics underlying the VHE radiation.