A Luminous and Highly Variable Gamma-Ray Flare Following the 2017 Periastron of PSR B1259–63/LS 2883

The PSR B1259−63/LS 2883 (Stephenson & Sanduleak 1971) system consists of a rotation-powered pulsar with a 48 ms spin period (first detected as a radio pulsar by Johnston et al. 1992) and a Be star in a 3.4 year (Negueruela et al. 2011), highly eccentric orbit (e = 0.87; Shannon et al. 2014). Pulsations have not been reported at any other wavelength. The Be star is surrounded by a disk-like outflow (inclined by 10°–40° from the orbit plane; Melatos et al. 1995) which the pulsar crosses through just before and just after periastron, leading to enhanced emission across the electromagnetic spectrum (e.g., Aharonian et al. 2005; Chernyakova et al. 2014). No evidence of accretion by the pulsar as it passes through the outflow has ever been reported.

Each time the pulsar crosses through the disk-like outflow of the Be star, the interaction between the pulsar wind and the outflow material leads to emission at radio, X-ray, and TeV wavelengths, producing a double-peaked light curve centered on periastron with the post-periastron hump being more luminous (see Figure 1 of Chernyakova et al. 2014, for a multiwavelength light curve). The maxima of these peaks occur approximately two weeks before and two weeks after periastron, with the X-rays peaking just before (after) the pre-(post-)periastron radio peak. In optical, the Hα and He i lines, integrated over the disk-like outflow, become increasingly activated and extended as the periastron interactions advance. The peak X-ray energy flux for the 2010 periastron passage was ∼4 × 10⁻¹¹ erg cm⁻² s⁻¹ (Chernyakova et al. 2014), while the TeV energy flux was ∼3 × 10⁻¹² erg cm⁻² s⁻¹ (Aharonian et al. 2005).

Launched in 2008 June, the Fermi Gamma-ray Space Telescope has now been operational during three periastron passages of this system. The Large Area Telescope (LAT, one of two instruments on board Fermi; Atwood et al. 2009) has now detected significant enhancements in emission above 0.1 GeV from this system following each periastron passage. Following the 2010 and 2014 periastron passages, Fermi LAT observations established that the GeV portion of the spectrum dominates the spectral energy distribution, with a peak energy flux above 0.1 GeV of ∼9 × 10⁻¹⁰ erg cm⁻² s⁻¹, and is also highly variable with weak emission, at best, before periastron and dramatic flaring after (Abdo et al. 2011; Tam et al. 2011, 2015; Chernyakova et al. 2014; Caliandro et al. 2015).

The GeV light curve during the 2014 periastron passage was broadly similar to that observed in 2010. The first detectable emission was observed 30 and 32 days post-periastron, in 2010 and 2014 respectively, rising quickly to a peak flux and then decaying smoothly, with small rises, over the next 40–50 days. In both events, the spectra of the emission were best characterized with power-law shapes falling with energy E approximately as E⁻³.

The 2010 and 2014 events present several challenges to emission models. Any viable model must be able to produce a gamma-ray luminosity (${L}_{\gamma }$) near 100% of the pulsar spin-down power ($\dot{E}=8.2\times {10}^{35}$ erg s⁻¹) during the flares (assuming a distance of d = 2.3 kpc; Negueruela et al. 2011). Recent interferometric observations of the system have resulted in an increased distance of ${2.70}_{-0.31}^{+0.41}$ kpc (Miller-Jones et al. 2018), resulting in ${L}_{\gamma }$ > $\dot{E}$ during previous periastron passages. Models must also explain the lack of a matching GeV flare occurring at approximately the same number of days before periastron. The absence of corresponding flares at other wavelengths is also a puzzle.

The GeV observations near periastron can be explained using models that invoke either inverse Compton (IC) or Doppler-boosted synchrotron emission. The model of Khangulyan et al. (2012) generates GeV gamma-rays via IC interactions between particles in the unshocked pulsar wind, leading to no contemporaneous radio or X-ray flare, and photons from the Be star and/or the disk-like outflow. The timing of the flare, ∼30 days after periastron, is a result of the pulsar wind zone increasing drastically, along the line of sight to Earth, as the pulsar exits the outflow of the Be star. Before periastron, the geometry is such that the pulsar wind zone expands in the opposite direction, explaining the lack of a pre-periastron GeV flare. Khangulyan et al. (2012) note that the photon field from the Be star is insufficient to lead to the Compotonization of the unshocked pulsar wind quickly enough to account for the near 100% efficiency of converting spin-down power into gamma-rays. They postulate that the outflow might be heated sufficiently from interactions with the pulsar, before the flare, to sufficiently increase the infrared photon density and match the observed ${L}_{\gamma }$.

