GeV Detection of HESS J0632+057

Gamma-ray binaries are binary systems producing most of their electromagnetic output in gamma-rays above 1 MeV (for a review, see Dubus 2015). They show orbitally modulated emission at essentially all frequencies. There are only a handful of gamma-ray binaries known: five in the Galaxy (PSR B1259-63, Aharonian et al. 2005a; Abdo et al. 2011; Caliandro et al. 2015; LS I +61 303, Albert et al. 2006; Abdo et al. 2009c; Hadasch et al. 2012; LS 5039, Aharonian et al. 2005b, 2006; Abdo et al. 2009a; Hadasch et al. 2012; Collmar & Zhang 2014; 1FGL J1018.6-5856, Li et al. 2011b; Ackermann et al. 2012b; Abramowski et al. 2015; HESS J0632+057, Aharonian et al. 2007; Bongiorno et al. 2011; Aliu et al. 2014) and one in the Large Magellanic Cloud (CXOU J053600.0-673507, Corbet et al. 2016). Cyg X-1 (McConnell et al. 2000; Albert et al. 2007; Sabatini et al. 2010) and Cyg X-3 (Abdo et al. 2009d; Tavani et al. 2009; Corbel et al. 2012) have also been detected in gamma-rays. However, their spectral energy distributions (SEDs) peak at X-ray energies, and their gamma-ray emission is not recurrent in every orbit. The currently known gamma-ray binaries are all high mass X-ray binary systems, hosting a massive O or Be star and a compact object. Except for PSR B1259-63, hosting a 48 ms pulsar, the nature of the compact objects in such binaries is unknown. Pulsar/stellar wind interaction (e.g., Maraschi & Treves 1981; Dubus 2006), pulsar wind zone processes (e.g., Sierpowska-Bartosik & Torres 2008; Bednarek 2011; Bednarek & Sitarek 2013), a transitioning pulsar scenario (e.g., Zamanov et al. 2001; Papitto et al. 2012; Torres et al. 2012), and microquasar jets (see, e.g., Bosch-Ramon & Khangulyan 2009 for a review) have been proposed as the origin of the gamma-ray emission for one or several gamma-ray binaries.

HESS J0632+057 was discovered as an unidentified TeV point source close to the rim of the Monoceros supernova remnant (SNR) and was proposed to be associated with the B0Vpe star MWC 148 (Aharonian et al. 2007). Follow-up XMM-Newton observations of HESS J0632+057 revealed a bright X-ray source, XMMU J063259.3+054801, positionally coincident with HESS J0632+057 and MWC 148 (Hinton et al. 2009). The low probability of a random coincidence between sources like HESS J0632+057 and MWC 148 (∼10⁻⁴, Aharonian et al. 2007), or between sources like MWC 148 and XMMU J063259.3+054801 (∼10⁻⁶, Hinton et al. 2009) strengthens the argument for a physical association. Since an isolated star is unlikely to accelerate particles to very high energy (VHE) (≫1 TeV), Hinton et al. (2009) proposed MWC 148 to be part of a binary system, concurrently classifying HESS J0632+057 as a new gamma-ray binary. Subsequent VERITAS observations of HESS J0632+057 did not yield any detection above 1 TeV (Acciari et al. 2009), implying a significant flux variability. Swift/XRT observations confirmed this flux variability in X-rays, from which a lower limit to the orbital period was estimated as ≥54 days (Acciari et al. 2009; Falcone et al. 2010). A similar constraint (>100 days) was obtained by Aragona et al. (2010) via optical spectroscopy.

The mass and radius of MWC 148 were estimated to be in the range of 13.2–19.0 M_⊙ and 6.0–9.6 R_⊙, respectively (Aragona et al. 2010). By fitting the SED, Aragona et al. (2010) proposed it to be at a distance between 1.1 and 1.7 kpc. The radio counterpart of HESS J0632+057 was detected both at 1280 MHz with the Giant Metrewave Radio Telescope and at 5 GHz with the Very Large Array (VLA; Skilton et al. 2009). The radio properties are consistent with established gamma-ray binary systems. With additional, years long Swift/XRT monitoring of HESS J0632+057, an orbital period of 321 ± 5 days was revealed, establishing the binary nature of HESS J0632+057 (Bongiorno et al. 2011). The orbital period was further refined to be ${315}_{-4}^{+6}\,\mathrm{days}$ by Aliu et al. (2014), and the eccentricity of the binary orbit was estimated as 0.83 ± 0.08 with a mass of the compact object in the range of 1.3–7.1 M_⊙ (Casares et al. 2012).

