GeV–TeV Counterparts of SS 433/W50 from Fermi-LAT and HAWC Observations

SS 433 is a microquasar in the supernova remnant W50 (see Margon 1984; Fabrika 2004 and references therein). It is likely composed of a ∼20 M_⊙ black hole orbiting a ∼30 M_⊙ supergiant companion with a 13.1 day period. The exotic system is located at a distance of 5.5 kpc (Blundell & Bowler 2004; see discussions about other distance measures in, e.g., Marshall et al. 2013) and about 2° below the Galactic plane. It produces two remarkable jets with kinetic power L_kin ∼ 10³⁹ erg s⁻¹. The jets are heavily loaded with baryons and move at a speed of 0.26c while precessing with a period of 162 days. The angle between the jets and the axis is ∼20°. Other periods are measured but the dynamics is poorly understood (Eikenberry et al. 2001).

Extended X-ray jets are observed on the eastern and western sides (in Galactic coordinates) of SS 433, as shown by the white contours in Figure 1 (Safi-Harb & Ögelman 1997). They interact with and distort the shell of the W50 nebula (Watson et al. 1983; Gregory et al. 1996). which is shown by the gray contours. A set of emission regions, denoted as e1, e2, and e3 centered at 24', 35' and 60' east of SS 433, and w1 and w2 centered at 18' and 31' west of SS 433. have been investigated in detail (Safi-Harb & Ögelman 1997). A bright knot is seen in soft X-rays at e2 (Safi-Harb & Ögelman 1997; Brinkmann et al. 2007), and emission from e1, e2 and w1, w2 is observed in hard X-rays (Safi-Harb & Petre 1999; Moldowan et al. 2005). The X-ray emission can be explained by synchrotron radiation of ∼100–200 TeV electrons in a ∼10 μG magnetic field.

Zoom In Zoom Out Reset image size

Very-high-energy (VHE) γ-ray emission has recently been detected from the SS 433 lobes by the High Altitude Water Cerenkov (HAWC) observatory (HAWC Collaboration et al. 2018). In a data set based on 1017 days of measurements, photons with energies of at least 25 TeV are observed. The TeV hotspots are located close to e1, e2, and w1, with spatially unresolved emission profiles. The flux can be explained by the inverse Compton emission of the same electron population whose synchrotron emission is observed by X-ray telescopes. On the other hand, 40–80 hr observations with the Major Atmospheric Gamma Imaging Cerenkov telescopes and High Energy Spectroscopic System (H.E.S.S.) report no evidence of γ-ray emission between a few hundred GeV and a few TeV from the jet termination regions, or from the central binary (MAGIC Collaboration et al. 2018). A similar upper limit is reported by the Very Energetic Radiation Imaging Telescope Array System (VERITAS; Kar & VERITAS Collaboration 2017).

Searches (Bordas et al. 2015; Rasul et al. 2019; Sun et al. 2019; Xing et al. 2019) have been made for a GeV counterpart in the data observed by the Fermi Large Area Telescope (LAT) (Atwood et al. 2009). Analysis of LAT data in this region faces two complications. A point source FL8Y J1913.3+0515 from the preliminary LAT 8 yr point source list⁴ (FL8Y) is tagged as possibly associated with W50. It is no longer a source in the 4FGL catalog due to different emission models (The Fermi-LAT Collaboration 2019). Different conclusions about detection of SS 433 have been reached depending on whether this source is included in the background model. In addition, analysis of the emission profile and spectrum of the SS 433 region is heavily impacted by the nearby pulsar PSR J1907+0602 (4FGL J1907.9+0602), which is not suppressed via selection on the rotational phase in the abovementioned works. Due to these difficulties, whether the jets of SS 433 shine at GeV energies is unknown.

Since the source is marginally significant in γ-rays yet close to the bright Galactic plane, it is difficult to study solely with the GeV or the TeV measurements. Here we jointly analyze an ROI surrounding SS 433 observed by LAT and HAWC. A simultaneous fit to the 100 MeV–100 TeV data directly addresses the question of whether γ-ray emission over six decades can be produced by a common cosmic-ray population inside the SS 433 lobes. We explain the methods in Section 2, present the results of the LAT analysis in Section 3.1 and that of the joint analysis in Sections 3.2 and 3.3. We discuss immediate implications of this analysis in Section 4.

