FERMI/LAT OBSERVATIONS OF LS 5039

LS 5039, PSR B1259−63 and LS I +61°303 are the only binaries with high-mass companions long known to be spatially coincident with sources detected at energies greater than 100 MeV, e.g., those listed in the Third Energetic Gamma-Ray Experiment (EGRET) catalog (Hartman et al. 1999). The latter binary was detected by the Large Area Telescope (LAT) on the Fermi mission, confirming it as a GeV gamma-ray source and finding variability that is consistent with modulation on the binary period of 26.6 ± 0.5 days (Abdo et al. 2009b). This constituted the first detection of orbital periodicity in high-energy (HE) gamma rays (20 MeV–100 GeV). In this Letter, we present the results of Fermi observations of LS 5039.

LS 5039 is one of a handful of X-ray binaries that have been detected recently at very high-energy γ-rays; results from ∼70 hr of observations distributed over many orbital cycles have been presented by Aharonian et al. (2005b, 2006). These observations yielded a modulation of the very high energy (VHE, >100 GeV) gamma-ray flux with a period of 3.9078 ± 0.0015 days (Aharonian et al. 2006), consistent with the orbital period as determined by Casares et al. (2005) from optical spectroscopy. Short timescale variability displayed on top of this periodic behavior, both in the flux and the spectrum, was also reported.

The nature of the LS 5039 compact object is unknown: a black hole and a neutron star have been invoked as possible compact object companions in a slightly eccentric (e ∼ 0.35), 3.90603 ± 0.00017 day orbit around the O6.5V star (Casares et al. 2005). The discovery of a jet-like radio structure in LS 5039 coincident with an X-ray and EGRET source prompted a microquasar interpretation (Paredes et al. 2000a). Variability in the EGRET source could not be established precisely (Torres et al. 2001; Nolan et al. 2003). Recently, it has been shown that the X-ray flux is modulated on the orbital period (Takahashi et al. 2009).

Ribó et al. (2008) provided Very Long Baseline Array (VLBA) radio observations of LS 5039 with morphological and astrometric information at milliarcsecond scales that cannot easily be explained by a microquasar scenario. Martocchia et al. (2005) assessed the low X-ray state, showing the absence of accretion features in the X-ray spectra. Thus, measurements at radio and VHE γ-rays in the cases of LS I +61°303 (Dhawan et al. 2006; Albert et al. 2008) or PSR B1259−63 (Aharonian et al. 2005a), whose overall spectral energy distribution is similar to that of LS 5039, gave the impression that all three systems are different realizations of the same scenario: a pulsar-massive star binary (Dubus 2006b).

Theoretical computations of the gamma-ray emission in both compact object scenarios have been made, with gamma rays produced by inverse Compton (IC) scattering of the stellar light by VHE electrons accelerated in the vicinity of the compact object. In the case of the black hole companion, HE and VHE emission would result from particles accelerated in the jet (Bednarek 2007a; Böttcher 2007; Khangulyan et al. 2008). Alternatively, it would involve the relativistic wind of a young, rotation-powered pulsar, either as a result of particle acceleration in the wind interaction region (Dubus 2006b) or by processes within the pulsar wind (Sierpowska-Bartosik & Torres 2008a, 2008b; Cerutti et al. 2008).

The LAT onboard Fermi 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).

The analysis data set spanned 2008 August 4 through 2009 June 22. The data were reduced and analyzed using the Fermi Science Tools v9r15 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 ("Pass 6 diffuse") event class is used. Throughout the analysis, the "Pass 6 v3 Diffuse" (P6_V3_DIFFUSE) instrument response functions (IRFs) are applied.

Time periods when the region around LS 5039 was observed at a zenith angle greater than 105° and for observatory rocking angles of greater than 43° were also excluded to avoid contamination from the Earth albedo photons. With these cuts, a photon count map of a 10° region around the binary is shown in Figure 1.

Zoom In Zoom Out Reset image size