By analyzing Chandra and XMM-Newton observations, significant flux and spectral variability between the high and low X-ray states of HESS J0632+057 were reported by Rea & Torres (2011). No pulsed emission from HESS J0632+057 was found, leading to a 3σ upper limit on the X-ray pulsed fraction of ∼30% (Rea & Torres 2011). This should be compared to the upper limits on the pulsed fraction of LS 5039 (15%), or of LS I +61 303 (10%; Rea et al. 2010, 2011).

HESS J0632+057 showed aligned orbital light curves in X-ray and TeV with an apparent peak in the orbital phase range 0.2–0.4 (Aleksić et al. 2012; Aliu et al. 2014). The X-ray peak of HESS J0632+057 is ∼0.3 of the orbit after periastron, similar to the case of LSI +61 303 (Torres et al. 2010; Zhang et al. 2010; Li et al. 2011a). HESS J0632+057 was detected as an extended radio source with a projected size of ∼75 au by the European VLBI Network at 1.6 GHz (Moldón et al. 2011). Its morphology, size, and displacement on au scales are similar to those found in other gamma-ray binaries.

All Galactic gamma-ray binaries have been detected in the high energy (HE; >100 MeV) and VHE (>100 GeV) range except for HESS J0632+057, which is a bright TeV source detected down to 136 GeV (Aliu et al. 2014 and references therein), but remained undetected in the GeV range (Caliandro et al. 2013). The latter authors carried out a search for HESS J0632+057 in the 0.1–100 GeV range using 3.5 years of Fermi-LAT data, which led to a 95% CL flux upper limit of 3 × 10⁻¹¹ erg cm⁻² s⁻¹. Recently, Malyshev & Chernyakova (2016) reported the detection of HESS J0632+057 at ∼5σ significance in the highest energy band of Fermi-LAT (200–600 GeV), at orbital phase 0.2–0.4 and 0.6–0.8. We discuss these results in detail below. In this paper, we report on a detailed search for gamma-ray emission from HESS J0632+057 in the GeV energy range, using nearly nine years of Fermi-LAT data.

The Fermi-LAT data included in this paper cover the period from 2008 August 4 to 2017 April 2. The analysis of the Fermi-LAT data was performed using the Fermi Science Tools,⁹ 11-05-02 release. Photons from the "P8 Source" event class (evclass = 128) and "FRONT+BACK" event type (evtype = 3) were selected.¹⁰ The "Pass 8 R2 V6" instrument response functions (IRFs) were used in the analysis. All photons in the energy range of 0.1–300 GeV and within a circular region of interest (ROI) of 10° radius centered on HESS J0632+057 were considered. A larger ROI of 15° radius leads to consistent results. To reject contaminating gamma-rays from the Earth's limb, only events with zenith angle <90° were selected.

The gamma-ray flux and spectral results presented in this work were calculated by performing a binned maximum likelihood (Mattox et al. 1996) fit using the tool gtlike. The spectral–spatial model constructed to perform the likelihood analysis includes Galactic and isotropic diffuse emission components ("gll_iem_v06.fits", Acero et al. 2016, and "iso_P8R2_SOURCE_V6_v06.txt", respectively¹¹ ) as well as known gamma-ray sources within 15° of HESS J0632+057, based on a preliminary seven-year source list. The spectral parameters of these sources were fixed at the source list values, except for those within 3° of our target, for which all the spectral parameters were left free. Due to the presence of the bright gamma-ray pulsar PSR J0633+0632 in the vicinity of HESS J0632+057, photons within a specific pulsar spin phase interval are selected, as explained in more detail in Section 3. The Test Statistic (TS) was employed to evaluate the significance of the gamma-ray fluxes coming from the sources. It is defined as $\mathrm{TS}=-2\mathrm{ln}({L}_{\max ,0}/{L}_{\max ,1})$, where L_max,0 is the maximum likelihood value for a model in which the source studied is removed (the "null hypothesis") and L_max,1 is the corresponding maximum likelihood value for the full model. The larger the value of TS, the less likely the null hypothesis is correct (i.e., a significant gamma-ray excess lies on the tested position) and the square root of the TS is approximately equal to the detection significance of a given source. A TS value greater than 25 was required for the inclusion in the preliminary seven-year source list. TS maps in this paper are produced with the pointlike analysis package (Kerr 2011). The systematic errors have been estimated by repeating the analysis using modified IRFs that bracket the effective area¹² (Ackermann et al. 2012a) and artificially changing the normalization of the Galactic diffuse model by ±6% (Abdo et al. 2013). The first (second) uncertainty shown in the paper corresponds to statistical (systematic) error.