Dubus & Cerutti (2013) also explained the GeV flare as the result of IC interactions, but their model uses electron–positron pairs generated near the pulsar, with a narrow energy distribution, to upscatter X-rays generated in the tail of the shocked pulsar wind. The X-rays are backscattered toward the Earth after periastron and in the opposite direction before. Their model does not rely on a Be star for target photons and matches the general shape of the GeV flare well, though it fails to produce the time delay between the flares at different wavelengths.

While it is difficult for models invoking IC emission to generate ${L}_{\gamma }$ ∼ $\dot{E}$, Doppler-boosted synchrotron emission would naturally explain this result. Kong et al. (2012) proposed that Doppler boosting of X-rays from the tail of the bow shock between the pulsar and Be star winds could explain the GeV flare. Before periastron, no boosted emission is seen and the synchrotron spectrum should peak below 100 MeV, in agreement with LAT upper limits. After periastron, a Doppler factor of ∼2 produces a spectrum that matched the LAT spectrum well. This model does not explain the lack of corresponding emission at other wavelengths during the LAT flare nor the timing of the flare ∼30 days after periastron, but a possible correlation of the X-ray spectrum with the LAT flare was reported by Tam et al. (2015). Characterizing the pattern of post-periastron GeV activity, via observing multiple orbital periods, is an important step toward understanding the emission mechanisms and dynamics of this intriguing system.

The most recent periastron passage of the PSR B1259−63/LS 2883 system occurred on MJD 58018.143 (UTC 2017 September 22 03:25:55.2) and it was quickly apparent that the behavior of the >0.1 GeV emission was different than in 2010 and 2014. He et al. (2017) reported a detection using LAT data during the month before periastron. Detections on daily timescales from 9 to 11 days post-periastron were reported by Johnson et al. (2017a). Following this early emission, no significant gamma-ray flux was detected until 40 days post-periastron. This renewed activity was more intense and displayed variability on timescales of less than six hours (Corbet 2017; Johnson et al. 2017b). Tam et al. (2018) later reported variability on 3 hr timescales. The emission lasted longer than either of the previous events observed with LAT, including a detection 70 days after periastron with a flux level of ∼1 × 10⁻⁵ cm⁻² s⁻¹ in an 8 hr interval (Johnson et al.2017c).

We report on LAT observations surrounding the 2017 periastron passage of the PSR B1259−63/LS 2883 system in more detail. We confirm the reports of rapid variability, show that the flux changes on timescales as short as ∼1.5 minutes, and report evidence for spectral curvature during the most intense emission.

We selected LAT Pass 8 data recorded between MJD 57920.143 and 58129.043 (2017 June 16 03:25:55.2 and 2018 January 11 01:01:55.2 UTC) with reconstructed directions within 15° of (R.A., decl.; J2000.0) = (195699, −63836) the radio timing position of PSR B1259−63 (Shannon et al. 2014), and energies between 0.1 and 300 GeV. This time frame provides enough pre-periastron time to obtain an acceptable fit of background sources and lasts until the flux light curve analysis (Section 2.1) gave two weeks of no detections with ≥3σ significance when fitting in one-week and one-day bins. We kept only those events with measured zenith angles ≤90°, to reduce contamination from gamma-rays associated with the limb of the Earth, and belonging to the SOURCE class, as defined under the P8R2_SOURCE_V6 instrument response functions. Finally, we filtered the data to include only time intervals when the data were marked as good and the observatory was in normal science operations mode.

We created a model of the region including all point and extended sources from the LAT four-year source catalog (3FGL; Acero et al. 2015) within 25° of PSR B1259−63 and the current Pass 8 Galactic and isotropic diffuse background models,⁸ "gll_iem_v06.fits" and "iso_P8R2_SOURCE_V6_v06.txt" (Acero et al. 2016). PSR B1259−63 does not have a 3FGL counterpart, so we added a point source at its position.