LS 5039 is detected at a level of 28.5σ (see Section 4). The gtfindsrc tool finds a best-fit location for LS 5039 of R.A. = 18^h26^m247, decl. = −14°48'39''(J2000) with a 95% error of 0.054° (including a 20% systematic error derived from the internal Fermi 11-month catalog). The nominal position of LS 5039 at R.A. = 18^h26^m151, decl. = −14°50'543 (J2000) (Zacharias et al. 2004) lies just on the Fermi 95% contour; the nominal position was used throughout the analysis.

LAT light curves were extracted using aperture photometry. The LAT point-spread function (PSF) is strongly energy dependent and, particularly since LS 5039 is located in the Galactic plane, there is also significant contribution to 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 0925 in the energy range 100 MeV–10 GeV. The time binning of the light curve was 1000 s. Exposures were calculated using gtexposure assuming a power-law spectrum with a photon index of Γ = 2.5.

A search was made for periodic modulation by calculating the weighted periodogram of the light curve (Lomb 1976; Scargle 1982; Corbet & Dubois 2007). 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 is shown in Figure 2. The arrow marks the Casares et al. (2005) orbital period and a highly significant peak is detected at this period; the false alarm probability is ∼10⁻¹⁰. The significance levels marked are for a "blind" search with 5000 independent frequency steps, however, the effects of the tuning of the aperture radius and energy range are not taken into account. The period error was estimated using a Monte Carlo approach: light curves were simulated using the observed LS 5039 light curve and randomly shuffling the data points within their errors, assuming Gaussian statistics. The corresponding periodogram was then calculated and the location of the peak at ∼3.9 days recorded. From ∼200,000 simulations, we calculate an error estimate of the orbital period of 3.903 ± 0.005 days (1σ).

Zoom In Zoom Out Reset image size

The gtlike likelihood fitting tool was used to perform the spectral analysis, wherein a spectral-spatial model containing point and diffuse sources is created and the parameters obtained from a simultaneous maximum likelihood fit to the data. The 10° region around the source was modeled for Galactic and isotropic diffuse contributions and 22 additional significant point sources (taken from the internal Fermi 11-month catalog) were included; point sources were modeled with simple power laws with the exception of the bright, nearby pulsar (PSR J1826−1256) for which a power law with an exponential cutoff was used. The flux contribution of PSR J1826−1256 to the region was minimized by excluding events which arrive during the peaks in the pulsar phase cycle (0.175 <ϕ< 0.3 and 0.625 <ϕ< 0.775 were excluded; see Abdo et al. 2009a). A scaling factor of 1/0.725 is applied to measured fluxes to account for livetime loss due to this cut.

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.

In analyzing the emission of LS 5039, we used models for the Galactic diffuse emission (gll\_iem\_v02.fit) and isotropic backgrounds currently recommended by the LAT team.⁵⁹ The model for the Galactic diffuse emission was developed using spectral line surveys of H i, CO (as a tracer of H₂) to derive the distribution of interstellar gas in Galactocentric rings.

The model of the diffuse gamma-ray emission was then constructed by fitting the gamma-ray emissivities of the rings in several energy bands to the LAT observations. The fitting also required a model of the IC emission calculated using GALPROP (Strong 2007) and a model for the isotropic diffuse emission. The latter was fitted to the LAT data using an analysis of the sky above |b| = 30° and includes the significant contribution of residual (misclassified) cosmic rays at high energies for the current IRFs.

