FERMI LAT OBSERVATIONS OF LS I +61°303: FIRST DETECTION OF AN ORBITAL MODULATION IN GeV GAMMA RAYS

The high-mass X-ray binary LS I +61°303 (=V615 Cas) has long been plausibly associated with a high-energy (HE, 20 MeV–100 GeV) gamma-ray source, although never before confirmed. The discovery of the COS B source 2CG 135+01 (Hermsen et al. 1977) quickly brought attention to this binary system's Be star localized within its error box, because of its unusual periodic radio emission (Gregory et al. 1979) and its X-ray emission (Bignami et al. 1981). 2CG 135+01 was to remain one of the brightest sources known in the HE gamma-ray sky, with a flux of ∼10⁻⁶ ph s⁻¹ cm⁻² above 100 MeV (Swanenburg et al. 1981). In the 1990s, EGRET detected the source with high confidence at the same average flux level and derived a power-law photon index of Γ = 2.05 ± 0.06 (Kniffen et al. 1997). Although there are no other objects of note (radio-loud active galactic nuclei, or pulsars) coinciding with the 3EG source (Hartman et al. 1999), its positional uncertainty was not small enough to firmly associate the gamma-ray source with the binary. Variability in the EGRET light curve could be neither firmly established nor related to variability at other wavelengths (Tavani et al. 1998; Nolan et al. 2003). Recently, AGILE has reported detecting the source at the same flux level (Pittori et al. 2009).

LS I +61°303 is an unusual binary system exhibiting strong variable emission from the radio to X-ray and TeV energies. At radio wavelengths, the source has been shown to exhibit radio outbursts that are modulated on an orbital period of 26.4960 ± 0.0028 days (Taylor & Gregory 1982; Gregory 2002). The phase of radio maximum has also been shown by Gregory (2002) to vary with a super-orbital period of 1667 ± 8 days. Observations of orbital modulation in the optical place constraints on the binary system parameters. The binary has an eccentric orbit (e = 0.55–0.72) and the Be star radial velocity is consistent with a neutron star companion or, if the orbital inclination is ⩽25°, with a ⩾3M_☉ black hole (Hutchings & Crampton 1981; Casares et al. 2005). Significant uncertainty still exists in key parameters of the orbital solution of the system (Grundstrom et al. 2007; Aragona et al. 2009).

Behavior in the X-ray band is much more complicated. Orbital modulation has been reported with the peak of emission appearing at phases 0.6–0.7 (Paredes et al. 1997; Esposito et al. 2007). However, the modulation is not smooth, with short-timescale flares and very strong orbit-to-orbit variability (Smith et al. 2009). Broadband spectral analysis of XMM-Newton and INTEGRAL data by Chernyakova et al. (2006) reveal LS I +61°303 to be well fitted by a simple absorbed power law with a hard photon index, Γ≃1.5, in the 0.5–100 keV band.

The MAGIC telescope detected a variable very high-energy (VHE, >100 GeV) gamma-ray source coincident with LS I +61°303 (Albert et al. 2006); a result that has been independently confirmed by the VERITAS collaboration (Acciari et al. 2008). More recently, the MAGIC collaboration has further reported that the VHE emission is periodic at the 26.5 day orbital period of the system (Albert et al. 2009). The VHE emission is consistently highest close to apastron, when the compact object is farthest from the Be star, and remains undetected at periastron. Like LS 5039 and PSR B1259 − 63 (Aharonian et al. 2005a, 2005b), and contrary to Cyg X-1 (Albert et al. 2007), LS I +61°303 is a gamma-ray binary with its spectral energy distribution (SED) peaking in HE gamma rays (for a full SED, see Dubus 2006; Chernyakova et al. 2006).

Calculations of the theoretical expectations of the gamma-ray emission from LS I +61°303 go back almost three decades (Maraschi & Treves 1981), and there has been a recent burst of activity following the MAGIC detection. Two scenarios have been put forward involving either the relativistic wind of a young, rotation-powered pulsar (Dubus 2006; Sierpowska-Bartosik & Torres 2008, 2009), or the relativistic jet of an accreting black hole or neutron star (Romero et al. 2005; Bednarek 2006; Gupta & Böttcher 2006; Bosch-Ramon et al. 2006). In light of the orbital modulations seen in radio, X-ray, and VHE gamma rays, a detailed light curve in the HE gamma-ray domain (where most of the energy is output) is an essential piece to identify the main radiative process at work and model the source.