2.1. Gamma-Ray Spectrum and Light Curve

We performed a binned maximum likelihood analysis of a 20° × 20° square region centered on PSR B1259−63, using the spatial and spectral model described previously, over the entire time span of the data set. We allowed the spectral parameters of PSR B1259−63 and sources detected in 3FGL with an average significance ≥10σ and within 6° of PSR B1259−63 to vary. For those sources flagged as significantly variable in 3FGL and within 8° of PSR B1259−63, we allowed the normalization parameters to vary if they did not pass the previous criteria. The spectrum of the Galactic diffuse component was modified by a power law, normalization, and index free, while the isotropic diffuse component had only a free normalization. Because the 3FGL catalog was derived from a shorter time span and using less-sensitive Pass 7 reprocessed data, we constructed a 10° × 10° test-statistic (TS; Mattox et al. 1996) map centered on PSR B1259−63 to look for unmodeled sources. No new sources were found in this TS map, suggesting that the use of the 3FGL sources is sufficient.

Using the results of our initial fit, we constructed a new model, fixing the parameters of free point and extended sources to the best-fit values but allowing normalizations to vary. We also fixed the spectral parameters of the Galactic diffuse emission, but left the normalization of the isotropic diffuse component free. Then, we defined one-week-long intervals such that periastron, which occurred on MJD 58018.143 (2017 September 22), was at the center of a bin. We performed similar maximum likelihood fits in each one-week interval, using the new region model, with the spectrum of PSR B1259−63 modeled as a power law, ${dN}/{dE}=K{(E/{E}_{0})}^{-{\rm{\Gamma }}}$, with the normalization K and photon index Γ both free to vary while the energy scale value was fixed at E₀ = 1 GeV (Figure 1, top panel). Our single, pre-periastron detection on one-week timescales corresponds to the time frame when He et al. (2017) first saw emission rise above TS of 6 in five-day bins. Our flux light curves agree reasonably well with those of Tam et al. (2018). We do not see corresponding detections in the two bins just before periastron for our one-week bin light curve, matching their five-day bins, with our TS values being 2 and 6. This may be due to slight differences in the analysis and our choice to place periastron at the center of a bin as opposed to the start/end of a bin.

Zoom In Zoom Out Reset image size

We defined an active time when emission was seen as MJD 57923.643 to 58126.643 and tested for spectral curvature by fitting the spectrum of PSR B1259−63 as both a power law and an exponentially cutoff power law, ${dN}/{dE}=K{(E/{E}_{0})}^{-{\rm{\Gamma }}}\exp \{-(E/{E}_{{\rm{C}}})\}$, where E_C is the cutoff energy. Following Abdo et al. (2013), we calculated TS_cut, testing the significance of the curvature with respect to a power-law shape, and found a preference for the cutoff over the power law with TS_cut = 19.6 (>4σ). The best-fit spectral values for both are given in Table 1, as well as the photon (F) and energy (G) fluxes integrated from 0.1 to 300 GeV.

Table 1.  Best-fit Spectral Parameters

Parameter	Power Law	Cutoff Power Law
Active Time:
Γ	2.96 ± 0.06	2.41 ± 0.24
EC (GeV)	⋯	0.83 ± 0.40
F (10−7 cm−2 s−1)	5.99 ± 0.35	5.86 ± 0.37
G (10−10 erg cm−2 s−1)	1.95 ± 0.12	1.91 ± 0.11
Main Flare:
Γ	2.78 ± 0.06	2.15 ± 0.17
EC (GeV)	⋯	0.81 ± 0.24
F (10−7 cm−2 s−1)	11.7 ± 0.7	11.2 ± 0.8
G (10−10 erg cm−2 s−1)	4.26 ± 0.26	4.02 ± 0.24

Note. Gamma-ray spectral parameters averaged over the active time (top) and main flare (bottom) assuming a power law (column 2) and an exponentially cutoff power law (column 3) for the spectrum of PSR B1259−63. The photon F and energy G flux values are integrated over 0.1–300 GeV.

Download table as:  ASCIITypeset image