Our analysis uses 10.5 yr of Fermi-LAT data and 1017 days of HAWC data. Details on the LAT and HAWC analyses, as well as background sources in each band, are presented in Appendices A and B. The setup of a joint analysis of the LAT and HAWC data based on the 3ML framework (Vianello et al. 2015) is presented in Appendix C. Here we describe the procedure for the joint analysis.

We first build a source model to describe the broadband γ-ray emission of SS 433. Three types of models are considered.

• 1.  
  γ-rays follow a power-law spectrum, ${dN}/{{dE}}_{\gamma }={K}_{\gamma }{\left({E}_{\gamma }/{E}_{\gamma ,\mathrm{piv}}\right)}^{-{\alpha }_{\gamma }}$.
• 2.  
  γ-rays follow a log parabola spectrum, ${dN}/{{dE}}_{\gamma }={K}_{\gamma }{\left({E}_{\gamma }/{E}_{\gamma ,\mathrm{piv}}\right)}^{-{\alpha }_{\gamma }-{\beta }_{\gamma }\mathrm{log}({E}_{\gamma }/{E}_{\gamma ,\mathrm{piv}})}$.
• 3.  
  Electrons are injected with a rate $d\dot{N}/{{dE}}_{e}={Q}_{e}\,{\gamma }_{e}^{-{\alpha }_{e}}\exp \left(-{E}_{e}/{E}_{e,\max }\right)$, where γ_e = E_e/m_e c² is the Lorentz factor of an electron. They up-scatter the cosmic microwave background (CMB) and infrared photons in W50 to γ-rays through the inverse Compton process, and produce synchrotron emission in a magnetic field B.

Models I and II are simple descriptions of γ-ray spectral shapes. Model III is physically motivated. The cooling time of VHE electrons, ${t}_{e,\mathrm{cool}}=2.5{(B/10\mu {\rm{G}})}^{-2}{({\gamma }_{e}/{10}^{8})}^{-1}$ kyr, is much less than the source age, ${t}_{\mathrm{age}}\sim 30\,\mathrm{kyr}$ (Fabrika 2004). In order to take into account the effects of cooling, we solve a transport equation for each set of parameter values. Details about the cooled electron spectrum can be found in Appendix D. Once a steady-state cosmic-ray spectrum is obtained, the γ-ray flux is calculated using the radiative functions of the naima package (Khangulyan et al. 2014; Zabalza 2015).

The models are then converted into data space and compared to observations through the joint analysis framework (see Appendix C). Since the GeV and TeV observations are carried out independently, a total log likelihood is evaluated by summing the log likelihoods from the GeV and TeV analyses. The total likelihood is then maximized by adjusting model parameters to obtain the best-fit source model. Finally, the likelihood test statistic (TS) of a target source is computed as twice the difference of the log likelihoods of the data given the models with and without the source.

3.1. LAT Analysis Results

Here we present results from the LAT-only analysis. The method is detailed in Appendix A. The main difference between our analysis and previous works is that we use a data set for which the PSR J1907+0602 is gated off, that is, the arrival times of photons are phase-folded, and the photons that arrive during the pulsar's pulse peak are removed (see J. Li et al. 2019, in preparation for details). Throughout the work we use the 4FGL catalog and the corresponding diffuse emission models to model background sources.

The significance of the residual γ-ray excess from the SS 433/W50 region in the LAT data between 100 MeV and 300 GeV is shown in the left panel of Figure 1. The most statistically significant excess is near the location of FL8Y J1913.3+0515. We call this excess J1913+0515 to differentiate it from FL8Y J1913.3+0515. J1913+0515 is at the boundary of W50, well outside the extended X-ray jets. When describing the spectral energy distribution (SED) with a power-law function, we obtain ${K}_{\gamma }=1.5\times {10}^{-12}\,{\mathrm{MeV}}^{-1}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$, ${E}_{\gamma ,\mathrm{piv}}=0.9\,\mathrm{GeV}$, and α_γ = 2.4. The best-fit location is very close but slightly different from the location listed in the FL8Y catalog. The test statistic of the source is TS = 32.8 using a regular likelihood and TS = 25.7 using a weighted likelihood (The Fermi-LAT Collaboration 2019) that takes into account estimated systematic uncertainties in the diffuse emission models. As the results obtained by the two methods are similar, we use a regular (i.e., unweighted) likelihood in the rest of our analyses.

