HAWC and Fermi-LAT Detection of Extended Emission from the Unidentified Source 2HWC J2006+341

2HWC J2006+341 was discovered by the High Altitude Water Cerenkov observatory (HAWC; Abeysekara et al. 2017a) in the Cygnus Region of the Galaxy at the Galactic coordinates (l, b) = (7133, 116). No other TeV sources are known to exist nearby, and no supernova remnants (SNR) have been detected within 1 degree of 2HWC J2006+341 (Green 2014). The pulsar PSR J2004+3429, located 04 away from the position of 2HWC J2006+341, is the nearest pulsar found in the ATNF catalog (Manchester et al. 2005). With a characteristic age of 18 kyr, the estimated distance and spindown power of this pulsar are 10.8 kpc and 5.8 × 10³⁵ erg s⁻¹, respectively (Barr et al. 2013).

In the GeV range, two point sources are found near 2HWC J2006+341 in the Fermi Large Area Telescope (LAT) fourth source catalog (4FGL; Abdollahi et al. 2020). 4FGL J2004.3+3339, ∼064 away, is associated with a binary system (powering the nebula G70.7+1.2, Kulkarni et al. 1992). 4FGL J2005.8+3357, with a detection significance in the 4FGL catalog of 5.8σ, is located about 02 from 2HWC J2006+341 and has no association at other wavelengths.

After the discovery of 2HWC J2006+341 the MAGIC and Fermi-LAT collaborations analyzed observations of the region to search for gamma rays (Ahnen et al. 2019). Only upper limits were placed on the emission. Given the initially reported point-source morphology for 2HWC J2006+341 by HAWC, no search for extended sources with radii above 02 was carried out by these authors.

This work expands our preliminary studies presented by Araya & HAWC Collaboration (2019) with a more detailed analysis of HAWC and Fermi-LAT data. Motivated by the possible presence of extended TeV emission in the 2HWC J2006+341 region in more recent and deeper HAWC maps, we searched for a corresponding signal in the GeV range with publicly available LAT data. The data from both observatories reveal emission with similar morphologies and consistent spectra, thus we believe that the GeV and TeV photons are produced by a single unidentified source and we discuss several scenarios for its nature.

2.1. Fermi-LAT

The LAT on board the Fermi satellite is a converter/tracker telescope capable of detecting gamma rays in the energy range between 20 MeV and ≳1 TeV (Atwood et al. 2009). We analyzed Pass 8 LAT data from 2008 August to 2020 August with fermitools version 1.2.23 and fermipy version 0.19.0. Events with good quality in the SOURCE class were selected (with the parameter evclass = 128), including both back and front converting event types combined (with the parameter evtype = 3). The corresponding detector response P8R3_SOURCE_V2 was used. To reduce the gamma-ray background produced by the Earth's limb, only events having zenith angles lower than 90° were considered. Events were binned with a spatial scale of 005 per pixel and using 10 bins per decade in energy for exposure calculation. The region of interest (ROI) analyzed had a radius of 15° and it was centered at the coordinates R.A. = 3015, decl. = 340.

The energy range considered for spectral analysis was 1–500 GeV while the analysis of the morphology was carried out in the 5–500 GeV range to take advantage of the improved LAT resolution at higher energies. Sources found in the 4FGL catalog, located within 20° of the center of the ROI, were included in the analysis (162 in total). The diffuse Galactic emission and the isotropic emission (including the residual cosmic-ray background) were modeled with the standard files gll_iem_v07.fits and iso_P8R3_SOURCE_V2_v1.txt, respectively, provided by the LAT team.³⁴ The normalizations of the spectra of the 53 sources located within 10° of the center of the ROI were left free to vary while their spectral shapes were fixed to their catalog values. All the spectral parameters of sources located farther away were fixed. The best-fit values of the free parameters and the optimization of the spectrum and morphology of the sources were obtained through the maximum likelihood technique (Mattox et al. 1996), which maximizes the probability for the model to explain the data. This procedure allows for the estimation of the detection significance of a new source by calculating the test statistic (TS), defined as −2 × log $({{ \mathcal L }}_{0}/{ \mathcal L })$, with ${ \mathcal L }$ and ${{ \mathcal L }}_{0}$ being the maximum likelihoods for models with and without the additional source, respectively. The weak source 4FGL J2005.8+3357, located in the region of 2HWC J2006+341 and having no association, was removed from the model in order to carry out a more detailed study of the emission.

In the first part of the analysis above 5 GeV, a search for extended emission in the region was carried out using two different morphological templates: a uniform disk and a symmetric Gaussian. A simple power-law spectral shape of the form