HESS J0632+057 is located in a complicated region. Within 3° of its location, there are several gamma-ray point sources from the preliminary seven-year source list, the Monoceros Loop SNR, and the Rosette Nebula, which are both known to be extended gamma-ray sources (Katagiri et al. 2016). Located ∼1° away from the source of interest, there is PSR J0633+0632, a bright, radio-quiet gamma-ray pulsar discovered in the first six months of Fermi-LAT observations (Abdo et al. 2009b). To minimize contamination from this pulsar, we gate off the pulsed emission from PSR J0633+0632, following a method similar to that used in The Second Fermi Large Area Telescope Catalog of Gamma-Ray Pulsars (Abdo et al. 2013, 2PC hereafter). We selected photons from PSR J0633+0632 within a radius of 06 and a minimum energy of 500 MeV, which maximized the H-test statistics (de Jager et al. 1989; de Jager & Büsching 2010). The current timing ephemeris for PSR J0633+0632¹³ has been extended to cover the Fermi-LAT data considered in this paper using the method described by Ray et al. (2011). Adopting the updated ephemeris, we assigned pulsar rotational phases to each gamma-ray photon that passed the selection criteria, using Tempo2 (Hobbs et al. 2006) with the Fermi plug-in (Ray et al. 2011). The pulse profile of PSR J0633+0632 is shown in Figure 1. To define off-peak intervals, we have deconstructed the pulsed light curve into simple Bayesian Blocks using the same algorithm described in the 2PC, details of which can be found in Jackson et al. (2005) and Scargle et al. (2013). The off-peak phases are defined as ϕ = 0.213–0.510 and 0.620–1.0 and are shown in Figure 1.

Zoom In Zoom Out Reset image size

The analysis of the region surrounding HESS J0632+057 was performed using the data in the off-peak phases of PSR J0633+0632 (Figure 1). To account for the off-peak phase selection, the prefactor parameter of all sources were scaled by 0.677. To account for the gamma-ray emission of the Rosette Nebula, beyond the LAT standard diffuse model, we adopted a spatial template based on the CO line emission, similar to that used in Katagiri et al. (2016; the spatial model is shown as green contours in Figure 2). A LogParabola spectral model was adopted, also following Katagiri et al. (2016). Point sources from the preliminary seven-year source list located within the spatial template were not included. Similarly, to account for the gamma-ray emission of the Monoceros Loop, we adopted the Gaussian emission profile and a LogParabola spectral model reported in Katagiri et al. (2016). PSR J0632+0646 and PSR J0633+0632 are included in the spatial model following Katagiri et al. (2016). Other point sources from a preliminary source list based on seven years of LAT data located within the central region of the Gaussian profile were not included. An additional point source modeled by a simple power law was added to the spatial model in the Monoceros Loop region (Figure 2). The best position of the additional source was determined with pointlike as R.A. = 9929 ± 0.07, decl. = 621 ± 0.06. The likelihood analysis of the new point source yields a TS = 33, a photon index of 2.45 ± 0.05 and an energy flux of (0.71 ± 0.15) × 10⁻¹¹ erg cm⁻² s⁻¹ in the 0.1–300 GeV range. We also tested alternative spatial modeling of the Monoceros Loop region. For instance, we used point sources from the preliminary seven-year source list plus a collection of a few additional point sources (following the method described by Caliandro et al. 2013), which yields consistent results on HESS J0632+057.

Zoom In Zoom Out Reset image size