A sub-threshold (TS < 25) excess is evident at the northeastern side of J1913+0515. Because it is spatially close to the TeV excess in the eastern lobe, we refer to it as the "eastern hotspot." It is not significant in the LAT data and has TS = 5.0 when J1913+0515 is included in the background model. The excess is due to several high-energy photons at ∼20–50 GeV, as shown by the SED in the top panel of Figure 2.

Zoom In Zoom Out Reset image size

In the western lobe, a sub-threshold excess is found between w1 and w2 (which we refer to as the "western hotspot" below). The excess region partially overlaps with the X-ray jets and touches the boundary of W50. We find a TS of 16.1 for the western hotspot when adding it to the baseline model. Its spectrum can be described by a power law of index 2.3, as shown in the bottom panel of Figure 2.

When including the potential sources in the baseline model simultaneously, we obtain TS ∼ 5 and TS ∼ 10 for the eastern and western hotspots, respectively. The fit results are summarized in Table 2 in Appendix A. Neither of the "hotspots" is statistically significant in the LAT observations but they are evident in the joint analysis, as will be shown in Section 3.3.

3.2. J1913+0515 and the TeV Emission

3.3. Joint Analysis Results

Because of its proximity and exotic structure, SS 433 has been one of the most observed Galactic high-energy sources for over 40 yr. Nonetheless, no consensus has been reached about what happens inside the SS 433/W50 complex. The detection of multiple tens of TeV photons from the object confirms the existence of particles at extreme energies, but deepens the question of why lower-energy γ-rays have not been observed. By jointly analyzing an ROI measured by both the Fermi-LAT and the HAWC Observatory, we find common sites of GeV and TeV γ-ray emission inside the SS 433 lobes. The SED is consistent with inverse Compton emission of an electron population accelerated by the jets but quickly cooled due to synchrotron radiation in a magnetized environment. We use a data set that suppresses emission by a nearby pulsar that highly impacts previous analyses. Our joint analysis concludes that the GeV point source J1913+0515 located at the boundary of W50 is an unlikely counterpart to the TeV emission. This addresses the dilemma encountered by Xing et al. (2019), Sun et al. (2019), and Rasul et al. (2019).

This is the first joint ROI analysis across γ-ray observatories, to our knowledge. Using a framework built on individual data analysis toolkits from γ-ray observatories, we have shown that such an approach is feasible. The joint analysis is designed to study shared properties of sources of γ-rays over a very broad spectrum. It maximizes the usage of data, including sub-threshold information, and is more powerful than simply combining results from each experiment. Future data from HAWC, especially with refined angular resolutions (HAWC Collaboration et al. 2019), will help improve the understanding of γ-ray emission from the western lobe. Future observations by IACTs, as well as by X-ray and radio telescopes, at the revised source locations will help further constrain the emission models.

The GeV–TeV spectrum can be explained as inverse Compton scattering by X-ray synchrotron-emitting ∼100 TeV electrons. Three scenarios can be entertained to account for the acceleration. The first is to invoke direct, diffusive shock acceleration of the electrons at the termination shocks of the precessing jets launched by the accretion disk. If the post-shock field strength is ∼10 μG then acceleration to these energies is possible, though only a small electron power ∼10³⁴ erg s⁻¹ is needed to account for the γ-rays. Second, ∼5 PeV protons may also be shock-accelerated. Their primary radiative loss could be due to Bethe–Heitler pair production on ∼2 eV optical photons from SS 433 with a cross section of a few millibarn. Maximum electron or positron energies ∼100 TeV are just possible. However, in order to account for the γ-ray power, a proton power ∼10³⁸ erg s⁻¹ is necessary. The third possibility is that a hitherto unobserved, ultrarelativistic, electromagnetic jet is formed by the spinning black hole. Such a jet can create an EMF ∼ 100(L_jet/10³⁹ erg s⁻¹)^1/2 PV that suffices to accelerate the emitting particles. Finally, if the gas density in the lobes is high ($\gtrsim 1\,{\mathrm{cm}}^{-3}$), pion production can also make a contribution to the γ-ray flux. High-energy neutrinos would be produced simultaneously and could be measured by IceCube (IceCube Collaboration et al. 2019).