$$\begin{eqnarray*}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}\cdot {\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-\gamma },\end{eqnarray*}$$

where E₀ is a (constant) scale factor, was assumed (and later proved to be appropriate, see below). The location and size of the extended templates, as well as the spectral parameters, were fitted to maximize the likelihood. In order to compare models and since the models used are nonnested (it is not possible to go from one model to the other with a variation of the parameters) we applied the Akaike Information Criterion (AIC; Akaike 1974). This is calculated as AIC = $2k-2\mathrm{ln}({ \mathcal L })$, where k is the number of parameters and ${ \mathcal L }$ the maximum likelihood. Given a set of models, the one that better describes the data is the one that minimizes the AIC.

Once the best-fit morphology of the emission is determined, the spectra of the sources was obtained in the 1–500 GeV energy range. The spectral shape of the source of interest was changed for comparison and the new free parameters were fitted. The alternative spectral shape used was a log-parabolic power-law function (or log-parabola) of the form

$$\begin{eqnarray*}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}\cdot {\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-\alpha -\beta \mathrm{ln}(E/{E}_{0})}.\end{eqnarray*}$$

The difference in TS values was used to indicate if there are significant deviations from a pure power law. In order to investigate the existence of additional sources in the region we carried out a systematic search for new point sources having a TS > 25 and added them to the model. This was done before both the morphological and spectral analyses.

2.2. HAWC

HAWC is a ground-based air shower detector array located at a latitude of ∼19° N in Mexico. Employing 300 water tanks of 7.3 m in diameter and 4.5 m in height each, it can scan about two-thirds of the sky (8.4 sr, from −26° to 64° decl.) detecting very-high-energy gamma rays with energies from hundreds of GeV to >100 TeV and a >95% duty cycle Abeysekara et al. (2017a, 2017b, 2019). The angular resolution depends on the source location in the sky and the energy of the events, and varies from 1° to 02 (Abeysekara et al. 2017b), which is similar to the resolution of the LAT (Abdollahi et al. 2020). The data used, with a livetime of 1038.7 days, were binned using the fraction of detectors triggered by an air shower, and the ground parameter (as defined in Abeysekara et al. 2019) was used to estimate the energy of the events. Recently optimized gamma/hadron separation cuts and reconstruction were used, which were found to improve the energy and angular resolutions with respect to previous analyses (Abeysekara et al. 2019). The Multi-Mission Maximum Likelihood framework³⁵ (3ML; Vianello et al. 2015) was used to perform the likelihood fits and determine the morphology and spectrum of the source. The HAWC accelerated likelihood³⁶ (HAL) plugin was used.

The ROI analyzed was a circular region with a radius of 7° around the cataloged position of 2HWC J2006+341 (Abeysekara et al. 2017a). Besides the source of interest, the model contained the pulsar wind nebula (PWN) 2HWC J2019+367, located ∼36 from the center of the ROI, and a background template of uniform brightness covering the entire ROI to account for emission from any possible unresolved sources. Changing the model of the gamma-ray background in the region would affect the parameters of the source of interest, and thus we considered an additional source of systematic uncertainty associated to this model. We calculated the deviations of the spectral parameters obtained when replacing the background template with two alternative models and optimizing the likelihood. In one alternative no gamma-ray background was included and for the other one we used a 2D Gaussian background centered at the Galactic plane. We added the differences between the resulting parameter values and the original ones in quadrature to the other systematic errors. Similarly to the analysis of LAT data, maximum likelihood fits were done using a point-source hypothesis as well as the uniform disk and symmetric Gaussian templates, assuming a simple power-law spectrum for the emission. Once the best-fit morphology is found, a fit using a log-parabolic spectrum was done to search for deviations from the simple power-law (spectral curvature). Figure 1 shows a significance map of the 2HWC J2006+341 region calculated with the maximum likelihood ratio and an extended source with a power-law spectrum with an index γ = 2.7, which results from the initial fit. The results of the morphological analysis are also shown in the figure.

Zoom In Zoom Out Reset image size