Figure 2 shows the TS map calculated with the Rosette Nebula, the Monoceros Loop, and the new point source included in the model. A previously unknown gamma-ray source appears, which is spatially coincident with HESS J0632+057. Using pointlike, the best-fit position of this gamma-ray source above 100 MeV is R.A. = 9825, decl. = 581, with a 95% confidence error circle radius of 008 (we shall refer to this source as Fermi J0632.6+0548). HESS J0632+057 is only 21 arcsec away from Fermi J0632.6+0548 and is well within its 95% confidence error circle, which hints at a possible association. By using the best-fit position and assuming a power-law spectral shape (dN/dE = N₀(E/E₀)^−Γ cm⁻² s⁻¹ GeV⁻¹), the gtlike analysis of Fermi J0632.6+0548 resulted in a TS value of 63. We also modeled Fermi J0632.6+0548 by a power law with an exponential cutoff (dN/dE = N₀(E/E₀)^−Γexp(−E/E₀) cm⁻² s⁻¹ GeV⁻¹). The two models are compared using the likelihood ratio test (Mattox et al. 1996). The ΔTS¹⁴ between the two models is less than 9, which indicates that a cutoff is not significantly preferred. The best-fit spectral parameters and corresponding TS values are listed in Table 1, while the SED¹⁵ along with the best-fit power-law model are shown in Figure 3. The flux level of Fermi J0632.6+0548 in the 0.1–300 GeV band is (0.92 ± 0.16 ± 0.08) × 10⁻¹¹ erg cm⁻² s⁻¹, which is consistent with the 3 × 10⁻¹¹ erg cm⁻² s⁻¹ flux upper limit of HESS J0632+057 set by Caliandro et al. (2013).

Zoom In Zoom Out Reset image size

To identify whether Fermi J0632.6+0548 is the GeV counterpart of HESS J0632+057, we carried out an orbital phase-resolved analysis. We adopted the same orbital phase definition of HESS J0632+057 as in Aliu et al. (2014): MJD₀ = 54857 and period P = 315 days. Aliu et al. (2014) reported detailed X-ray/TeV orbital light curves of HESS J0632+057 with aligned enhanced activity in orbital phase 0.2–0.4. However, because of the low statistics, we could not reach the same orbital phase refinement. Thus, in order to search for orbital variability of Fermi J0632.6+0548, we have carried out a binned likelihood analysis in two broad orbital phases, 0.0–0.5 and 0.5–1.0. The two panels of Figure 4 show the TS maps of Fermi J0632.6+0548 in the orbital phases 0.0–0.5 and 0.5–1.0, respectively. Fermi J0632.6+0548 is significantly detected in the phase interval 0.0–0.5 (Figure 4, left panel) with a TS value of 57, an energy flux of (1.43 ± 0.26 ± 0.20) × 10⁻¹¹ erg cm⁻² s⁻¹ and a photon index of 2.55 ± 0.04 ± 0.05 in the 0.1–300 GeV range (Table 1). In the orbital interval 0.5–1.0, the detection of Fermi J0632.6+0548 is less significant (Figure 4, right panel), yielding TS = 23, an energy flux of (0.50 ± 0.21 ± 0.09) × 10⁻¹¹ erg cm⁻² s⁻¹ and a photon index of 2.08 ± 0.12 ± 0.07 in 0.1–300 GeV (Table 1). The flux in the orbital interval 0.0–0.5 is larger than that in the interval 0.5–1.0 at the 98% confidence level, and is also consistent with the orbital variation in X-rays and TeV (Aliu et al. 2014). The spectrum in the orbital interval 0.0–0.5 is steeper than that in the interval 0.5–1.0 at the 99.7% confidence level. A similar steeper-when-brighter behavior was also observed in other gamma-ray binaries (e.g., LSI +61 303, Hadasch et al. 2012; LS 5039, Abdo et al. 2009a), strengthening the association between Fermi J0632.6+0548 and HESS J0632+057.

Zoom In Zoom Out Reset image size

The orbital variations of the flux and spectra are good arguments for a physical association between Fermi J0632.6+0548 and HESS J0632+057, albeit with the caveat of dealing with a dim source that in smaller orbital bins does not reach the detection threshold. We have checked that an orbital light curve produced with a smaller binning (i.e., a binning of 0.1 in phase) yields no significant variation. Finally, adopting the best-fit spatial and spectral model derived from the orbital phase-averaged analysis in Section 4, we calculated the probability of photons coming from Fermi J0632.6+0548 within a radius of 3° using gtsrcprob. A weighted 30-day-binned light curve was produced based on them, and each time bin was exposure-corrected. In order to search for the orbital periodic signal in the light curve, we used the Lomb–Scargle periodogram method (Lomb 1976; Scargle 1982). Power spectra were generated for the light curve using the PERIOD subroutine (Press & Rybicki 1989). No significant periodic signal was discovered in the light curve.

Using nearly nine years of Fermi-LAT data, we have carried out a detailed search for gamma-ray emission from HESS J0632+057, leading to the discovery of a previously unknown gamma-ray source, Fermi J0632.6+0548.

