ar X iv :1 60 3. 06 97 3v 2 [a st ro -p h. H E ] 21 A pr 2 01 6 Astronomy & Astrophysicsmanuscript no. Super-orbital˙variability˙LSI˙Final˙Version-Publisher c© ESO 2016 April 22, 2016 Super-orbital variability of LS I +61◦303 at TeV energies M. L. Ahnen1, S. Ansoldi2, L. A. Antonelli3, P. Antoranz4, A. Babic5, B. Banerjee6, P. Bangale7, U. Barres de Almeida7,25, J. A. Barrio8, J. Becerra González9,26, W. Bednarek10, E. Bernardini11,27, B. Biasuzzi2, A. Biland1, O. Blanch12, S. Bonnefoy8, G. Bonnoli3, F. Borracci7, T. Bretz13,28, S. Buson14, A. Carosi3, A. Chatterjee6, R. Clavero9, P. Colin7, E. Colombo9, J. L. Contreras8, J. Cortina12, S. Covino3, P. Da Vela4, F. Dazzi7, A. De Angelis14, B. De Lotto2, E. de Oña Wilhelmi15, C. Delgado Mendez16, F. Di Pierro3, A. Domı́nguez8,29, D. Dominis Prester5, D. Dorner13, M. Doro14, S. Einecke17, D. Eisenacher Glawion13, D. Elsaesser13, A. Fernández-Barral12, D. Fidalgo8, M. V. Fonseca8, L. Font18, K. Frantzen17, C. Fruck7, D. Galindo19, R. J. Garcı́a López9, M. Garczarczyk11, D. Garrido Terrats18, M. Gaug18, P. Giammaria3, N. Godinović5, A. González Muñoz12,, D. Gora11, D. Guberman12, D. Hadasch20,⋆, A. Hahn7, Y. Hanabata20, M. Hayashida20, J. Herrera9, J. Hose7, D. Hrupec5, G. Hughes1, W. Idec10, K. Kodani20, Y. Konno20, H. Kubo20, J. Kushida20, A. La Barbera3, D. Lelas5, E. Lindfors21, S. Lombardi3, F. Longo2, M. López8, R. López-Coto12, A. López-Oramas12,34,⋆, P. Majumdar6, M. Makariev22, K. Mallot11, G. Maneva22, M. Manganaro9, K. Mannheim13, L. Maraschi3, B. Marcote19, M. Mariotti14, M. Martı́nez12, D. Mazin7,30, U. Menzel7, J. M. Miranda4, R. Mirzoyan7, A. Moralejo12, E. Moretti7, D. Nakajima20, V. Neustroev21, A. Niedzwiecki10, M. Nievas Rosillo8, K. Nilsson21,31, K. Nishijima20, K. Noda7, R. Orito20, A. Overkemping17, S. Paiano14, J. Palacio12, M. Palatiello2, D. Paneque7, R. Paoletti4, J. M. Paredes19, X. Paredes-Fortuny19, G. Pedaletti11, M. Persic2,32, J. Poutanen21, P. G. Prada Moroni23, E. Prandini1,33, I. Puljak5, W. Rhode17, M. Ribó19, J. Rico12, J. Rodriguez Garcia7, T. Saito20, K. Satalecka8, C. Schultz14, T. Schweizer7, S. N. Shore23, A. Sillanpää21, J. Sitarek10, I. Snidaric5, D. Sobczynska10, A. Stamerra3, T. Steinbring13, M. Strzys7, L. Takalo21, H. Takami20, F. Tavecchio3, P. Temnikov22, T. Terzić5, D. Tescaro9, M. Teshima7,30, J. Thaele17, D. F. Torres24,⋆, T. Toyama7, A. Treves2, V. Verguilov22, I. Vovk7, J. E. Ward12, M. Will 9, M. H. Wu15, R. Zanin19, (the MAGIC Collaboration), J. Casares9, and A. Herrero9 (Affiliations can be found after the references) Received ... , 2015; accepted ... , 2015 ABSTRACT Context. The gamma-ray binary LS I+61◦303 is a well-established source from centimeter radio up tovery high energy (VHE; E> 100 GeV). The broadband emission shows a periodicity of∼26.5 days, coincident with the orbital period. A longer (super-orbital) period of 1667± 8 days was proposed from radio variability and confirmed using optical and high-energy (HE; E ¿ 100 MeV) gamma-ray observations. In this paper, we report on a four-year campaign performed by MAGIC together with archival data concentrating on a search for a long-timescale signature in the VHE emission from LS I+61◦303. Aims. We focus on the search for super-orbital modulation of the VHE emission, similar to that observed at other energies, and on the search for correlations between TeV emission and an optical determination of the extension of the circumstellar disk. Methods. A four-year campaign has been carried out using the MAGIC telescopes. The source was observed during the orbital phases when the periodic VHE outbursts have occurred (φ = 0.55 – 0.75, one orbit= 26.496 days). Additionally, we included archival MAGIC observations and data published by the VERITAS collaboration in these studies. For the correlation studies, LS I+61◦303 has also been observed during the orbital phases where sporadic VHE emission had been detected in the past (φ = 0.75 – 1.0). These MAGIC observations were simultaneous with optical spectroscopy from the LIVERPOOL telescope. Results. The TeV flux of the periodical outburst in orbital phasesφ = 0.5 – 0.75 was found to show yearly variability consistent with the long-term modulation of∼4.5 years found in the radio band. This modulation of the TeV flux can be well described by a sine function with a best-fit period of 1610±58 days. The complete data, including archival observations, span two super-orbital periods. There is no evidence for acorrelation between the TeV emission and the mass-loss rate of the Be star, but this may be affected by the strong, short-timescale (as short as intra-day) variation displayed by the Hα fluxes. Key words. binaries: general – gamma rays: observations –gamma rays: binary – stars: individual (LS I+61◦303) 1. Introduction LS I +61◦303 (= V615 Cas ) is a gamma-ray binary com- posed of a rapidly rotating Be star of spectral type B0Ve ⋆ Corresponding authors: A. López-Oramas, e- mail: alopez@ifae.es, D. Hadasch, e-mail: hadasch@icrr.u-tokyo.ac.jp, D. F. Torres, e-mail: dtorres@ice.csic.es (Hutchings & Crampton 1981) with a circumstellar disk and a compact object of unknown nature. The compact object, ei- ther a neutron star (NS) or a stellar-mass black hole (BH), has an eccentric orbit (e = 0.54 ± 0.03) with a period of 26.4960(28) days, determined from radio observations (Gregory 2002). The periastron passage occurs at phaseφ = 0.23 ± 0.03, although its precise timing depends on the orbital solu- 1 http://arxiv.org/abs/1603.06973v2 M. L. Ahnen et al.: Multi-year study of LS I+61◦303 tion (Gregory 2002; Casares et al. 2005; Grundstrom et al. 2007; Aragona et al. 2009).1 LS I +61◦303 is one of the few binary systems detected from radio to VHE gamma rays. The origin of the broad- band emission is still under debate. LS I+61◦303 was pro- posed as a microquasar, based on its extended jet-like radio- emitting structures, see Massi et al. (2004), for instance,and more recently, Massi & Torricelli-Ciamponi (2014). However, images obtained with the VLBA during a complete orbital cy- cle showed a rotating tail-like morphology of overall size 5– 10 mas (Dhawan et al. 2006), consistent with a pulsar wind sce- nario (Maraschi & Treves 1981). Similar phase-resolved struc- tures were later observed by Albert et al. (2008d). The source ap- pears extended in X-rays (Paredes et al. 2007; Kargaltsev etal. 2014). Deep observations searching for X-ray pulsing have yielded only upper limits (Rea et al. 2010). LS I +61◦303 was first detected in the VHE regime by MAGIC in 2006 (Albert et al. 2006). During the follow- ing years, more observations have been performed by the MAGIC (Albert et al. 2008d; Anderhub et al. 2009; Albert et al. 2009; Aleksić et al. 2012a) and VERITAS collaborations (Acciari et al. 2008). The VHE emission region is unresolved(< 0.1 deg). The system is also a HE gamma-rayFermi-LAT source (Abdo et al. 2009). Long-term monitoring shows that the HE emission undergoes periodic outbursts slightly after the perias- tron passage, around phasesφ ∼ 0.3 – 0.45 (Hadasch et al. 2012). In the orbit from MJD 53752.7 until MJD 53779.2 (January - February 2006), a TeV peak was first detected at phasesφ ∼ 0.6 – 0.7 (which correspond to the phases next to the apastron) at a level of∼16% of the Crab Nebula flux for energies above 400 GeV. However, in October 2009 – January 2010, the source was in a low flux state (Aleksić et al. 2012a) and the emission peak was above 300 GeV with a flux of only 5.4% of the Crab Nebula. Despite the state of the system, the highest TeV flux also occurred at phasesφ = 0.6 – 0.7, consistent with previous observations. The spectral fit parameters agreed with our previ- ous determinations (Albert et al. 2008d; Anderhub et al. 2009; Albert et al. 2009). Significant emission has also occasionally been observed at phasesφ = 0.8 – 1, where MAGIC detected the system at∼4% Crab Units (C.U). (Albert et al. 2009). Emission around periastron has been reported once, by VERITAS, in late 2010 at phaseφ = 0.081 (Acciari et al. 2011) when the compact star was near superior conjunction. Orbit-to-orbit variability has been associated with a super- orbital modulation. This was first proposed based on centimeter radio variations that show approximately sinusoidal modulation over 1667± 8 days (Paredes 1987; Gregory 2002). A similar long-term behavior has recently been suggested for X-rays (3 – 30 keV observed withRXTE; Li et al. (2012); Chernyakova et al. (2012)), hard X-rays (18 – 60 keV observed withINTEGRAL; Li et al. (2014)) and HE gamma-rays (100 MeV – 300 GeV ob- served with Fermi/LAT; Ackermann et al. (2013)). In the LAT energy range, the super-orbital variability is almost invisible around the periastron, where the compact object is inside (or highly affected by) the Be circumstellar disk, but appears around apastron. Two short (< 0.1 s) and very luminous (> 1037 erg s−1) bursts have been detected from the direction of LS I+61◦303 by the Swift Burst Alert Telescope (Barthelmy et al. 2008; Burrows et al. 2012). A detailed analysis of these events waspre- sented by Torres et al. (2012), where their similarities to those 1 Here, the reference for the periastron passage is taken atT0 = 43366.275 MJD (Gregory 2002). observed from magnetars were analyzed. If the bursts are asso- ciated with the system, as it seems, it requires that the system harbors a neutron star and not a black hole. The idea of a transitioning system, that is, a system that swings from one state to another at every orbital period, has been suggested by Zamanov et al. (2001). Torres et al. (2012) and Papitto et al. (2012) extended this to the super-orbitalbe- havior and analyzed orbitally induced variability. In this so-called flip-flop scenario, the system changes from propeller regime accretion at periastron, where the pulsarmag- netosphere is disrupted, to an ejector regime (rotational-powered pulsar) at apastron, when particles are accelerated to TeV ener- gies in the inter-wind shock formed at the collision region be- tween the neutron star and the stellar wind. Despite the system differences, this would make LS I+61◦303 similar to the transi- tional pulsars in redbacks, for example, Archibald et al. (2009) and Papitto et al. (2013). Changing the Be star mass-loss rate can cause the switching from one state to the other to vary in orbital phase, modulated along the super-orbital timescale as measured by the equivalent width (EW) of the Hα line, for instance. The sizes of the stel- lar disks of Be stars correlate with the EW of the Hα emission line (e.g., see Porter & Rivinius (2003) and Reig (2011) for a review). Grundstrom & Gies (2006) also showed that the EW is correlated with the radius of the circumbinary disk, and therefore it can be used as a proxy of the latter. For times of high mass-loss rate, at which the influence of the disk matter upon the compact object can be felt along a greater portion of the orbit, the pro- peller regime can be operative even close to apastron. If this is the case, the TeV emission of LS I+61◦303 would be reduced (Torres et al. 2012). 2. Observations A campaign using the MAGIC telescopes together with the op- tical telescope was carried out over four years. Both instruments are located on the island of La Palma, in the Canary Islands, Spain, at the observatory of El Roque de Los Muchachos (28◦N, 18◦W, 2200 m above the sea level). The MAGIC telescopes are two Imaging Atmospheric Cherenkov Telescopes (IACTs) with diameters of 17 m, each one with a pixelized camera, cov- ering a field of view of∼ 3.5◦. The current sensitivity of the stereoscopic MAGIC telescopes is 0.71%± 0.02% of the Crab Nebula flux in 50 h of observation for energies above 250 GeV (Aleksic et al. 2015). The angular resolution at these energies is . 0.1◦ (68% containment radius) and the energy resolution is ∼ 18%. For monoscopic observations (also referred to as mono- observations) the integral sensitivity above 280 GeV is about 1.6% of the Crab Nebula flux in 50 hours (Aliu et al. 2009). The observations were carried out in wobble mode, pointing at two different symmetric regions situated 0.4◦ away from the source position to simultaneously evaluate the background (Fominet al. 1994). The data were analyzed using the standard MAGIC anal- ysis and reconstruction software, MARS (Zanin et al. 2013). The recorded images were calibrated, cleaned, and parametrized (Hillas 1985; Albert et al. 2008a). The background rejection and the estimation of the gamma-direction was performed us- ing the Random Forest (RF) method (Albert et al. 2008b). The energy of each event was estimated using look-up tables gen- erated by Monte Carlo simulations (Aleksić et al. 2012b). For mono-observations, the event direction and energy of the pri- mary gamma ray were also reconstructed with the RF method. 2 M. L. Ahnen et al.: Multi-year study of LS I+61◦303 Table 1. Observations of LS I+61◦303 performed by MAGIC. The first column denotes the orbit number, the second column indicates the dates of the observation, the third and fourthcol- umn display the orbital and super-orbital phases, while thefifth and sixth columns indicate the integration time and its distribu- tion, respectively. Orbit MJD φorbital φsuper−orbital Time Number number range range hours of days 1 55415.2 0.75 0.23 1.14 1 2 55441.2 - 55444.2 0.73 - 0.84 0.25 3.98 3 3 55471.1 0.86 0.27 0.76 1 4 55486.1 - 55500.1 0.42 - 0.95 0.28 3.63 4 5 55512.0 0.40 0.29 1.92 1 6 55543.0 0.57 0.30 2.06 1 7 55568.9 - 55574.0 0.55 - 0.74 0.32 10.81 6 21 55944.0 - 55945.0 0.70 - 0.74 0.55 2.56 6 22 55969.8 - 55977.8 0.68 - 0.99 0.56 3.91 6 32 56242.0 - 56243.0 0.95 - 0.99 0.72 2.20 2 33 56266.9 - 56267.9 0.89 - 0.93 0.74 2.10 2 34 56295.9 - 56296.8 0.99 - 0.01 0.77 4.04 2 44 56549.1 - 56550.1 0.54 - 0.58 0.91 5.67 2 45 56576.1 - 56579.1 0.56 - 0.67 0.92 7.90 4 46 56602.0 - 56607.1 0.54 - 0.73 0.94 9.90 5 48 56656.9 - 56663.9 0.61 - 0.87 0.98 15.65 8 57 56900.1 0.79 0.12 2.22 1 58 56920.1 - 56930.1 0.54 - 0.92 0.13 20.72 10 The LIVERPOOL robotic telescope is an optical telescope that is also located on La Palma. Its instrument, FRODOspec, provides spectra with R∼ 5500 resolution simultaneously in the blue and red spectral ranges. The red spectrum includes the Hα line, which we used as the prime indicator of the Be circumstel- lar disk through the equivalent width (EW), full width at half maximum (FWHM), and centroid velocity, the latter two ob- tained through a single Gaussian fit to the emission profile. LS I +61◦303 was observed between August 2010 and September 2014. All data were obtained with the MAGIC stereoscopic system, except for January 2012, when MAGIC- I was inoperative. The source was observed during the orbital rangeφ = 0.5 – 0.75 to observe the complete trend of the pe- riodical outburst peak of the TeV emission, with the aim of detecting a putative long-term modulation. Contemporaneous observations with MAGIC and LIVERPOOL were performed during orbital phases 0.75 – 1.0, which are the phases where sporadic VHE emission had been detected and which does not seem to present yearly periodical variability of the flux level (Aleksić et al. 2012a). Since the fluxes in this orbital period are not influenced by the long-term modulation, changes in the rel- ative optical and TeV fluxes are larger and easier to measure. The aim of these contemporaneous observations is to search for (anti-)correlation between the mass-loss rate of the Be star and the TeV emission. The details of the observations in an orbit-to- orbit basis are summarized in Table 1. 3. Results 3.1. Spectral stability LS I +61◦303 has shown variability on timescales of years in the strength of its periodic outburst peaks. To understand if the mechanism producing gamma rays is the same, independent of the flux of the source and its super-orbital state, it is interesting to search for spectral variability as a function of time. TheVHE Table 2. Daily integrated flux for energies above 300 GeV of LS I +61◦303 measured by MAGIC from 2010 to 2014. The nights with an asterisk are those where simultaneous optical data taken by LIVERPOOL is available. Horizontal lines indicatedif- ferent orbits. MJD φ Significance Integral flux Timee f f (Li&Ma) (E ¿ 300 GeV) [Days] [σ] [10−12cm−2s−1] [hours] 55415.2* 0.75 1.5 3.3± 1.6 1.14 55441.2 0.73 6.4 5.8± 1.3 1.94 55442.1* 0.76 1.8 4.9± 1.7 1.19 55444.2* 0.84 -0.7 -1.2± 0.9 0.85 55471.1* 0.86 0.6 -0.6± 1.2 0.76 55486.1 0.42 3.4 3.9± 1.4 1.12 55498.1* 0.87 0.6 0.4± 0.7 1.24 55499.1 0.91 -0.3 -0.5± 1.4 0.38 55500.1 0.95 0.9 0.2± 1.4 0.89 55512.0 0.40 1.6 1.5± 1.2 1.92 55543.0 0.57 2.3 3.1± 1.0 2.06 55568.9 0.55 12.3 9.4± 1.5 1.87 55569.9 0.59 2.4 2.3± 1.1 1.61 55571.0 0.62 7.0 5.8± 1.1 2.62 55572.0 0.66 6.0 10.0± 2.3 1.11 55573.0 0.70 2.9 1.7± 1.2 1.13 55574.0 0.74 2.2 0.4± 0.3 2.47 55944.0 0.70 4.0 12.0± 3.2 1.70 55945.0 0.74 0.7 1.0± 4.2 0.86 55969.8 0.68 -0.4 -1.1± 1.1 0.92 55970.8 0.72 1.7 2.5± 1.9 0.35 55975.8 0.91 1.3 1.3± 0.8 1.66 55976.8 0.95 1.1 1.1± 1.9 0.42 55977.8* 0.99 1.3 1.3± 1.9 0.57 56242.0 0.95 0.7 2.6± 1.6 0.99 56243.0 0.99 1.7 4.5± 1.8 1.24 56266.9* 0.89 2.3 6.6± 1.9 1.09 56267.9* 0.93 5.0 7.0± 2.0 1.02 56295.9* 0.99 3.0 4.1± 1.4 1.71 56296.8* 0.01 -0.6 -1.4± 0.9 2.33 56549.1 0.54 3.1 6.1± 1.6 2.69 56550.1 0.58 0.6 0.3± 1.0 2.98 56576.1 0.56 0.2 0.8± 1.8 1.09 56577.0 0.60 -2.7 -1.8± 1.1 2.82 56578.0 0.63 -2.7 -1.8± 1.0 2.69 56579.1 0.67 0.2 0.4± 1.5 1.30 56602.0 0.54 0.7 1.6± 1.3 2.10 56603.0 0.58 2.6 1.3± 1.0 2.93 56604.0 0.61 1.0 2.3± 1.2 2.72 56606.0 0.69 1.6 1.8± 1.6 0.98 56607.1 0.73 3.6 5.3± 1.9 1.17 56656.9 0.61 0.9 1.2± 1.5 1.68 56657.9 0.65 -0.1 -0.2± 1.1 2.93 56658.9 0.69 2.0 2.7± 1.5 2.37 56659.9 0.72 3.0 4.1± 1.6 1.75 56660.9* 0.76 1.3 1.7± 1.4 2.20 56661.9* 0.80 0.1 0.1± 1.5 1.14 56662.9* 0.84 0.7 1.0± 1.5 1.95 56663.9* 0.87 4.6 5.0± 1.4 1.63 56900.1 0.78 -0.21 0.9± 1.1 2.22 56920.1 0.54 4.92 4.6± 1.4 1.74 56921.1 0.58 9.27 12.1± 1.6 2.27 56922.1 0.61 1.43 2.2± 1.0 2.93 56923.1 0.65 6.48 7.2± 1.2 2.92 56924.1 0.69 0.88 0.5± 0.9 2.93 56925.1 0.73 0.21 1.7± 1.3 1.49 56927.1 0.80 0.12 1.3± 1.9 0.62 56928.1 0.84 -0.83 -0.5± 1.0 1.94 56929.1 0.88 -1.49 -0.6± 0.9 1.93 56930.1 0.92 -1.42 -1.8± 0.9 1.94 spectrum for the complete observed data set (see Table 2) canbe described as dNγ dAdtdE = N0 ( E E0 )α , (1) with N0 = (4.4 ± 0.1stat ± 0.2sys) × 10−13TeV−1cm−2s−1, α = −2.4± 0.2stat ± 0.2sys and the normalization energy E0 = 1 TeV. 3 M. L. Ahnen et al.: Multi-year study of LS I+61◦303 Table 3. Spectral results for the different observational cam- paigns of LS I+61◦303. The first column indicates the cam- paign name, the second the orbital range where the SED was computed, the third the average super-orbital phase of the cam- paign, and the fourth the spectral index with its corresponding statistical and systematic uncertainties. Campaign Orbital Super-orbital Spectral interval phase index (Albert et al. 2006) 0.4 – 0.7 ∼0.22 −2.6± 0.2± 0.2 (Anderhub et al. 2009) 0.4 – 0.7 ∼0.59 −2.6± 0.2± 0.2 (Albert et al. 2009) 0.6 – 0.7 ∼0.41 −2.6± 0.2± 0.2 (Aleksić et al. 2012a) 0.6 – 0.8 ∼0.08 −2.5± 0.5± 0.2 2011 (this work) 0.6 – 0.8 ∼0.28 −2.2± 0.1± 0.2 2012 (this work) 0.7 ∼0.55 −2.7± 0.5± 0.2 2013 (this work) 0.5 – 0.8 ∼ 0.92 −2.5± 0.5± 0.2 2014 (this work) 0.5 – 0.8 ∼ 0.13 −2.5± 0.1± 0.2 Table 3 collects the spectral indices for each of the observa- tional campaigns. They are all similar within systematic errors. In 2011, a hardening was observed, but with low statistical sig- nificance. It was not confirmed by later observations Therefore, we can conclude that the VHE spectrum is consistent with a sin- gle power-law during different epochs. The dependence ofα on phase with the 1667 day super- orbital period is shown in Fig. 1. The data can be fit to a constant value of 2.43± 0.04 (χ2/dof= 8.9/7). We divided the dataset into intervals of high (defined as the flux being at 5 – 10% of the Crab Nebula flux) or low (flux at 2 – 5% of the Crab Nebula flux) activity. We also separated the data depending on the orbital interval of the emission into a period- ical (φ = 0.5 – 0.75) interval, where the periodic TeV outburst occurs, and a sporadic (φ = 0.75 – 1.0) interval, where signifi- cant emission has been detected sporadically in the past. Innei- ther case does the spectral index vary significantly. Fitting the data with a constantα results in an average value of 2.4 ± 0.1 (χ2/do f = 0.90/2) with a probability of 0.83. Superorbital Phase 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 α 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 Albert et al. 2006 Anderhub et al. 2009 Albert et al. 2009 Aleksic et al. 2012a 2011 (this work) 2012 (this work) 2013 (this work) 2014 (this work) Fig. 1. Super-orbital dependence of the spectral indexα for all MAGIC campaigns of LS I+61◦303, considering a 1667-day period. The blue line corresponds to the average value. 3.2. Super-orbital modulation We performed a search for indications of a super-orbital flux modulation using the peak flux measurements at orbital phases φ = 0.5 – 0.75. However, some of the system’s emission peaks may have been missed, for instance, because observations can only be performed once every 24 hours and because some nights were lost because of bad weather or technical problems. To min- imize this effect, we only selected those measurements for which at least three consecutive observations in the requisite phase in- terval were successful and for which the middle time shows a higher flux than the others. We justified this approach by simulating exactly this situa- tion. Since observations can only be performed once every 24 hours, a maximum in emission from the binary system may be missed. Additionally, observations were also lost during inter- vals of bad weather or technical problems. To minimize the ef- fects of missing the maximum emission during any cycle, we only considered those orbits for which observations occurred during at least three consecutive nights in the specific phase in- terval. We also required that the middle observation showedthe highest flux. By further assuming a shape (Landau distribution) similar to the one observed in Anderhub et al. (2009), we then calculated the percentage deviation between the actual maxi- mum of the emission and the assigned value of the maximum. We checked that at one sigma, the difference between the esti- mated maximum flux level and the actual maximum was smaller than< 35%, while this value increased to 90% if the maximum was used, even requiring three consecutive nights. This estimate allowed us to study the dependence of the amplitude on the or- bital periodic emission observed at VHE with time. All archival data of LS I+61◦303 recorded by MAGIC since its detection in 2006 (Albert et al. 2008d; Anderhub et al. 2009; Albert et al. 2009; Aleksić et al. 2012a) and the data from the observing campaigns presented were folded onto the superor- bital period of 1667 days. We considered statistical and sys- tematic uncertainties in the integrated flux: 12% systematics for stereo data, according to Aleksić et al. (2012b), and 15% sys- tematics for mono data, estimated from Albert et al. (2008c). We augmented our sample using the many observations made by VERITAS of LS I +61◦303, that is, Acciari et al. (2008, 2009, 2011) and Aliu et al. (2013), using the same procedure to iden- tify the peak of emission. We assumed 20% systematic uncer- tainties (Griffin 2014) that were added quadratically to the cor- responding statistical uncertainties. The temporal evolution of the peak-integrated flux above 300 GeV for each orbit is presented in Fig. 2. Folded light-curve data were fit by first assuming a hypothesis of a constant flux, and secondly assuming a sinusoidal function. The fit to a sine function yielded a period of 1610± 58 days with a probability of 6% (χ2/dof= 26.6/17). Since the fit probability for a constant fit is much lower, with 4.5× 10−12 (χ2/dof = 114.8/20), we can conclude that the TeV flux is variable along the superorbit and that a sinusoidal behavior with a period of 1610 days is preferred over a constant one. Figure 3 shows the same data, but folded onto the super- orbital period of 1667 days. Here the data were fit not only with a constant and a sinusoidal, but also to two-emission levels(step function): f (x) = { a if p1 < x < p2, b if x < p1 & x > p2, wherea, b, p1 , andp2 are the free parameters of the fit. The results of the different fits are presented in Table 4. The prob- 4 M. L. Ahnen et al.: Multi-year study of LS I+61◦303 MJD 54000 54500 55000 55500 56000 56500 57000 ] -1 s -2 c m -1 2 In te gr al F lu x P ea k (E > 30 0 G eV ) [1 0 10− 5− 0 5 10 15 20 25 30 35 40 MAGIC VERITAS Constant Sinusoidal Fig. 2. Temporal evolution of the peak of the VHE emission for each orbital period in terms of the MJD. MAGIC (magenta dots) and VERITAS (blue squares) data within orbital phases 0.5 – 0.75 have been considered. The gray dashed line represents 10% of the Crab Nebula flux, the gray solid line the zero level. Thefit with a sinusoidal is plotted in red and the fit with a constant in blue. Superorbital Phase 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 ] -1 s -2 c m -1 2 In te gr al F lu x P ea k (E > 30 0 G eV ) [1 0 10− 5− 0 5 10 15 20 25 30 35 40 MAGIC VERITAS Constant Sinusoidal Step function Fig. 3. Peak of the VHE emission in terms of the super-orbital phase defined in radio. Each data point represents the peak flux emitted in one orbital period during orbital phases 0.5 – 0.75 and is folded into the super-orbit of 1667 days known from ra- dio observations (Gregory 2002). MAGIC (magenta dots) and VERITAS (blue squares) points have been used in this analysis. The fit with a sinusoidal (solid red line), with a step function (solid green line), and with a constant (solid blue line) arealso represented. The gray dashed line represents 10% of the Crab Nebula flux, the gray solid line the zero level for reference. ability for a constant flux is negligible, 4.5×10−12. Assuming a sinusoidal signal, the fit probability reaches 8% (χ2/dof = 27.2/18). The fit to a step function resulted in a fit probability of 7% (χ2/dof = 26.4/17). We furthermore quantified the prob- ability that the improvement found when fitting a sinusoid ora step function instead of a constant is produced by chance. To obtain this probability, we considered the likelihood ratio test (Mattox et al. 1996). In both cases this chance probability is < 2.5× 10−10, which is low. This shows that the observed inten- sity distribution can be described by a high and a low state and Table 4.Fitting functions with the corresponding fit probabilities for MAGIC + VERITAS data of LS I+61◦303 folded into the superorbit (Fig. 3). Function Fit probability χ2/dof Constant 4.5× 10−12 114.8/20 Step function 0.07 26.4/17 Sinusoidal 0.08 27.2/18 Table 5.Correlations between the TeV flux obtained by MAGIC and theHα parameters (EW, FWHM, and vel) measured by LIVERPOOL for the extended orbital interval 0.75 – 1.0. Only TeV data with a significance higher than 1σ have been consid- ered. The first column indicates the level of simultaneity ofthe observations, the second column shows the parameters we used to search for a correlation, and the third and fourth columnsgive the Pearson correlation coefficient and the associated probability for a non-correlation. Simultaneity Parameters r Prob Nightly TeV - EW -0.51 0.04 Nightly TeV - FWHM -0.22 0.27 Nightly TeV - vel -0.38 0.11 Contemporary (hourly) TeV - EW -0.14 0.37 Contemporary (hourly) TeV - FWHM -0.44 0.16 Contemporary (hourly) TeV - vel -0.21 0.35 with a smoother transition. We conclude that there is a super- orbital signature in the TeV emission of LS I+61◦303 and that it is compatible with the 4.5-year radio modulation seen in other frequencies. 3.3. Simultaneous optical-TeV observations The correlation between the TeV flux measured by MAGIC and the Hα parameters measured with the telescope (EW, FWHM, profile centroid velocity) were determined including statistical and systematic uncertainties and the weighted Pearson correla- tion coefficient (Soper et al. 1917). To search for correlations, we classified the data into three different categories: – simultaneous: The optical observations were performed pre- cisely during the period when MAGIC was performing its observations. Only three such points were obtained under this condition. Because of the scarcity of such data we do not provide any correlation coefficient for this group. – contemporary (hourly): Data have a time difference of three hours at most. – nightly: Data were obtained during the same night. Only TeV data points with a significance higher than 1σ were considered. The data used for the correlation are given in Table 5. No statistically significant correlation was found for the sample at orbital phasesφ = 0.75 – 1.0. A hint of a correlation is ob- served, but its significance is low. A stronger correlation might be blurred as a result of the fast variability of the optical param- eters on short timescales compared to the long exposure times required by MAGIC, and as a result of the relatively large uncer- tainties and small number of data points used for this analysis. 5 M. L. Ahnen et al.: Multi-year study of LS I+61◦303 0 2 4 6 8 10 E W [A ] 8 9 10 11 12 13 14 15 16 17 Nightly Contemporary Simultaneous 0 2 4 6 8 10 F W H M [k m /s ] 650 700 750 800 850 ]-1 s-2 cm-12TeV flux (E>300GeV) [10 0 2 4 6 8 10 ve l H al ph a [k m /s ] 20− 10− 0 10 20 30 40 Fig. 4. Correlations between the TeV flux obtained by MAGIC and the Hα parameters (from top to bottom: EW, FWHM, and centroid velocity (vel)) measured by the LIVERPOOL Robotic Telescope for the orbital interval 0.75 – 1.0. Only TeV data points with a significance higher than 1σ have been con- sidered. Each data point represents a ten-minute observation in the optical and a variable integration in the TeV regime: nightly (blue), contemporary (red), and strictly simultaneous data (green). Black error bars represent statistical uncertainties, while systematic uncertainties are plotted as magenta error bars. Figure 4 shows the Hα measurements plotted against the TeV flux. 4. Discussion and conclusion The main conclusions from this multi-year analysis of LS I +61◦303 TeV observations are listed below. 1. We achieved a first detection of super-orbital variability in the TeV regime. Using the new VHE data and the MAGIC and VERITAS archival data, we found that the super-orbital signature of LS I+61◦303 is consistent with the 1667-day ra- dio period within 8%. The TeV period obtained when fitting the data to a sine function with a free frequency instead of performing the fit as before with a fixed super-orbital period known from radio, is 1610± 58 days with a fit probability of 6%. 2. There is no statistically significant intra-day correlation be- tween Hα line properties and TeV emission, nor is there an obvious trend connecting the two frequencies. The flip-flop model (Zamanov et al. 2001; Torres et al. 2012; Papitto et al. 2012) considers LS I+61◦303 as a pulsar-Be star binary that changes accretion states from a propeller regime dur- ing periastron to an ejector regime at apastron. The change of state is thought to be driven by the influence of matter. The higher the pressure of this matter on the (hypothesized) neutron star magnetosphere, the more likely the latter will be compressed and disrupted. When this occurs, the pulsar wind is affected or disappears, and the inter-wind shock, which is thought to con- tribute to the multi-frequency non-thermal emission, disappears as well. Even when the neutron star passage at periastron is able to cut the disk off, mass accumulation will be higher in periods in which the decretion rate is higher, reaching farther out in the orbit. A larger EW of Hα may be assumed as a proxy for such a situation. If this is the case, the optical emission should be anti- correlated with the TeV flux. The detection of TeV gamma-ray super-orbital emission conforms to the predicted long-term be- havior of the flip-flopping model (Torres et al. 2012). The source was found in high and low states, when expected. This result ex- tends the earlier indications for this phenomenology foundusing smaller samples of TeV data (Li et al. 2012). However, we do not see this (anti-) correlation in the intra-day scales we tested. It may be that the EW of Hα is not the best tracer for decre- tion disk size or mass accumulation, or simply that the fast and extreme changes (by up to a factor of several in the same day) in the Hα data of the source and the vastly different integration times (minutes as compared with several hours) needed to claim a detection in both frequencies prevent us from measuring any possible trend. Acknowledgements. We would like to thank the Instituto de Astrofı́sica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The support of the German BMBF and MPG, the Italian INFN, the Swiss National Fund SNF, and the ERDF funds under the Spanish MINECO is gratefully acknowledged. This work was also supported by the CPAN CSD2007-00042 and MultiDark CSD2009-00064 projects of the Spanish Consolider-Ingenio 2010 programme, by grant 127740 of the Academy of Finland, by the Croatian Science Foundation (HrZZ) Project 09/176, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, and by the Polish MNiSzW grant 745/N-HESS-MAGIC/2010/0. J.C., and D.F.T. ac- knowledge support by the Spanish Ministerio de Economı́a y Competividad (MINECO) under grant AYA2010-18080, and by MINECO and the Generalitat de Catalunya under grants AYA2012-39303 and SGR 2014-1073,respectively. References Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 701,L123 Acciari, V. A., Aliu, E., Arlen, T., et al. 2011, ApJ, 738, 3 Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, ApJ, 700, 1034 Acciari, V. A., Beilicke, M., Blaylock, G., et al. 2008, ApJ,679, 1427 Ackermann, M., Ajello, M., Ballet, J., et al. 2013, ApJ, 773,L35 Albert, J., Aliu, E., Anderhub, H., et al. 2009, ApJ, 693, 303 Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771 Albert, J., Aliu, E., Anderhub, H., et al. 2008a, Nuclear Instruments and Methods in Physics Research A, 594, 407 Albert, J., Aliu, E., Anderhub, H., et al. 2008b, Nuclear Instruments and Methods in Physics Research A, 588, 424 Albert, J., Aliu, E., Anderhub, H., et al. 2008c, ApJ, 674, 1037 Albert, J., Aliu, E., Anderhub, H., et al. 2008d, ApJ, 684, 1351 Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012a,ApJ, 746, 80 Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012b,Astroparticle Physics, 35, 435 Aleksic, J., Ansoldi, S., Antonelli, L. A., et al. 2015, ApJ,927, 6505 Aliu, E., Anderhub, H., Antonelli, L. A., et al. 2009, Astroparticle Physics, 30, 293 Aliu, E., Archambault, S., Behera, B., et al. 2013, ApJ, 779,88 6 M. L. Ahnen et al.: Multi-year study of LS I+61◦303 Anderhub, H., Antonelli, L. A., Antoranz, P., et al. 2009, ApJ, 706, L27 Aragona, C., McSwain, M. V., Grundstrom, E. D., et al. 2009, ApJ, 698, 514 Archibald, A. M., Stairs, I. H., Ransom, S. M., et al. 2009, Science, 324, 1411 Barthelmy, S. D., Baumgartner, W., Cummings, J., et al. 2008, GRB Coordinates Network, 8215, 1 Burrows, D. N., Chester, M. M., D’Elia, V., et al. 2012, GRB Coordinates Network, 12914, 1 Casares, J., Ribas, I., Paredes, J. M., Martı́, J., & AllendePrieto, C. 2005, MNRAS, 360, 1105 Chernyakova, M., Neronov, A., Molkov, S., et al. 2012, ApJ, 747, L29 Dhawan, V., Mioduszewski, A., & Rupen, M. 2006, in VI Microquasar Workshop: Microquasars and Beyond Fomin, V. P., Stepanian, A. A., Lamb, R. C., et al. 1994, Astroparticle Physics, 2, 137 Gregory, P. C. 2002, ApJ, 575, 427 Griffin, S. 2014, in CAP Congress Grundstrom, E. D., Caballero-Nieves, S. M., Gies, D. R., et al. 2007, ApJ, 656, 437 Grundstrom, E. D. & Gies, D. R. 2006, ApJ, 651, L53 Hadasch, D., Torres, D. F., Tanaka, T., et al. 2012, ApJ, 749,54 Hillas, A. M. 1985, International Cosmic Ray Conference, 3,445 Hutchings, J. B. & Crampton, D. 1981, PASP, 93, 486 Kargaltsev, O., Rangelov, B., Hare, J., & Pavlov, G. G. 2014,Astronomische Nachrichten, 335, 301 Li, J., Torres, D. F., & Zhang, S. 2014, ApJ, 785, L19 Li, J., Torres, D. F., Zhang, S., et al. 2012, ApJ, 744, L13 Maraschi, L. & Treves, A. 1981, MNRAS, 194, 1P Massi, M., Ribó, M., Paredes, J. M., et al. 2004, A&A, 414, L1 Massi, M. & Torricelli-Ciamponi, G. 2014, A&A, 564, A23 Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396 Papitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517 Papitto, A., Torres, D. F., & Rea, N. 2012, ApJ, 756, 188 Paredes, J. M. 1987, PhD thesis, , University of Barcelona, (1987) Paredes, J. M., Ribó, M., Bosch-Ramon, V., et al. 2007, ApJ,664, L39 Porter, J. M. & Rivinius, T. 2003, PASP, 115, 1153 Rea, N., Torres, D. F., van der Klis, M., et al. 2010, MNRAS, 405, 2206 Reig, P. 2011, Ap&SS, 332, 1 Soper, H. E., Young, A. W., Cave, B. M., Lee, A., & Pearson, K. 1917, Biometrika Torres, D. F., Rea, N., Esposito, P., et al. 2012, ApJ, 744, 106 Zamanov, R., Marti, J., & Marziani, P. 2001, in The Second National Conference on Astrophysics of Compact Objects, 50 Zanin, R., Carmona, E., & Sitarek, J. 2013, in Proceedings ofthe ICRC 2013, Internatinal Cosmic Ray Conference 1 ETH Zurich, CH-8093 Zurich, Switzerland 2 Università di Udine, and INFN Trieste, I-33100 Udine, Italy 3 INAF National Institute for Astrophysics, I-00136 Rome, Italy 4 Università di Siena, and INFN Pisa, I-53100 Siena, Italy 5 Croatian MAGIC Consortium, Rudjer Boskovic Institute, University of Rijeka, University of Split and University ofZagreb, Croatia 6 Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Salt Lake, Sector-1, Kolkata 700064, India 7 Max-Planck-Institut für Physik, D-80805 München, Germany 8 Universidad Complutense, E-28040 Madrid, Spain 9 Inst. de Astrofı́sica de Canarias, E-38200 La Laguna, Tenerife, Spain; Universidad de La Laguna, Dpto. Astrofı́sica, E-38206 La Laguna, Tenerife, Spain 10 University of Łódź, PL-90236 Lodz, Poland 11 Deutsches Elektronen-Synchrotron (DESY), D-15738 Zeuthen, Germany 12 Institut de Fisica d’Altes Energies, The Barcelona Institute of Science and Technology (IFAE-BIST), Campus UAB, 08193 Bellaterra (Barcelona), Spain 13 Universität Würzburg, D-97074 Würzburg, Germany 14 Università di Padova and INFN, I-35131 Padova, Italy 15 Institute for Space Sciences (CSIC/IEEC), E-08193 Barcelona, Spain 16 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas, E-28040 Madrid, Spain 17 Technische Universität Dortmund, D-44221 Dortmund, Germany 18 Unitat de Fı́sica de les Radiacions, Departament de Fı́sica, and CERES-IEEC, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Spain 19 Universitat de Barcelona, ICC, IEEC-UB, E-08028 Barcelona, Spain 20 Japanese MAGIC Consortium, ICRR, The University of Tokyo, Department of Physics and Hakubi Center, Kyoto University,Tokai University, The University of Tokushima, KEK, Japan 21 Finnish MAGIC Consortium, Tuorla Observatory, Universityof Turku and Department of Physics, University of Oulu, Finland 22 Inst. for Nucl. Research and Nucl. Energy, BG-1784 Sofia, Bulgaria 23 Università di Pisa, and INFN Pisa, I-56126 Pisa, Italy 24 ICREA and Institute for Space Sciences (CSIC/IEEC), E-08193 Barcelona, Spain 25 now at Centro Brasileiro de Pesquisas Fı́sicas (CBPF/MCTI), R. Dr. Xavier Sigaud, 150 - Urca, Rio de Janeiro - RJ, 22290-180, Brazil 26 now at NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA and Department of Physics and Department of Astronomy, University of Maryland, College Park, MD 20742, USA 27 Humboldt University of Berlin, Institut für Physik Newtonstr. 15, 12489 Berlin Germany 28 now at Ecole polytechnique fédérale de Lausanne (EPFL), Lausanne, Switzerland 29 now at Department of Physics & Astronomy, UC Riverside, CA 92521, USA 30 also at Japanese MAGIC Consortium 31 now at Finnish Centre for Astronomy with ESO (FINCA), Turku, Finland 32 also at INAF-Trieste 33 also at ISDC - Science Data Center for Astrophysics, 1290, Versoix (Geneva) 34 now at Laboratoire AIM, Service d’Astrophysique, DSM\IRFU, CEA\Saclay FR-91191 Gif-sur-Yvette Cedex, France 7 1 Introduction 2 Observations 3 Results 3.1 Spectral stability 3.2 Super-orbital modulation 3.3 Simultaneous optical-TeV observations 4 Discussion and conclusion