Initially a simple power law, E^−Γ, was fitted to the data from all orbital phases yielding a photon index of Γ ∼ 2.54. However, as can be clearly seen in Figure 3, the energy spectrum appears to turn over at energies above ∼2 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.6 ×10⁻¹⁶. The maximum likelihood exponential cutoff spectral model has a likelihood test statistic (Mattox et al. 1996) value of ∼814.6, or 28.5σ. The photon index is Γ = 1.9 ± 0.1(stat) ± 0.3(syst); the 100 MeV–300 GeV flux is (4.9 ± 0.5(stat) ± 1.8(syst)) × 10⁻⁷ photon cm⁻² s⁻¹ and the cutoff energy is 2.1 ± 0.3(stat) ± 1.1(syst)GeV (see below for a discussion of systematics). The correlation between the photon index and cutoff energy was explored by fitting a family of models over a grid of indices and cutoff energies centered on the best-fit parameter values. The 1σ, 2σ, and 3σ error contours are shown in Figure 4; both parameters are bounded and well constrained. A total of 359,789 photons were found in the 10° region. Evaluating the fit parameters, 3992 ± 63 photons were observed from LS 5039 above 100 MeV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A number of effects are expected to contribute to the systematic errors. The largest is uncertainty in the diffuse modeling as evidenced by the intense swath of photons along the Galactic plane shown in Figure 1. A reasonable range of shape difference has been explored using the GALPROP model of the region. Both models give reasonable residuals maps and show differences of 10%, 37%, and 50%, respectively, for index, cutoff energy and flux. In all diffuse models tested, the exponential cutoff model is a significant improvement over the power law.

The impact of systematic uncertainties due to the IRFs are 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%, and 20%) at log(E/MeV) of (2, 2.75, and 4), respectively. The variation on the index was 15%; the other parameter variations were small compared to those due to the diffuse modeling. Other effects considered are fitting technique, cuts applied (minimum and maximum energies), but they are all within the ranges defined by diffuse modeling and bracketing effective area variations.

4.1. Phase-Resolved Spectral Analysis

We also searched for orbital dependence of the spectral shape. gtlike fits were performed for each phase interval of Δϕ = 0.1 in the same way as for the phase-averaged spectral analysis. Figure 5 shows the hardness ratio, F_{1−100 GeV}/F_{0.1−1 GeV}, as a function of orbital phase. Due to limited statistics, bin widths of Δϕ of 0.2 are used for the phase interval ϕ = 0.3–0.6, while Δϕ = 0.1 are employed for the other phase intervals. The spectral shape varied such that the spectrum is softer around periastron and is harder around apastron.

Zoom In Zoom Out Reset image size

Aharonian et al. (2006) define two broad phase intervals, inferior conjunction (0.45 <ϕ< 0.9) and superior conjunction (0.9 <ϕ<0.45), both H.E.S.S. spectra being shown in Figure 3. Taking the same phase intervals with the LAT data, we find a power-law spectrum with Γ = 2.25 ± 0.11 at inferior conjunction; an energy cutoff was not statistically justified. At superior conjunction, a power law with an exponential cutoff was preferred with Γ = 1.91 ± 0.16 and a cutoff energy of 1.9 ± 0.5 GeV. Aharonian et al. (2006) report spectral variability in the source between superior and inferior conjunction, however, they do not see any indication of long-term variability suggesting that it may be reasonable to compare non-contemporaneous observations.

The initial association of LS 5039 with the EGRET source 3EG J1824−1514 (Paredes et al. 2000b) had remained tentative due to the large EGRET error circle and the lack of timing signatures. The association was bolstered by the discovery of point-like, modulated gamma-ray emission above 100 GeV (Aharonian et al. 2005b, 2006). The Fermi observations enabled the detection of an orbital modulation, indicating that the binary is also a source of gamma rays above 100 MeV. The periods determined independently from the Fermi-LAT and H.E.S.S. data are self consistent and compatible with the binary period obtained from radial velocity measurements of the companion (Aragona et al. 2009). This is only the second high-mass X-ray binary after LS I +61°303 (Abdo et al. 2009b) to become a confirmed emitter in the HE gamma-ray domain.

The short orbital period means that the compact object passes within a stellar radius of the surface of the kT_⋆ ≈ 3 eV, R_⋆ ≈ 10 R_☉ companion (Casares et al. 2005). Gamma rays emitted in the vicinity of the compact object with energies above the 30 GeV threshold inevitably pair produce with stellar photons (see, e.g., Protheroe & Stanev 1987; Moskalenko 1995; Böttcher & Dermer 2005; Dubus 2006a). Emission in the Fermi range is largely unaffected by absorption and should allow easier identification of the intrinsic spectrum and variability of the source. However, it can be affected by cascading of higher-energy photons.