In the LAT data no point sources having a TS > 25 were found in the region of 2HWC J2006+341 above 5 GeV. However, significant emission was detected using extended templates that gave TS > 40 for both the disk and Gaussian above 5 GeV. The likelihood ratio between the best-fit extended model and the best-fit point-source model was such that TS_ext > 20 in both cases. Similarly, in HAWC data the extended templates result in significant detections. Table 1 gives the fitted extensions, ΔAIC and TS values of the fits for each analysis. Here, ΔAIC = AIC_m − AIC_min is the difference in AIC between each model m and the one that minimizes the AIC (ΔAIC = 0 for the best available model). In the analysis of HAWC data the Gaussian template provided the best description compared to the others. The fit with a Gaussian morphology also produced smaller residuals in HAWC data. For these reasons this template was chosen to represent the emission in the region of 2HWC J2006+341 at TeV energies. The locations of the Gaussian centers that maximized their respective likelihoods (with their respective 1σ-level statistical uncertainties) are given by R.A. = 30155 ± 018, decl. = 3435 ± 016 in the HAWC analysis, and R.A. = 30141 ± 015, decl. = 3433 ± 016 in the fits to the LAT data. The locations are compatible within uncertainties.

Table 1.  Results of the Independent Morphological Analyses of the Fermi-LAT and HAWC Data Sets

Data Set	Spatial Model	Fitted Sizea (°)	R.A. (°)	Decl. (°)	ΔAIC	TSb
Fermi-LAT
	Disk	${0.85}_{-0.05}^{+0.04}$	301.50 ± 0.06	34.47 ± 0.08	0	44.6
	Gaussian	${0.59}_{-0.11}^{+0.10}$	301.41 ± 0.15	34.33 ± 16	3.6	41.0
HAWC
	Point source	⋯	301.1 ± 0.07	34.16 ± 0.05	23.1	24.0
	Disk	0.88 ± 0.10	301.6 ± 0.15	34.5 ± 0.15	7.2	41.9
	Gaussian	0.72 ± 0.14	301.55 ± 0.18	34.35 ± 0.16	0	49.0

Notes.

^aRadius for the disk and σ for the Gaussian. ^bThe LAT TS values were obtained with events above 5 GeV.

Download table as:  ASCIITypeset image

In the analysis of LAT data both extended morphologies tested represented a similar improvement with respect to the null hypothesis. The Gaussian template was chosen for the rest of the analysis for a better direct comparison with the HAWC results. The effect of changing the morphological model was considered a source of systematic uncertainty. The other source of systematic uncertainty in the LAT data was estimated by propagating the uncertainty in the effective area onto the spectral index and flux. A set of bracketing response functions were used for this purpose as recommended by Ackermann et al. (2012). As a consistency check, we repeated the entire analysis after replacing the diffusion model with a different version, given by the file gll_iem_v06.fits, and obtained similar results. The source extension was also confirmed with TS_ext > 28 after artificially varying the best-fit value of the diffuse Galactic normalization by ±6%, an estimate of the systematic uncertainties in the diffuse emission model used in the past by the LAT team (e.g., Abdo et al. 2010). Systematic uncertainties in the HAWC data result from uncertainties in the modeling of the detector. For the spectral parameters measured by HAWC, these uncertainties were estimated as explained in Abeysekara et al. (2019). New fits were done with alternative instrument response files changing key parameters used to model the instrument. The differences between the resulting spectral parameters and the nominal values were used to estimate the systematic uncertainties.

3.1. Spectrum of the HAWC Source

3.2. Spectrum and Significance of the LAT Source in the 1–500 GeV Range

3.3. GeV–TeV Spectral Energy Distribution

Since the morphologies of the LAT and HAWC sources are consistent with each other we attribute the emission to a single new extended gamma-ray source. Figure 2 shows the spectral energy distributions (SED) obtained from both data sets and the models discussed in Section 4.

Zoom In Zoom Out Reset image size

The analysis of the LAT data shows that there is an extended source of GeV emission in the region of 2HWC J2006+341. The size and location of this source are consistent with those of the TeV emission revealed in the detailed analysis of the HAWC data presented here. This result and the fact that the SEDs measured by both instruments show comparable fluxes, as can be seen in Figure 2, are evidence that the particles responsible for the GeV and TeV gamma rays come from the same source. Furthermore, the overall TeV spectrum is not an extrapolation of the power-law spectrum seen with the LAT, which is harder, indicating the presence of spectral softening.

Regarding the origin of the gamma rays, both an SNR and a PWN could in principle be responsible, as they typically produce extended gamma-ray emission (e.g., Acero et al. 2016; H.E.S.S. Collaboration et al. 2018a; Abdollahi et al. 2020). In order to account for the SED seen in Figure 2, simple emission mechanisms known to exist within these objects can be explored to get an idea of the energetics in the particles producing the radiation. For this purpose, fits to the data points shown in Figure 2 were done with the Markov Chain Monte Carlo fitting routines of naima, a package for the calculation of nonthermal emission from relativistic particles (Zabalza 2015). Using a particle energy distribution that is a power law with an exponential cutoff, the data were fit to fluxes resulting from both leptonic and hadronic scenarios. In the leptonic scenario the particles interact with the cosmic microwave background (CMB) bringing these photons to gamma-ray energies through inverse Compton (IC) scattering (for details on the implementation of the calculation, see Khangulyan et al. 2014). In the hadronic scenario, relativistic protons produce gamma rays through collisions with ambient protons (see Kafexhiu et al. 2014). The best-fit models obtained in both scenarios are plotted in Figure 2. The resulting AIC values are 13.2 (for the leptonic model) and 15.4 (for the hadronic scenario).