We then investigated the week with significant emission prior to periastron, the three weeks immediately after periastron, and the main flare from 40 to 75 days after periastron. During the main flare, a cutoff in the spectrum was significantly preferred over a power-law shape with TS_cut = 26.5 (∼5σ). The main flare best-fit parameters are reported in Table 1. The cutoff energy for this interval is consistent with the fit of the entire active time, and the photon index is smaller, though compatible within the uncertainties. The best-fit cutoff energy and photon index point to a need for observations in the MeV energy range. In the week before periastron, we found no evidence for curvature in the spectrum (TS_cut = 0), which is not surprising given that PSR B1259−63 is found just at our TS limit of 12 in this interval. Just after periastron, the total emission is more significant, but we did not find a strong preference for a cutoff with TS_cut = 8.8, slightly less than 3σ.

With the expectation that the spectral curvature would not be detectable on one-day timescales, we used the best-fit power-law parameters given in Table 1 to construct a one-day flux light curve. We only allowed the normalization parameters of PSR B1259−63, the isotropic diffuse emission, and sources within 3° of PSR B1259−63 to vary. The resulting flux light curve is shown in the bottom panel of Figure 1.

The pre-periastron detection, on one-week timescales, of Figure 1 is temporally consistent with the marginally significant (5 < TS < 25) pre-periastron detections, on five-day timescales, found in 2010 and 2014 (Tam et al. 2015), though our spectrum is significantly softer. In the 2010 and 2014 events, the emission started rising earlier than in 2017, reached a maximum over a few days, and then decayed, with a few lower flux peaks as the emission faded. As can be seen in the bottom panel of Figure 1, the GeV flare after the 2017 periastron consists of separate peaks, the first of which does not correspond to the highest flux level. Caliandro et al. (2015) found maximum one-day flux levels, 36 and 38 days post-periastron, of approximately 30 × 10⁻⁷ and 20 × 10⁻⁷ cm⁻² s⁻¹ for the 2010 and 2014 events, respectively. Figure 1 shows a maximum one-day flux level of (42 ± 6) × 10⁻⁷ cm⁻² s⁻¹ 56 days after periastron. A more detailed comparison of the gamma-ray and multiwavelength light curves of these three periastron passages is deferred to a future paper.

2.2. Rapid Variability

Following reports of rapid variability (Corbet 2017; Johnson et al. 2017b), we performed an aperture photometry analysis spanning nearly the entire Fermi mission in 12 hr bins (similar to the analysis in Abdo et al. 2009). We used the same data selection criteria described in Section 2, but increased the maximum zenith angle to 105° and restricted events to within 1° of PSR B1259−63. The more strict zenith angle cut described in Section 2 is necessary for a spectral analysis, but tests have shown that periodic signals are stronger with a more lenient selection and any extra contamination in such a small region of interest is acceptable (Corbet et al. 2016). No barycenter correction was applied and no background was subtracted. The resulting light curve is shown in Figure 2, and all three post-periastron gamma-ray flares are evident. The aperture photometry suggests that similar, rapid variability was not present during the 2010 and 2014 events, in agreement with the 3 hr aperture photometry of Tam et al. (2018).

Zoom In Zoom Out Reset image size

To characterize any rapid variability, we first created flux light curves on timescales of the Fermi spacecraft orbit for each one-day bin with significant emission between 40 and 75 days post-periastron. We used the spacecraft file to calculate when the zenith angle of PSR B1259−63 reached a maximum. We then defined an orbit to be from one maximum to the next. This definition has the benefit of not breaking up exposure in an orbit, depending on where in the orbital precession period the spacecraft is.

In preparation for unbinned maximum likelihood fits, we included the instrument azimuth angle dependence when calculating the exposure files, as we cannot guarantee that the position of the source can be treated as the average over this coordinate. To model the region, we used the results from fitting the active time and fixed the parameters of all sources except the normalization and photon index of PSR B1259−63 (using a power-law spectral shape), the normalization of the isotropic diffuse emission, and the normalization of 3FGL J1305.7−6241, the brightest point source within 3° of PSR B1259−63. We only analyzed data for orbits in which there were at least 10 events recorded within 15° of PSR B1259−63. This last requirement was made to avoid erroneous results giving unrealistically high flux levels, based on experience analyzing the data in these short intervals.