Future multi-messenger observations of the γ-ray emission regions have the potential to discriminate between these scenarios.

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. This work was performed in part under DOE Contract DE-AC02-76SF00515.

We thank Henrike Fleischhack, Colas Riviére, and Giacomo Vianello for their help with the usage of the 3ML and HAWC-HAL packages.

We analyze 10.5 yr of Pass 8 data⁵ taken between 2008 August 4 15:43:36 UTC and 2019 January 28 00:00:00 UTC using version 0.17.4 of fermipy⁶ and version ScienceTools-11-04-00 of Fermitools.⁷ We define the ROI as the 15° × 15° region in Galactic coordinates centered at SS 433 (l = 39.69, b = −2.24). γ-ray events with energies between 100 MeV and 300 GeV are selected. Other event selection criteria include a 90° zenith cut and a filter expression of "DATA_QUAL > 0 && LAT_CONFIG == 1" which are standard quality criteria recommended by the Fermi Science Support Center. We use the P8R2_SOURCE event selection, and the corresponding P8R2_SOURCE_V6 LAT instrument response functions. Unlike previous works analyzing the SS 433 region (Rasul et al. 2019; Sun et al. 2019; Xing et al. 2019), here we use the latest LAT 8 yr Point Source Catalog (The Fermi-LAT Collaboration 2019), together with the corresponding Galactic diffuse model gll_iem_v07.fits and the isotropic diffuse model. The 4FGL catalog and the updated diffuse emission model turn out to considerably impact the analysis of 100–300 MeV photons in this region, comparing to the FL8Y catalog.

As the nearby pulsar PSR J1907+0602 is very bright in the GeV band, we use the same method from J. Li et al. (2019, in preparation) to suppress the pulsar emission. The same pulsar ephemeris is adopted in pulsar gating, which amounts to 44% of the observing time. The exposure is scaled accordingly.

There are 33 4FGL sources within 5° of SS 433 and 61 4FGL sources within 8°. The baseline ROI analysis is performed using the fermipy.job sub-package fermipy-analyze-roi. The optimized model is referred to as "baseline model."

A significance map of the residual γ-ray excess is shown in the left panel of Figure 1. The color scale corresponds to the square root of the TS when there is a new point source at a given location, in addition to known sources from the 4FGL catalog, the Galactic diffuse emission and the isotropic diffuse emission. The test point source is assumed to have an E⁻² spectrum and the TS is evaluated for each location on a grid with 01 × 01 spacing.

After setting up the baseline model, we add J1913+0515, the eastern and western hotspots, to the background model and refit the new model to the data. The new model is refit to the data using GTAnalysis.optimize, which fits sources in order of their fluxes. The fit returns TS values of 26.1, 5.0, and 9.6 for the three candidate sources, respectively. We also tested an alternative fitting method, where we fixed the parameters of background sources to their best-fit values, and vary only the normalization of new sources using GTAnalysis.fit. This approach returned TS values of 28.0, 4.9, and 10.4 for the three. Since the difference of the results from the two fitting methods is minor while the latter is much more efficient in computation time, we use the second approach to calculate the LAT likelihood in a joint analysis. The fit results are summarized in Table 2.

Although MGRO J1908+06 is one of the brightest TeV sources, an extended source at its location is not significant in the LAT data. We thus do not include it in our background model.

The data points for the SED were obtained by binning the spectrum with 2 bins per decade in energy and performing a likelihood analysis in each energy bin.

In the TeV band, we analyze the public data⁸ from the HAWC Observatory (HAWC Collaboration et al. 2018). The data set contains 1017 days of γ-ray events collected between 2014 November 26 and 2017 December 20. The reconstruction of the arrival direction of primary γ-rays is based on the relative arrival times of photoelectron hits detected by the photomultipliers inside the water Cerenkov detectors. (This kind of reconstruction is referred to as the nhit method.) Angular resolution from the nhit analysis ranges from around 1° below 1 TeV to <02 above 10 TeV.

We adopt the same ROI as in HAWC Collaboration et al. (2018), which is defined to be a semicircular region with a radius of 25 centered on the position of MGRO J1908+06 (as shown in Extended Data Figure 1 of HAWC Collaboration et al. 2018). By masking the sources close to the Galactic plane, the contamination from the Galactic diffuse emission is significantly reduced.