PWNe are believed to produce gamma rays mainly from electrons and positrons interacting with soft background light, and therefore because these particles cool off efficiently as they propagate away from the source the size of PWNe resulting from IC scattering usually decreases with increasing energy. This prediction is consistent with observations of HESS J1825-137 and Geminga (Aharonian et al. 2006; Di Mauro et al. 2019; Principe et al. 2019). However, the results of the morphological analysis of 2HWC J2006+341 presented here reveal extensions at GeV and TeV energies that are consistent within uncertainties. This could disfavor the PWN origin for the emission and deeper observations will be required to confirm this. In the leptonic IC-CMB scenario, the resulting spectral index and cutoff energy of the particle distribution, as well as the total energy in the particles above 1 GeV (and their corresponding 1σ statistical uncertainties) are ${2.38}_{-0.19}^{+0.10}$, $({42}_{-15}^{+19})$ TeV, and ${({8.8}_{-5.1}^{+6.6})\times {10}^{47}(d/1\mathrm{kpc})}^{2}$ erg, respectively. Here, d is the source distance. As noted in Section 1 the only pulsar known near the gamma-ray source is PSR J2004+3429, with an estimated distance of 10.8 kpc. Assuming the same distance to 2HWC J2006+341 the required energy in the leptons would be in the range 10⁴⁹–10⁵⁰ erg, considering the uncertainty. A rough estimate of the energy injected by the pulsar can be obtained by the product of the spindown power with its characteristic age, which yields ∼4 × 10⁴⁷ erg. Assuming a braking index for the pulsar of 2, rather than the default value of 3, could change this estimate by a factor of a few (Gaensler & Slane 2006), but the energy in the particles would be two orders of magnitude above this value, making it difficult to reconcile the energetics. If additional photon fields are considered in the calculation of the IC fluxes the required energy in the leptons is still very large. For temperatures and densities of far-infrared dust emission and near-infrared stellar emission seen in the solar system the total required energy decreases by a factor of ∼3 with respect to the IC-CMB scenario, which is not enough.

Assuming the same distance of 10.8 kpc as PSR J2004+3429 for the gamma-ray source, the luminosity in the 2–13 TeV range would be ∼2 × 10³⁵ erg s⁻¹, or about 7 times the luminosity of the Crab Nebula in the same energy range. This TeV luminosity for 2HWC J2006+341 and the characteristic age of pulsar PSR J2004+3429 would be similar to the corresponding values for the firmly identified PWN HESS J1825-137 and its associated pulsar. However, the spindown power of PSR J2004+3429 is almost 5 times lower than that of the pulsar associated to HESS J1825-137. This would make 2HWC J2006+341 fall outside of the trends seen for other PWNe regarding their TeV luminosities (H.E.S.S. Collaboration et al. 2018b). Furthermore, the feature that would make 2HWC J2006+341 a very unusual PWN would be its large physical size. An extension of ∼17 in the sky corresponds to a size of more than 300 pc for a distance to PSR J2004+3429 of 10.8 kpc. It then becomes difficult to propose PSR J2004+3429 as the source of the gamma-ray emission from 2HWC J2006+341. Although it is possible in principle for 2HWC J2006+341 to be a PWN produced by an unknown pulsar located at a closer distance, the similar sizes at GeV and TeV energies found here could be in conflict with the PWN scenario. Deeper observations with more statistics and a detailed study of the environment at the location of 2HWC J2006+341 are necessary to confirm or reject this hypothesis.

The discovery by HAWC of extended TeV emission around the Geminga and PSR B0656+14 pulsars (Abeysekara et al. 2017c) points to the possible existence of a population of middle-aged pulsars producing "TeV halos" in the surrounding interstellar medium (e.g., Linden et al. 2017). The possibility that a previously undetected pulsar is responsible for the gamma rays from 2HWC J2006+341 cannot be discarded. However, at GeV energies the likely counterpart of the Geminga halo is much more extended than the TeV emission (Di Mauro et al. 2019), which is not seen for 2HWC J2006+341.