Fermi J0632.6+0548 is spatially coincident with HESS J0632+057, and has a flux level that is consistent with the upper limit previously reported by Caliandro et al. (2013). Based on the orbital phase definition of HESS J0632+057 (Aliu et al. 2014), we searched for orbital variability, finding a flux and spectral change in two broad phase intervals (0.0–0.5 and 0.5–1.0). This variability further hints at a physical association with HESS J0632+057. However, because of the low statistics, neither a significant flux variability in an orbital light curve built with smaller bins, nor the 315-day orbital period in the power spectrum could be detected, leaving the association as likely, but conservatively unconfirmed.

Malyshev & Chernyakova (2016) recently reported a 200–600 GeV detection of HESS J0632+057 at the ∼5σ level during the orbital phases of 0.2–0.4 and 0.6–0.8. For the sake of comparison, we carried out a Fermi-LAT data analysis in the 10–600 GeV range without gating off PSR J0633+0632, similar to what was done by Malyshev & Chernyakova (2016). In the 200–600 GeV range, we confirm that two photons at energies 223 GeV (arrived at mission elapsed time (MET) 301884864, MJD 55404.04) and 578 GeV (arrived at MET 347664434, MJD 55933.89) are spatially consistent with HESS J0632+057. However, no detection of HESS J0632+057 was made during orbital phases 0.2–0.4 and 0.6–0.8 in 200–600 GeV, which is inconsistent with Malyshev & Chernyakova (2016). The inconsistency may be due to the different orbital phase definition adopted: in Malyshev & Chernyakova's work, the orbital phases for the above-mentioned two photons are reported as 0.70 (223 GeV photon) and 0.36 (578 GeV photon). In fact, these authors are using the orbital phase definition from Bongiorno et al. (2011; MJD₀ = 54857, period P = 321 days). These two photons yield the detection of HESS J0632+057 at the ∼5σ level during orbital phases 0.2–0.4 and 0.6–0.8, in the 200–600 GeV range. On the other hand, in our analysis, we used the orbital phase definition from Aliu et al. (2014), which has the same MJD₀ but a refined period (P = 315 days). Correspondingly, the orbital phase of these two photons are calculated as 0.74 (223 GeV photon) and 0.42 (578 GeV photon). Thus, there is only one photon located in these orbital phases, which may explain the non-detection. The different spatial-spectral models used may also lead to the inconsistency: The preliminary seven-year source list was adopted in our analysis together with additional extended templates accounting for gamma-ray contributions from the Rosette Nebula and Monoceros Loop, while Malyshev & Chernyakova (2016) used the second catalog of hard Fermi-LAT Sources (2FHL; Ackermann et al. 2016).

For a constraint on the spectral turnover from the VHE to the HE range, Malyshev & Chernyakova (2016) modeled Fermi-LAT data with a broken power law during the orbital phases 0.2–0.4 and 0.6–0.8 over 10–600 GeV. A 2σ (3σ) limit on the break energy (E_br) was put as E_br = 180–200 GeV (E_br = 140–200 GeV), with a corresponding photon index Γ < 1.2 (Γ < 1.6) below E_br. In the orbital phases 0.2–0.4 and 0.6–0.8, our analysis yielded non-detection, both in the 10–600 GeV range and in the sub energy ranges (10–200 GeV or 200–600 GeV). Thus, further spectral constrains are insignificant.

Fermi J0632.6+0548 is spatially coincident with 3FHL J0632.7+0550, which is a gamma-ray source detected in the Third Catalog of Hard Fermi-LAT Sources (3FHL, Fermi-LAT collaboration 2017). 3FHL J0632.7+0550 is proposed to be associated with HESS J0632+057 and is located within the 95% error circle of Fermi J0632.6+0548. Without gating off PSR J0633+0632, Fermi J0632.6+0548 is detected in the range of 10–600 GeV with TS = 25 and a photon index of 1.74 ± 0.41, which is consistent with the photon index of 3FHL J0632.7+0550, 1.86 ± 0.37, hinting at a possible association.