The combined orbit-by-orbit flux light curve and photon index variation are shown in Figure 3. Each point is centered on orbit noon, the point in each orbit when PSR B1259−63 reaches a minimum zenith angle, and the horizontal extent reflects the time period when PSR B1259−63 is within 60° of the LAT boresight, not the start and stop of the orbit, as this typically gives a better idea of when events could be recorded from the source.

Zoom In Zoom Out Reset image size

Figure 3 displays significant variation on timescales as short as one orbit, with flux levels up to ∼7 times what a one-day analysis would suggest. The photon index stays relatively stable across the orbits we analyzed.

We selected the five orbits with the highest TS values to probe variability on even shorter timescales. For each of these orbits, we created a 30 s time bin, starting at the time of the first event, and required at least 10 events in the 15° radius region or else we enlarged the time bin in 30 s steps until this criterion was met. Next, we performed an unbinned maximum likelihood fit, with only the normalizations of PSR B1259−63 and the isotropic diffuse component free to vary, and required that PSR B1259−63 be detected with TS ≥ 9 and at least 4 predicted counts, otherwise we enlarged the bin by 30 s and repeated the fit. Once those criteria were met, we recorded the best-fit information and started the process anew. If the last time bin did not meet the TS or predicted number of counts criteria, we compared the flux of the previous bin to that derived by merging that bin with the last bin. If these two flux values were comparable, we kept the merged bin, otherwise we report a 95% confidence-level upper limit for the last bin. The resulting intra-orbit light curves are shown as numbered insets in Figure 3.

For these orbits, the bins range in size from 1 to ∼20 minutes, with flux values doubling or tripling over timescales as small as 1–1.5 minutes. As a measure of the significance of the intra-orbit variability, we calculated $\delta \mathrm{TS}={\mathrm{TS}}_{{\rm{\Sigma }}}-{\mathrm{TS}}_{\mathrm{orb}}$ where TS_orb is the TS value found when fitting the entire orbit and TS_Σ is the sum of the TS values in each intra-orbit bin. This gave ΔTS = 11, 84, 56, 26, and 61 for insets 1 through 5, in order. Testing against the hypothesis that the flux does not vary over an orbit, the ΔTS values should follow a ${\chi }_{\mathrm{nbin}-1}^{2}$ distribution. Under this assumption we find significances of 1.1σ, 7.6σ, 5.8σ, 1.5σ, and 6.1σ, respectively. These results suggest that the intra-orbit light curves in insets 1 and 4 show only marginal variability, but this may be due to our method emphasizing the highest time resolution over the most significant binning.

For inset 3, following the prescription described previously resulted in the first bin having an extremely high flux. However, upon investigation, we discovered that this bin had unusually low exposure and all of the events had zenith angle values >80°. We therefore merged that bin with the one after it, resulting in a more reasonable flux value. We verified that a similar issue did not occur in the other orbits we analyzed in detail.

As discussed in Section 1, proposed emission models invoke either IC or Doppler-boosted synchrotron emission to explain the post-periastron GeV flares from the PSR B1259−63 system. IC models can be hard-pressed to generate emission with ${L}_{\gamma }$ ∼ $\dot{E}$, while synchrotron models predict contemporaneous flares at other wavelengths. The results presented in previous sections provide even more challenges for emission models.

A useful value when discriminating between different emission mechanisms is ${L}_{\gamma }$/$\dot{E}$. In Table 2, we consider the maximum luminosity on the different timescales we have analyzed, defining ${L}_{\gamma }$ = $4\pi {{Gd}}^{2}$ and assuming a distance of $d={2.70}_{-0.31}^{+0.41}$ kpc. On daily timescales, the observed gamma-ray luminosity pushes the limits of what the spin-down power, $\dot{E}$, can provide. The variability on even shorter timescales requires significant beaming or an additional source of energy.

Table 2.  Maximum Gamma-Ray Energetics on Different Timescales

Timescale	G	${L}_{\gamma }$	${L}_{\gamma }$/$\dot{E}$
	(10−10 erg cm−2 s−1)	(1035 erg s−1)