The HE modulation peaks at ϕ∼ 0.0–0.1, close to superior conjunction (ϕ = 0.06) while the VHE modulation peaks close to inferior conjunction (ϕ = 0.72). The phase difference can be explained mostly as a result of the competition between IC scattering on high-energy electrons and VHE pair production, assuming the companion star to be close to the compact object. The star provides target photons for both processes with the radiation density varying by a factor 4 along the eccentric orbit. IC scattering will vary with radiation density but, because the source of seed photons is anisotropic, the flux will also depend on the geometry seen by the observer in non-trivial ways (Khangulyan et al. 2008; Sierpowska-Bartosik & Torres 2008b). The emission is enhanced (reduced) when the highly relativistic electrons seen by the observer encounter the seed photons head-on (rear-on), i.e., at superior (inferior) conjunction. Inversely, VHE absorption due to pair production will be maximum (minimum) at superior (inferior) conjunction. The phases of minimum and maximum flux in Fermi as well as the anti-correlation with H.E.S.S. are consistent with these expectations, suggesting IC scattering is the dominant radiative process above 100 MeV with the additional effect of pair production above 30 GeV (Bednarek 2007b; Dubus et al. 2008; Sierpowska-Bartosik & Torres 2008a). It is, however, yet unclear whether the IC VHE emission is mainly produced already in the pulsar wind zone of the system, given the high opacity already found therein for particles accelerated at the pulsar site (Sierpowska-Bartosik & Torres 2008b) or as a result of particle acceleration at the shock formed in the wind collision region (Dubus 2006a) or even well beyond the system (Bosch-Ramon et al. 2008).

However, the extension of these pictures from the TeV into the GeV domain is undermined by the presence of an exponential cutoff at a few GeV in the Fermi spectrum. A cutoff at ∼6.3 GeV was also observed in LS I +61°303  indicating that this may be a common spectral feature in this class of source. The companion star in LS I +61°303 is a Be star with a dense equatorial disk. Passage of the compact object through this disk might have explained the exponential cutoff: for instance, this would crush a putative pulsar wind nebula closer to the neutron star, increasing synchrotron losses and introducing a strong dependence with orbital phase of the electron energy distribution (Dubus 2006b). But there is no such large, systematic contrast in the density of the stellar wind from the O6.5V star along the orbit in LS 5039. The presence of a similar cutoff in both systems argues against explanations related to the properties of the orbit or companion star. The cutoff seems to require that the radiative process or radiating electrons be different between the HE and VHE domains, in disagreement with the picture presented above.

An intriguing possibility is that the emission in the Fermi range from both LS 5039 and LS I +61°303 is magnetospheric emission as seen in the dozens of pulsars that have now been detected by Fermi. The typical Fermi pulsar emission has a hard power-law spectrum with a photon index ≈1.5 and an exponential cutoff at ≈2.5 GeV. In this case, the cutoff energy is thought to be set by the balance between acceleration and losses to curvature radiation. The emission should in such a case be pulsed, although the orbital motion makes it exceedingly difficult to find in the Fermi data with no prior knowledge of the spin period. Pulsar wind emission would dominate in the neighboring hard X-ray and VHE bands, supported by their similar orbital modulations (Hoffmann et al. 2009; Takahashi et al. 2009). It is to be noted that there is, however, an issue associated with having dominant magnetospheric emission in the HE band: magnetospheric emission is usually thought to be due to curvature radiation and has no obvious reason to be modulated with the orbital motion (although magnetospheric emission from pulsars in close binaries has hardly been modeled). The dense photon environment of the binaries may perhaps introduce differences compared to gamma-ray magnetospheric emission from isolated pulsars (for instance, in the pair multiplicity). Hence, the spectrum suggests a magnetospheric emission interpretation, which is hard to reconcile with the IC scattering interpretation suggested by the modulation. Future HE and VHE should lead to better constraints on the variability of the spectral parameters along the orbit. This will help resolve whether there are several components and what their relative amplitudes are.

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 is gratefully acknowledged from the Spanish CSIC and MICINN, the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

Facility: Fermi