The Fermi Gamma-ray Space Telescope was launched on 2008 June 11, from Cape Canaveral, Florida. The Large Area Telescope (LAT) is an electron–positron pair production telescope, featuring solid state silicon trackers and cesium iodide calorimeters, sensitive to photons from ∼20 MeV to >300 GeV (Atwood et al. 2009). Relative to earlier gamma-ray missions, the LAT has a large ∼2.4 sr field of view, a large effective area (∼8000 cm² for >1 GeV on axis) and improved angular resolution or point-spread function (PSF, better than 1° for 68% containment at 1 GeV). The Fermi survey mode operations began on 2008 August 4, after the conclusion of a flawless commissioning period. In this mode, the observatory is rocked north and south on alternate orbits to provide more uniform coverage so that every part of the sky is observed for ∼30 minutes every 3 hr. Thus, Fermi is ideally suited for long-term all-sky observations. The data set for this analysis spanned 2008 August 4 through 2009 March 24. Thus, LS I +61°303 was observed for approximately nine orbital periods.

The data were reduced and analyzed using the Fermi Science Tools v9r8 package.⁵⁸ The standard onboard filtering, event reconstruction, and classification were applied to the data (Atwood et al. 2009), and for this analysis the high-quality ("diffuse") event class is used. Time periods, when the region around LS I +61°303 was observed at a zenith angle greater than 105°, were also excluded to avoid contamination from Earth albedo photons. With these cuts, a photon count map of a 10° region around the binary is shown in Figure 1. The alignment of the LAT pointing direction with the celestial frame was calibrated using a large set of high-latitude gamma-ray sources to better than 10'' (Abdo et al. 2009b). The position of LS I +61°303 was found to be R.A. = 02^h40^m223, decl. = 61°13'30''(J2000) with a 95% error of 0069; in agreement with the accepted position (Dhawan et al. 2006).

Zoom In Zoom Out Reset image size

2.1. Spectral Analysis