One-week	7.3 ± 0.6	${6.4}_{-1.6}^{+2.0}$	0.8 ± 0.2
One-day	14 ± 2	${12}_{-3}^{+4}$	${1.5}_{-0.4}^{+0.5}$
One-orbit	70 ± 16	${61}_{-14}^{+18}$	${7.4}_{-1.7}^{+2.2}$
Intra-orbit	280 ± 100	${244}_{-56}^{+74}$	${29.8}_{-6.8}^{+9.0}$

Note. For the timescales listed during the 2017 periastron passage, this table provides the maximum energy flux (G), gamma-ray luminosity (${L}_{\gamma }$), and luminosity as a fraction of the spin-down power, $\dot{E}=8.2\times {10}^{35}$ erg s⁻¹ (${L}_{\gamma }$/$\dot{E}$). For the uncertainty on ${L}_{\gamma }$, we incorporate both the energy flux and distance uncertainties.

Download table as:  ASCIITypeset image

Our results for the 2017 periastron passage disfavor models that generate GeV emission primarily through the IC mechanism, which cannot easily produce ${L}_{\gamma }$ > $\dot{E}$ (e.g., Khangulyan et al. 2012). However, given the stark differences in the progression of the 2017 event from the 2010 and 2014 events, it is possible that a different mechanism was responsible for the emission in 2017, though it is currently unclear why this would be the case.

Following Tam et al. (2011), if we assume that the emission near periastron represents the unboosted synchrotron flux with daily energy flux values of ∼4 × 10⁻¹⁰ erg cm⁻² s⁻¹ and compare to our maximum intra-orbit energy flux, using Equation (2) of Dubus et al. (2010) the ratio of these fluxes should be ${D}^{3+\alpha }$, where D is the Doppler factor and α = Γ − 1 = 1.17 using the best-fit photon index from that orbit. This yields D ∼ 3 as compared to values between 1.5 and 2 found by Tam et al. (2011) for the 2010 event.

Kong et al. (2012) presented a model for the post-periastron GeV flares as Doppler-boosted synchrotron emission and estimated this would result in a cutoff energy of ∼0.236ζD GeV, where ζ is the electron acceleration efficiency. They chose a value of ζ = 0.36 so that the synchrotron emission would not contribute significantly above 0.1 GeV before periastron. If we use the best-fit cutoff energy from the main flare data and D = 3, we find ζ ∼ 1.1, implying a high acceleration efficiency.

Arguably the most exciting result from the 2017 periastron passage is the rapid variability. If we take the shortest rise/fall time of the intra-orbit variability (insets in Figure 3) to be the implied light crossing time, Δt ∼ 1.5 minutes, then following Ackermann et al. (2016) we can place an upper limit on the radius of the emission region of R ≲ DcΔt = 8 × 10⁷ km, where c is the speed of light, and using D = 3. Using the orbital inclination angle of 1534 derived by Miller-Jones et al. (2018), we calculate the distance between the pulsar and the system center of mass, close to the center of the Be star, to be 1.8 × 10⁸ km and 2.6 × 10⁸ km 40 and 72 days after periastron, respectively, when we observed intra-orbit variability. Our upper limit on the radius of the emission region thus varies from ∼40% to ∼30% of the distance between the pulsar and the Be star during the main flare. New models must therefore explain the flux variation on minute timescales and adhere to this limit on the emission region.

We have shown that the GeV gamma-ray emission associated with the 2017 periastron passage of the PSR B1259−63 system progressed in a wholly different manner than the 2010 and 2014 events. Significant emission was seen for ∼3 weeks immediately following periastron and then not again until 40 days after periastron when several flares, lasting a few days each and reaching daily integral photon fluxes as high as 4 × 10⁻⁶ cm⁻² s⁻¹, were observed. The flares demonstrated significant spectral curvature and variability on timescales shorter than the Fermi spacecraft orbital period with integral photon flux levels, on minute timescales, as high as ∼4 × 10⁻⁵ cm⁻² s⁻¹. Excluding gamma-ray bursts and solar flares, this is the fastest GeV variability observed in LAT data. Our results challenge the existing emission models of the PSR B1259−63/LS 2883 system, during the periastron passage and the few months after, and demonstrate the need for continued monitoring over additional orbital cycles.

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 the science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

Portions of this research performed at the Naval Research Laboratory are sponsored by NASA DPR S-15633-Y and Fermi Guest Investigator grant 16-Fermi10-0006.

Facility: Fermi(LAT). -