If the association between Fermi J0632.6+0548 and HESS J0632+057 posed in this paper is real, it will be the first detection of HESS J0632+057 in the HE GeV range, completing its radiation spectrum from radio to TeV. Adopting a distance of 1.4 kpc (Aragona et al. 2010; Casares et al. 2012), the GeV luminosity of HESS J0632+057 is ∼2 × 10³³ erg s⁻¹, about two orders of magnitude lower than those of known gamma-ray binaries (Ackermann et al. 2012b; Hadasch et al. 2012; Caliandro et al. 2013, 2015; Corbet et al. 2016). The radio, X-ray, and TeV luminosities of HESS J0632+057 are also dimmer than known galactic gamma-ray binaries (e.g., Paredes et al. 2007; Skilton et al. 2009; Aliu et al. 2014). Despite the different orbital parameters and multi-wavelength behavior, the companion stars in gamma-ray binaries HESS J0632+057 and LS I +61 303 are very similar. HESS J0632+057 has a B0Vpe star as a companion (MWC 148; Aragona et al. 2010), whereas the spectral type of the companion star in LS I +61 303 is B0Ve (Zamanov et al. 2016). The lower GeV luminosity can be due to a much larger orbital separation (at periastron, the system is twice the size of LS I +61 303, while at apastron it is about seven times bigger, Casares et al. 2012; Zamanov et al. 2016). MWC 148 has a similar radius and mass as LS I +61 303, but its circumstellar disk is about five times larger (Zamanov et al. 2016). The compact object in LS I +61 303 only passes through the outer part of the circumstellar disk at periastron. However, in HESS J0632+057 the compact object goes into the innermost parts and penetrates deeply in the disk during periastron passage (Zamanov et al. 2016), which may lead to large absorption/obscuration effects and explain the low GeV emission.

Detection of HESS J0632+057 with ground-based imaging atmospheric Cherenkov telescopes from hundreds of GeV to several TeV (Figure 3; Aliu et al. 2014) indicates that the VHE spectrum is not a simple extrapolation of the LAT spectra we detected, but likely a different spectral component. Thus, a spectral turnover should exist in Fermi-LAT spectrum. The spectral turnover could arise due to pair production on stellar photons for gamma-rays above ∼50 GeV (Dubus 2006; Sierpowska-Bartosik & Torres 2009), or distinct emission components for HE and VHE spectra. We modeled the HESS J0632+057 with a broken power law in the 0.1–300 GeV range. However, the likelihood ratio test indicates that a broken power law is not significantly preferred over a simple power-law model. Thus, the spectral turnover in Fermi-LAT spectrum could not be explicitly determined because of the low statistics. Based on the SEDs of HESS J0632+057 (Figure 3), we propose the spectral turnover to be above 10 GeV, which is consistent with the estimation by Caliandro et al. (2013). In the well-studied gamma-ray binaries LS 5039 and LS I +61 303, the GeV spectra are best represented by a power law with an exponential cutoff. These spectra do not extrapolate to the VHE range either (Hadasch et al. 2012). Thus, despite its low GeV flux, HESS J0632+057 resembles known gamma-ray binaries and hints at the authenticity of this gamma-ray association.

LS I +61 303 shows 1667-day multi-wavelength super-orbital modulation, which may be due to the quasi-periodic variation of the circumstellar disk (Chernyakova et al. 2012; Li et al. 2012, 2014; Ackermann et al. 2013; Ahnen et al. 2016; Saha et al. 2016). Hosting a similar companion, HESS J0632+057 may also have a multi-wavelength super-orbital modulation. However, its much longer orbital period than LS I +61 303 (26.496 days, Gregory 2002) makes the detection difficult. During the process of publication of this article, an updated orbital solution of HESS J0632+057 appeared (Moritani et al. 2017). With the updated orbital phase definition (MJD0 = 55075.5 and period P = 308 days), we produced an orbital light curve with a binning of 0.1 in phase, yielding no significant variation.

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 Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

We acknowledge the support from The National Key Research and Development Program of China (2016YFA0400800), the grants AYA2015-71042-P, SGR 2014-1073, and the National Natural Science Foundation of China via NSFC-11473027, NSFC-11503078, NSFC-11133002, NSFC-11103020, NSFC-11673013, XTP project XDA 04060604, and the Strategic Priority Research Program The Emergence of Cosmological Structures" of the Chinese Academy of Sciences, Grant No. XDB09000000, as well as the CERCA Programme of the Generalitat de Catalunya. J.L. acknowledges support by the Faculty of the European Space Astronomy Centre. K.S.C. is supported by the GRF Grants of the Government of the Hong Kong SAR under GRF 17302315. We acknowledge the assistance from Dr. M. Kerr on the gamma-ray ephemeris for PSR J0633+0632.