Another possibility for the origin of 2HWC J2006+341 is the SNR shell scenario. SNRs are mainly seen in radio observations. However, inspection of radio images from the 1.4 GHz NRAO VLA Sky Survey (NVSS; Condon et al. 1998) and the Westerbork Northern Sky Survey (WENSS; Rengelink et al. 1997) show no obvious hints of emission. The known SNR that is closest in the sky to 2HWC J2006+341 is ∼18 away (G69.7+1.0, Green 2019). The lack of a known counterpart at lower wavelengths is not necessarily evidence against the SNR scenario. There are several SNR shell candidates such as HESS J1614-518 that have only been detected at TeV energies (H.E.S.S. Collaboration et al. 2018c). Using the same simple leptonic scenario described above, the total energy in the particles is reasonable for a large range of possible distances to the source. The typical kinetic energy in the shocks of SNRs that is available for accelerating particles is of the order of E_SN = 10⁵¹ erg. The required total energy in the leptons above 1 GeV, ∼9 × 10⁴⁷(d/1 kpc)² erg, is not problematic for a large range of possible distances to the source. Given the extension of 2HWC J2006+341 and that the typical diameters of SNRs are of tens of parsecs, the maximum distance to the source would likely be several kiloparsecs. For a source distance of 1 kpc, for instance, the diameter of the gamma-ray emission region would be  30 pc and the energy content in the leptons (above a particle energy of 1 GeV) 10⁴⁷–10⁴⁸ erg, under the simple one-zone model. This energy is comparable to the energy content in leptons estimated, for example, in the SNR RCW 86 (H.E.S.S. Collaboration et al. 2018d).

The gamma rays could also be caused by hadronic interactions between energetic cosmic rays, accelerated in an SNR, and ambient matter. The fit to the data yielded a particle spectral index of ${1.69}_{-0.9}^{+0.06}$ and a cutoff energy of ${58}_{-16}^{+14}$ TeV. The required total energy in the hadrons, in terms of the distance to the source and the number density of the target material (n) is (2.2 ± 0.2) × 10⁴⁹ (d/1 kpc)² (1 cm⁻³/n) erg. For a range of plausible distances and typical gas densities the required total energy in the particles could be consistent with theory and observations (∼0.03–0.3 E_SN, e.g., Drury et al. 1994; Ackermann et al. 2013). However, the particle distribution derived from the data is harder than that predicted in the test-particle diffusive shock acceleration theory (see, e.g., Bell 1978), and considerably harder than inferred for typical hadronic SNRs (e.g., Ackermann et al. 2013). Thus the leptonic SNR scenario might be a more plausible explanation for the origin of the emission.

Deeper observations with HAWC and more detailed studies looking for counterparts at other wavelengths are necessary to understand the nature of the source 2HWC J2006+341. The work described here is a first step in this direction. For example, observations of this source in the X-rays could help constrain the maximum electron energies while lower energy gamma-ray observations could probe for the characteristic pion bump of the hadronic scenario. Follow-up observations of the GeV point source PS J2005.0+3419 could also reveal its nature and any possible association to the GeV–TeV emission. Radio observations could find a previously undetected SNR. The matching morphologies at GeV and TeV energies make 2HWC J2006+341 consistent with an SNR shell with no clear counterparts at lower energies, perhaps similar to the TeV sources HESS J1614-518 and HESS J1912+101, or to the extended GeV source G350.6-4.7 (with a spectral index of ∼1.7 and no known counterpart at lower energies; Araya 2018) and other similar objects (Ackermann et al. 2018), having a similar spectrum at GeV energies. Consistent models for a possible population of such objects have yet to be explored.

We thank the anonymous referee for helpful suggestions that improved this work. We acknowledge the support from the US National Science Foundation (NSF); the US Department of Energy Office of High-Energy Physics; the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory; Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051, 232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288, A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grants IG101320, IN111315, IN111716-3, IN111419, IA102019, IN112218; VIEP-BUAP; PIFI 2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at Los Alamos National Laboratory; Polish Science Centre grant, DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la Universidad Michoacana; Royal Society - Newton Advanced Fellowship 180385; Generalitat Valenciana, grant CIDEGENT/2018/034; Chulalongkorn University's CUniverse (CUAASC) grant; the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 690575; Universidad de Costa Rica grants B9171, B6509, B5198. Thanks to Scott Delay, Luciano Díaz, and Eduardo Murrieta for technical support.

Facilities: HAWC - , Fermi. -

Software: astropy (Astropy Collaboration et al. 2013), 3ML (Vianello et al. 2015), naima (Zabalza 2015), fermitools.