The gtlike likelihood fitting tool was used to perform the spectral analysis, with "Pass 6 v3" (P6_V3) instrument response functions (IRFs); the P6_V3 IRFs are a post-launch update to address gamma-ray detection inefficiencies that are correlated with trigger rate. The 10° region around the source was modeled for Galactic and extragalactic diffuse emission, and included one nearby point source at (R.A., decl.) of (02^h23^m12^s, 62°0'0''), too faint to be found in the 3 month Bright Source List (Abdo et al. 2009a). It is important to include this nearby source in the fitting model because at low energies the PSF is sufficiently wide that despite being ∼22 away the PSF wings extend across the location of LS I +61°303 contributing approximately 13% to the flux at this position. Simultaneous modeling of this source accounts for its contribution to the flux in this region and any uncertainty is folded into the statistical error of the flux of LS I +61°303 found by the likelihood fitting tool. The 10° region was chosen to capture the broad PSF obtained at 100 MeV. An alternate fitting method using energy-dependent regions of interest was used, yielding compatible results that were folded into the systematic errors.

The Galactic diffuse emission was modeled using GALPROP, described in Strong et al. (2004) and Strong (2007), updated to include recent H i and CO surveys, more accurate decomposition into Galactocentric rings, and many other improvements, including some from comparison with LAT data (Abdo et al. 2009b). The GALPROP run designation for our model is 54_59varh7S. The diffuse sources contribute ∼95% of the observed photons shown in Figure 1.

Initially, a simple power law, E^−Γ, was fit to the orbital phase-averaged data yielding a photon index of Γ ∼ 2.42. However, as indicated in Figure 2, the energy spectrum appears to turn over at energies ∼6 GeV. The possibility of an exponential cutoff was investigated, in the form E^−Γexp [−(E/E_cutoff)]. The chance probability to incorrectly reject the power-law hypothesis was found to be 1.1 × 10⁻⁹. The best-fit exponential cutoff returns a test statistic (Mattox et al. 1996) significance value of about 4770, or roughly 70σ. The photon index is Γ = 2.21 ± 0.04 (stat) ± 0.06 (syst); the flux above 100 MeV is (0.82 ± 0.03 (stat) ± 0.07 (syst)) ×10⁻⁶ ph cm⁻² s⁻¹ and the cutoff energy is 6.3 ± 1.1(stat) ± 0.4(syst) GeV (see below for a discussion of systematics). A total of 135,659 photons were found in the 10° region. Evaluating the fit parameters, 6467 ± 80 photons were observed from LS I +61°303 above 100 MeV. Figure 2 shows the best-fit cutoff power-law model as well as the fluxes fit per energy bin and archival data from MAGIC (Albert et al. 2009) and VERITAS (Acciari et al. 2008).

Zoom In Zoom Out Reset image size

A number of effects are expected to contribute to the systematic errors. Primarily, these are uncertainties in the effective area and energy response of the LAT as well as background contamination. These are currently estimated by using outlier IRFs that bracket our nominal ones in effective area. These are defined by envelopes above and below the P6_V3 IRFs by linearly connecting differences of (10%, 5%, 20%) at log(E/MeV) of (2, 2.75, 4), respectively. Other potential sources of systematic effects investigated are: fitting technique; cuts applied (zenith angle, minimum and maximum energies); and details of the diffuse modeling. The systematic errors estimated using the bracketing IRFs were found to be greater than these additional effects, hence the bracketing IRF results were quoted for the upper limits on the systematics.

2.2. Timing Analysis

LAT light curves were extracted using aperture photometry. The LAT PSF is strongly energy dependent and, particularly since LS I +61°303 is located in the Galactic plane, there is also significant contribution to the flux within an aperture from diffuse emission and point sources that depends on the aperture size and the energy range used. The aperture and energy band employed were independently chosen to maximize the signal-to-noise level. The optimum aperture radius was found to be approximately 24 in the energy range 100 MeV–20 GeV. The time resolution of the light curve was 11,478 s, equal to twice the Fermi orbital period. Exposures were calculated using gtexposure and used to determine the count rate in each time bin. In the exposure calculation, the spectral shape is assumed to be a power law with a photon index of 2.4. The 1 day binned light curve is shown in Figure 3. Contributions from the nearby source and Galactic and extragalactic diffuse backgrounds were estimated based on the spectral fit and subtracted from the light curve.

Zoom In Zoom Out Reset image size

A search was made for periodic modulation by calculating the periodogram of the light curve (Lomb 1976; Scargle 1982). Since the exposure of the time bins was variable, the contribution of each time bin to the power spectrum was weighted based on its relative exposure. The periodogram of the unbinned, unsmoothed light curve is shown in Figure 4. The vertical line marks the Gregory (2002) orbital period and a highly significant peak is detected at this period. The significance levels marked are for a "blind" search with 500 independent frequency steps, however, the effects of the tuning of the aperture radius and energy range are not taken into account. The period and its error from the LAT observations were estimated using a Monte Carlo approach: light curves were simulated using the observed LS I +61°303 light curve and randomly shuffling the data points within their statistical errors, assuming Gaussian statistics. The corresponding periodogram was then calculated and the location of the peak at ∼26.5 days was recorded. From ∼250,000 simulations, the distribution of values gives an estimation of the measured orbital period and its associated error of 26.6 ± 0.5 days (1σ).

Zoom In Zoom Out Reset image size

The binned LAT light curve folded on the Gregory (2002) period with zero phase at MJD43, 366.2749 (Gregory et al. 1979) is shown in Figure 5. The folded light curve shows a large modulation amplitude with maximum flux occurring slightly after periastron passage. The overall light curve can be fit reasonably well by a simple sine wave, yielding a reduced χ²_ν of 1.4 for 1682 degrees of freedom (dof). However, if we use the known orbital period and ephemeris of the system (Gregory 2002) to fit a sine wave to each of the individual nine orbits observed, then we find that the best-fit amplitude varies between 6.8 ± 0.9 and 2.2 ± 0.9 ×10⁻⁷ ph cm⁻² s⁻¹, which suggests some orbit-to-orbit variability.

Zoom In Zoom Out Reset image size

2.3. Phase-resolved Spectral Analysis

The Fermi data enable for the first time the detection of a modulation in GeV gamma rays at the orbital period of a binary system. The derived period is in excellent agreement with the radio and optical-based ephemeris (Gregory 2002). The COS B source 2CG 135+01 is now firmly identified as the gamma-ray counterpart to LS I +61°303, resolving a 30 year long suspicion that the two were associated. With the identification originally based on localization only, the detection of orbital-modulated VHE emission (>100 GeV) from LS I +61°303 by MAGIC and VERITAS (Albert et al. 2006, 2009; Acciari et al. 2008) had already provided very strong support in favor of this association.

LS I +61°303 is detected at a mean flux level above 100 MeV consistent with that seen by EGRET and AGILE. Averaged over the orbital modulation, the source persists as one of the brightest high-energy gamma-ray sources in the sky over a timescale of decades (see Abdo et al. 2009a, Bright Source List, in which this source is the 15th brightest). The folded Fermi light curve peaks around phase 0.3, which is compatible with periastron passage (when the compact object is closest to the Be star) according to the latest radial velocity studies (Aragona et al. 2009). This contrasts with the behavior at VHEs where peak flux occurs at phases 0.6–0.7 and detections are achieved only at phases ranging from 0.5 to 0.8, before or at apastron. In X-rays, LS I +61°303 also appears to peak at phases 0.6–0.7 (Paredes et al. 1997; Esposito et al. 2007), whereas the radio peak occurs over a wide range of phases depending upon a four-year super-orbital cycle (Gregory 2002).

The average Fermi and EGRET spectra have compatible power-law indices and fluxes taking into account systematics, but the Fermi spectrum also shows a cutoff at approximately 6 GeV. There is no evidence for a phase dependence of the spectral shape and hence, the index or cutoff energy. VERITAS reports upper limits during the only VHE observations that are contemporary with Fermi, covering only part of one orbit from phase 0 to 0.75 (up to 2008 November 9; Holder et al. 2009). The later phases have short exposure times. Moreover, the past VHE history of the source shows several non-detections at phases 0.6–0.7 (Acciari et al. 2008; Albert et al. 2009), perhaps due to variability from one orbital cycle to the other. The Fermi light curve displays signs of orbit-to-orbit variability superposed on the mean behavior, with the primary peak always around phase 0.3. Such variability could be attributed to changing conditions in the Be star wind, affecting the interaction with the pulsar wind or relativistic jet. Indeed, optical spectra show evidence for changes in wind emission with the orbit (Zamanov et al. 1999).

The obvious radiative process to invoke in the HE and VHE range is inverse Compton scattering of the abundant stellar photons into gamma rays by a population of electrons accelerated in the vicinity of the compact object (e.g., in a relativistic jet or in a pulsar wind). Then, all else being equal, the peak flux phase is determined by where the seed photon density is highest and by geometry; favorable when the high-energy electrons are seen behind the star by the observer, e.g., Dubus et al. (2008), Khangulyan et al. (2008), and Sierpowska-Bartosik & Torres (2008). Superior conjunction is close in phase to periastron passage in LS I +61°303 (ϕ_per − ϕ_sup = 0.07 to 0.17 depending on the orbital solution). Hence, having the Fermi flux peak close to periastron is consistent with inverse Compton emission from electrons located close to the compact object. The cutoff in the average spectrum could arise due to radiative losses (because of different accelerating conditions for electrons, because of the magnetic field amplitude in the relativistic jet or the pulsar wind along the orbit and/or because of the greater photon density at periastron), or due to a varying maximum energy for accelerated electrons or to pair production on stellar photons for gamma rays above ≈50 GeV (Dubus 2006; Sidoli et al. 2006; Cerutti et al. 2008; Sierpowska-Bartosik & Torres 2009). In the latter case, cascade emission might also be seen in the Fermi range. All these effects introduce phase-dependent spectral changes. Hadronic interactions related to crossings of the Be star's equatorial wind (disk) could also contribute (Chernyakova et al. 2006). This would provide an independently varying spectral component to explain why the HE and VHE emission peak at different phases vary with orbital cycle. The expectation is that hadronic interactions would result in two asymmetric peaks in the light curve whose amplitude depends upon the intercepted matter density during the crossings and occurring at phases a priori unrelated to periastron passage but on the orientation of the orbit of the compact object relative to the Be star disk.

Continued monitoring by Fermi combined with dedicated campaigns by pointed instruments is needed to better constrain spectral variability and establish the multiwavelength connections: how do orbit-to-orbit variations compare in different energy ranges? Are there separate HE and VHE spectral components?

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase from the following agencies is also gratefully acknowledged: the Spanish CSIC and MICINN and the Istituto Nazionale di Astrofisica in Italy.

We thank the anonymous referee for useful and constructive comments.

Facility: Fermi