Three sources remain in the ROI: MGRO J1908+06, the eastern and the western hotspots in the SS 433 lobes. Following HAWC Collaboration et al. (2018), we use the electron diffusion model to describe the spatial morphology of MGRO J1908+06. Other spatial models with Gaussian and power-law radial profiles lead to similar results.

Unlike the analysis in HAWC Collaboration et al. (2018), which is based on the HAWC analysis framework AERIE, here we redo the analysis using the HAWC Accelerated Likelihood (HAL) framework. HAL provides faster convolution with the detector response functions, which is needed for the joint analysis in our work. We have confirmed that the two analysis frameworks lead to results that are consistent at the 1% level.

The workflow of a joint analysis is diagrammed in Figure 3. The joint analysis is implemented in the Multi-Mission Maximum Likelihood framework (3ML; Vianello et al. 2015).⁹ 3ML is a data analysis architecture that converts emission models for an ROI into data spaces for specified instrument(s), and compares the model predictions to the corresponding data based on the likelihood formalism. The Fermi-LAT module of the package provides a wrapper of fermipy (Wood et al. 2017) and the Fermi Science Tools.¹⁰ The HAWC module links to the HAWC analysis tools including the Accelerated Likelihood (HAL) framework.¹¹

Zoom In Zoom Out Reset image size

Since an ROI analysis is not implemented in 3ML, we use fermipy to perform a "baseline" analysis externally (see Appendix A) and use that as a starting point for 3ML analysis. Meanwhile, the source model and a model of MGRO J1908+06 are passed to the HAWC plugin (see Appendix B). In this way the contribution of background sources is taken into account properly.

The cooling of relativistic electrons in the lobes of SS 433 can be described by a transport equation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {N}_{e}}{\partial t}+\displaystyle \frac{\partial }{\partial {\gamma }_{e}}\left[\dot{{\gamma }_{e}}\,{N}_{e}({\gamma }_{e},t)\right]={Q}_{e}({\gamma }_{e},t),\end{eqnarray} \tag{ 1 }$$

where ${\dot{\gamma }}_{e}=-4/3\,{\gamma }_{e}^{2}\,c{\sigma }_{T}({u}_{B}+{u}_{\gamma })/({m}_{e}\,{c}^{2})\equiv -\nu \,{\gamma }_{e}^{2}$ is the energy loss rate due to inverse Compton and synchrotron emission. σ_T is the Thomson cross section. N_e and Q_e are the spectrum and injection rate of electrons, respectively. u_B = B²/(8π) and ${u}_{\mathrm{CMB}}=0.26\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$ are the energy density of magnetic field and the CMB. We also adopt a FIR background at 20 K with ${u}_{\mathrm{FIR}}=0.3\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$ motivated by the dust emission in the solar neighborhood (Vernetto & Lipari 2016). Background photons with higher energies are not important due to the Klein–Nishina effect. Due to its much lower energy density, the synchrotron radio emission of W50 and the lobes is not expected to contribute significantly to the cooling of electrons or production of high-energy γ-rays. In Equation (1) we have ignored the diffusion of electrons, as it is a slower process than cooling for TeV electrons, and also because doing so saves computing time. Assuming that electrons are injected with a simple power-law spectrum constantly over time, ${Q}_{e}({\gamma }_{e},t)={Q}_{e,0}\,{\gamma }_{e}^{-\alpha }$, the solution of Equation (1) can be written as

$$\begin{eqnarray}&&N({\gamma }_{e},t)=\displaystyle \frac{{Q}_{e,0}}{{\gamma }_{e}^{2}}{\int }_{{t}_{\min }}^{t}{{dt}}_{i}{\left[{\gamma }_{e}^{-1}-\nu (t-{t}_{i})\right]}^{\alpha -2}\end{eqnarray} \tag{ 2 }$$

where ${t}_{\min }=\max [0,t-{\nu }^{-1}({\gamma }_{e}^{-1}-{\gamma }_{e,\max }^{-1})]$ is the earliest time that an electron with γ_e can be injected and still not cooled after time t, and γ_e,max is the maximum electron energy.