HAWC Study of the Ultra-high-energy Spectrum of MGRO J1908+06

MGRO J1908+06, which is in the Galactic plane, was originally discovered by the Milagro detector (Abdo et al. 2007) in the very-high-energy regime (median energy ∼20 TeV) and subsequently confirmed by other TeV observatories, including H.E.S.S, VERITAS, and ARGO (Aharonian et al. 2009; Aliu et al. 2014; Bartoli et al. 2012). Both imaging atmospheric Cherenkov telescopes (IACTs) report that the source is extended and has a fairly hard spectral index (≤2.2), while ARGO reports a slightly softer index (2.54 ± 0.36).

Emission from this region has recently been detected by the High Altitude Water Cherenkov (HAWC) Observatory and included in the collaboration's third catalog (3HWC) with the source name 3HWC J1908+063 (Albert et al. 2020). This emission is centered at (28705, 639) in the (R.A., decl.) J2000 coordinate system. The 1σ statistical uncertainty in this location is 006.

Recently, both HAWC and LHAASO reported that this region is one of only a handful emitting above 100 TeV (Abeysekara et al. 2020; Cao et al. 2021). Even at the highest energies, the source remains extended, with the HAWC collaboration reporting a Gaussian width of 052 above 56 TeV.

The high-energy emission makes this an intriguing source to study. The cosmic-ray spectrum contains a bump known as the "knee" around 1 PeV (Particle Data Group et al. 2020). Which Galactic sources are capable of accelerating particles to this energy is still an open question. Cosmic-ray interactions with their environment produce neutron pions, which then, in turn, decay to gamma rays. These gamma rays are approximately one order of magnitude less in energy compared to the primary cosmic ray. Therefore, in order to probe the knee of the cosmic-ray spectrum, studies of 100 TeV gamma rays are essential. These 100 TeV gamma rays can also be produced in other ways, such as inverse Compton scattering. In fact, leptonic mechanisms likely dominate the highest-energy sky (Albert et al. 2021; Breuhaus et al. 2021a, 2021b; Sudoh et al. 2021)

There are several objects in the region that could serve as counterparts to the TeV emission, including a supernova remnant (SNR) and several pulsars. The radio SNR G40.5-0.5 (Green 2009) has an estimated age of 20–40 kyr (Downes et al. 1980). The distance to this supernova remnant is quite uncertain, with distances between 3.4 kpc and 8.7 kpc appearing in the literature (Duvidovich et al. 2020). The SNR is offset from the TeV emission seen by VERITAS, which led the collaboration to conclude that it is not the main source of the TeV emission (Aliu et al. 2014). The SNR has an angular separation of 029 from the center of the 3HWC source.

There are also three pulsars nearby. The most energetic is PSR J1907+0602, a radio-faint pulsar that was discovered in the GeV energy range by the Fermi Large Area Telescope (Fermi-LAT) (Abdo et al. 2010). According to the ATNF pulsar database ³⁷ (Manchester et al. 2005), this pulsar has a high spin-down power ($\dot{E}$) of 2.8 × 10³⁶ erg s⁻¹, a characteristic age of 19.5 kyr, and is estimated to be 2.37 kpc from the Earth.

While it is possible that this pulsar could have been born in SNR G40.5-0.5, it seems somewhat unlikely. Given the present distance between them (∼28 pc, assuming both objects are the same distance from the Earth), a transverse velocity three times higher than typical due to a supernova kick is required (Aliu et al. 2014).

There has been X-ray emission observed near PSR J1907+0602 with the Chandra X-ray Observatory (Abdo et al. 2010). Originally, it was thought that this emission may be slightly extended and could be the X-ray pulsar wind nebula (PWN), but recently Li et al. (2021) have shown that the X-ray data is consistent with a point source and can be attributed to the pulsar. Additionally, a study using XMM-Newton data found no diffuse or extended emission in a 0.75 by 075 area (Pandel 2015) around the TeV source.

A second pulsar in the region, PSR J1907+0631, is not as energetic as PSR J1907+0602 ($\dot{E}$ of 5.3 × 10³⁵ erg s⁻¹). ³⁸ Its location is very close to the center of the SNR G40.5-0.5 and its estimated age is consistent with the estimated age of the SNR. It has been suggested that this is the pulsar formed in the supernova event (Lyne et al. 2017).

The centroid of the HAWC emission lies roughly in between these two pulsars, roughly ∼03 from both of them. Figure 1 shows the HAWC significance map of the region with the possible counterparts to the emission labeled.

Zoom In Zoom Out Reset image size

There is also a third pulsar in the region, PSR J1905+0600, but it is weaker ($\dot{E}$ = 5.1 × 10³² erg s⁻¹) and older (6 Myr) than the other two pulsars (Hobbs et al. 2004), so it is not expected to contribute to the observed TeV emission. This pulsar is ∼087 away from the centroid of the HAWC emission. Its distance is 8.8 kpc away from Earth.

If any of the emission is hadronic in nature, neutrinos are also expected. Due to the size and hard spectrum of the TeV gamma-ray emission, MGRO J1908+06 has been frequently suggested as a target for neutrino searches (Gonzalez-Garcia et al. 2009; Halzen et al. 2017). In an IceCube search for point-like sources in the astrophysical muon neutrino flux, this region had the second-best p value (highest for a Galactic source), although still consistent with background (Aartsen et al. 2019).

A radio study of the region (Duvidovich et al. 2020) resulted in the discovery of new molecular clouds to the eastern, southern, and western borders of the radio shell of the SNR. The authors of that paper hypothesized that the TeV emission actually consists of two separate components: a leptonic component powered by PSR J1907+0602 and a hadronic component produced by interactions between G40.5-0.5 and the newly discovered molecular gas. If SNR G40.5-0.5 is 8.5 kpc away, as some distance estimates suggest, it is much further away than PSR J1907+0602 and the source we see may actually consist of two superimposed sources. Crestan et al. (2021) also suggest that the emission is comprised of two populations.

Recent observations using Fermi-LAT (Li et al. 2021) have resulted in the detection of extended GeV gamma-ray emission in this area, said to be the GeV counterpart of the TeV emission. This emission contains two components: a soft, low-energy (<10 GeV) component and a harder (>10 GeV) component. The first component is attributed to molecular clouds surrounding the supernova remnant, while the second is likely leptonic in origin and originates from the PWN of PSR J1908+0602.

In this work, we use data from the HAWC Observatory to study 3HWC J1908+063. The HAWC detector consists of 300 water Cherenkov detectors, each instrumented with four photomultiplier tubes. It is designed to detect the byproducts of the extensive air showers that are induced when a gamma ray or a cosmic ray enters the Earth's atmopshere and interacts with particles there.

Located in the state of Puebla, Mexico, HAWC is sensitive to sources with declinations between −26° and +64°. It is capable of continuously monitoring the sky and has achieved a sensitivity of a few percent of the Crab flux over the last five years (Albert et al. 2020). More information on the design of HAWC can be found in Smith (2015) and Abeysekara et al. (2017a).

This paper uses a data set consisting of 1343 days of data collected between 2015 June and 2019 June. The data is binned using a 2D scheme of the estimated energy ($\hat{E}$) and the fraction of the HAWC array hit during an air-shower event, as described in Abeysekara et al. (2019). The estimated energy bins are each a quarter decade in width in log₁₀ space; the first bin starts at $\hat{E}$ = 1 TeV and the last bin ends at $\hat{E}$ = 316 TeV. The "ground parameter" energy estimator is used. This algorithm uses the fit to the lateral distribution function to measure the charge density 40 m from the shower core, along with the zenith angle of the air shower, to estimate the energy of the primary gamma ray. The standard quality cuts described in Abeysekara et al. (2019) are used.

The paper is organized as follows. In Section 2, we describe the diffusion model we use to fit data in the 3HWC J1908+063 region. Section 3 gives the best-fit results using this diffusion model. We also compare the results presented here to those obtained by other observatories. A potential spectral hardening feature at the highest energies is also discussed. In Section 4, we discuss possible models to describe the TeV emission from HAWC. In Section 5, we discuss implications of this model for detection by observatories operating at different wavelengths and with different messengers. In Section 6, we present the conclusions.

Figure 2 shows the HAWC significance map and the best-fit model for this region. Figure 3 shows the residual map (the significance map of the difference between the data and the model). The best-fit parameters can be seen in Table 1, and the HAWC spectral energy distribution (SED) can be seen in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Best-fit Values for the Continuous Injection Diffusion Model

Parameter	Best-fit Value	Statistical Uncertainty	Systematic Uncertainty
θd	178	±008	${}_{-0.28}^{+0.07}$
ϕ0	1.17 × 10−13 (TeV cm2 s)−1	±0.06 × 10−13 (TeV cm2 s)−1	${}_{-0.23}^{+0.10}$×10−13 (TeV cm2 s)−1
α	2.545	±0.026	${}_{-0.06}^{+0.01}$
β	0.134	±0.018	${}_{-0.03}^{+0.02}$
ϕSS433E	2.0 × 10−16 (TeV cm2 s)−1	${}_{-0.7}^{+1.0}$×10−16 (TeV cm2 s)−1	${}_{-0.1}^{+0.2}$×10−16 (TeV cm2 s)−1
ϕSS433W	3.0 × 10−16 (TeV cm2 s)−1	${}_{-0.8}^{+1.1}$×10−16 (TeV cm2 s)−1	${}_{-0.6}^{+0.2}$×10−16 (TeV cm2 s)−1

Note. The first four variables pertain to 3HWC J1908+063 while ϕ_SS433E and ϕ_SS433W are the ϕ₀ values for the east and west lobes of SS433, respectively (see Equation (1)). The pareameter θ_d is reported at the gamma-ray pivot energy of 10 TeV, which corresponds to electrons with energy ∼65 TeV. The column labeled "Systematic Uncertainty" contains the uncertainty from mis-modeling of the detector along with the uncertainty related to modeling of the Galactic diffuse emission. These two uncertainties are combined in quadrature. See Section 3.2 for a discussion of systematic uncertainties.

Download table as:  ASCIITypeset image

A test statistic (TS) for each source in the model is computed. TS is defined as twice the likelihood ratio:

$$\begin{eqnarray}&&{TS}=2\mathrm{ln}\left(\displaystyle \frac{{L}_{s}}{{L}_{n}}\right),\end{eqnarray} \tag{ 8 }$$

where L_s is the best-fit likelihood for the hypothesis including the source and L_n denotes the null hypothesis. When Wilks' theorem (Wilks 1938) is assumed, the square root of the TS can be equated to Gaussian significance. Wilks' theorem is valid for HAWC data (Abeysekara et al. 2017c).

The source 3HWC J1908+063 is significantly detected (>2σ) in each energy bin. Using the method presented in Abeysekara et al. (2019), this corresponds to a detection between 470 GeV and 213 TeV. The TS for this source is 1924.9. The east and west lobes of SS433 have TS values of 8.5 and 15.0, respectively. The TS values for the lobes of SS433 are lower than previously published because the analysis presented here uses an additional quality cut that removes all events whose reconstructed shower core is off the main HAWC array. This greatly increases the angular resolution of the detector, allowing for better morphological studies, but removes approximately half of the data.

The spectral index of 3HWC J1908+063, α, is 2.545 ± 0.026 ${}_{-0.04}^{+0.01}$, and the curvature parameter, β, is 0.134 ± 0.018 ${}_{-0.03}^{+0.01}$. The first set of uncertainties are statistical, and the second set is systematic uncertainties stemming from the mis-modeling of the detector (discussed further in Section 3.2). The flux normalization, ϕ₀ is (1.17 ± 0.06 ±0.10) × 10⁻¹³ (TeV cm² s)⁻¹.

The values reported here are different from those reported in the HAWC highest-energy (>56 TeV) catalog (Abeysekara et al. 2020); most notably the source has a slighter softer spectral index than previously reported. This can be attributed to two differences in the analysis. First, in this paper, we report the spectrum assuming that the reported 3HWC source location is the center of the source (R.A. = 28705 ± 006, decl. = 639 ± 006), while in Abeysekara et al. (2020), the peak of the >56 TeV emission (R.A. = 28691 ± 010, decl. =632 ± 009) was assumed to be the center of the source. These coordinates are 016 away from the 3HWC source location. Second, the high-energy catalog assumed a Gaussian morphology instead of the diffusion model presented here. These two models predict different fluxes, especially at the lower energy end of the spectrum.

The source 3HWC J1908+063 has a best-fit diffusion angle, θ_d, of 1.78 ± 0.08 ${}_{-0.02}^{+0.07}$ degrees, reported at the gamma-ray pivot energy of 10 TeV. The 10 TeV gamma rays imply ∼65 TeV electrons (Equation (6)). Assuming that the distance to the pulsar J1907+0602, 2.37 kpc, is the same as the distance to the source, this means that the 65 TeV electrons will have diffused ∼74 pc (see Equation (4)).

The cooling time of 65 TeV electrons is ∼14,000 yr assuming a magnetic field of 3 μG. This is slightly smaller than the age of the source, so it is used in Equation (5) to calculate the diffusion coefficient. The diffusion coefficient, D, for 65 TeV electrons is 2.92 × 10²⁸ cm² s⁻¹. While this is lower than the interstellar medium (ISM) value by approximately one order of magnitude, this value is not unrealistic. Lower diffusion coefficients have been observed before (Abeysekara et al. 2017b) and a suppression of the diffusion coefficient has been predicted in certain cases, such as near TeV halos (Evoli et al. 2018).

Figure 5 shows the radial profile of 3HWC J1908+063. As can be seen from the figure, the data matches the expectation for continuous injection of particles with diffusion away from the center of the source. Emission is seen out to several parsecs.

Zoom In Zoom Out Reset image size

The flux normalization, ϕ₀, for the East (West) lobe of SS433 is 2.0${}_{-0.7}^{+1.0}$ (3.0${}_{-0.8}^{+1.1}$) × 10⁻¹⁶ (TeV cm² s)⁻¹. These values agree with those reported in Abeysekara et al. (2018), within uncertainties.

A comparison between the spectrum reported here and results from other detectors is discussed in Section 3.3.

Systematic uncertainties related to the modeling of the HAWC detector are investigated as described in Abeysekara et al. (2019). The effects of the detector systematic uncertainties are shown as the pink band in Figure 4. These uncertainties are a function of energy. The width of the systematic uncertainties on the flux range from ${\sim }_{-11 \% }^{+14 \% }$ at 1 TeV to ${\sim }_{-10 \% }^{+22 \% }$ at 100 TeV.

The other main source of uncertainty is the effect of Galactic diffuse emission (GDE) on the spectrum. The base model neglects diffuse emission, as was done with past analyses of the region. The source is located 08 away from the Galactic plane but is quite extended, with some emission observed at b = 0. Here, we investigate whether diffuse emission has any effect on the reported fits.

The GDE is modeled as a Gaussian distribution in Galactic latitude, centered on the Galactic plane. The width of the Gaussian distribution, σ, is allowed to float in the fit. The emission is assumed to emit according to a power-law spectrum:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={\phi }_{0}{\left(\displaystyle \frac{E}{7\mathrm{TeV}}\right)}^{-\alpha }.\end{eqnarray} \tag{ 9 }$$

The value of α is fixed at 2.75, chosen to match the observed cosmic-ray spectrum around 10 TeV (Abdo et al. 2008; Zhou et al. 2017). The flux normalization, ϕ₀, is a free parameter in the fit, along with σ.

The best-fit values for the GDE are ϕ₀ = (1.9 ± 0.4) × 10⁻¹⁴ (TeV cm⁻² s)⁻¹ and σ = 064 ± 017.

The best-fit values for 3HWC J1908+063, with diffuse emission included, are θ_d = 150 ± 010; ϕ₀ = (0.96 ±0.08) × 10⁻¹³ (TeV cm⁻² s)⁻¹; α = 2.505 ± 0.032; and β =0.150 ± 0.022.

The estimated systematic uncertainty of the diffusion angle due to diffuse emission is −15.7%. The effect on the flux varies by energy, ranging from 31.2% at 1 TeV to 17.7% at 100 TeV. The effect is larger at lower energies.

A follow-up analysis investigating the morphology of 3HWC J1908+063 is in progress. This will include an in-depth study of the energy-dependent effect of the diffuse emission.

Figure 6 compares the HAWC result (black line) to those obtained using other detectors, including both IACTs and other air-shower arrays. HAWC measures a higher flux than IACT detectors. Differences in the field of view, angular resolution, and background estimation methods between IACTs and HAWC contribute to discrepancies in the measured flux (Abdalla et al. 2021).

Zoom In Zoom Out Reset image size

Additionally, as discussed previously in Albert et al. (2021), IACTs may extract their spectrum from a region that may be different than the measured morphology. This is a difference from HAWC, where the spectrum and morphology are fit simultaneously.

These effects combined lead to a systematic flux offset between HAWC and IACTs. To account for this, the VERITAS points (Aliu et al. 2014) have been scaled by a factor of 2.03 using the technique outlined in Brisbois (2019). The H.E.S.S. points, taken from the H.E.S.S. Galactic Plane survey (Abdalla et al. 2018), are not scaled as the H.E.S.S. collaboration took this effect into account.

It is apparent that, even with the scaling described above, the HAWC result includes more flux than the IACTs. The discrepancy is more prominent at the lower energies.

There are also differences between HAWC and other air-shower arrays. Note the Milagro result (magenta point). The ARGO spectrum (Bartoli et al. 2012) agrees with the HAWC spectrum at the very lowest energies but quickly begins to diverge as the energy increases. It is interesting to note that the LHAASO result (Cao et al. 2021) agrees with the HAWC result relatively well, especially at the highest energies.

The source of this difference can be attributed to differences between the morphologies assumed by the other observatories. HAWC is unique in using a diffusion model, which has a long tail. If the HAWC data is instead fit assuming a Gaussian morphology with the extent fixed to the IACT extent (here, 044 is used to match the VERITAS extent), the discrepancy disappears (see the blue curve in Figure 6). In this study, the R.A. and decl. are left fixed to the 3HWC coordinates. Due to their background estimation methods and morphological assumptions, IACTs are not sensitive to the long, extended tails of the gamma-ray emission. This may affect their physics conclusions.

The HAWC Gaussian flux points are consistent with the HAWC diffusion flux points above 100 TeV. This will be important in the Theoretical Modeling section (Section 4) where these data points are used to draw conclusions.

We look for evidence that the spectrum hardens at the highest energies. Such a discovery could indicate that there are hadronic contributions; if the source were entirely leptonic, curvature in the spectrum is expected at the highest energies due to Klein-Nishina effects (Moderski et al. 2005). This study is motivated by the last two flux points shown in Figure 4, which appear as if they are beginning to deviate from the best-fit log-parabola shape by exhibiting a slight flux enhancement over expectation. It is advantageous to have more flux points at the highest energies because this is where the deviations between leptonic and hadronic models will have the most prominent differences. To investigate this, the three highest-energy bins, which each span a quarter decade in log_10-energy space ($\hat{E}\gt 56$ TeV) are subdivided into six smaller bins, each spanning an eighth of a decade in width, before rerunning the spectral fit.

Figure 7 shows the spectrum of 3HWC J1908+063 with the highest-energy bins split into sub-bins. As in the nominal analysis, a log-parabola spectral shape and the diffusion morphology are used. The overall forward-folded fit results are compatible with those presented in Section 3.1: θ_d = 1.77 ±0.08, ϕ₀ = 1.16 × 10⁻¹³, α = −2.544 ± 0.026, and β = 0.133 ±0.018. All uncertainties are statistical only.

Zoom In Zoom Out Reset image size

A potential spectral hardening feature is visible in the flux points by eye, beginning around ∼70 TeV. The source is not significantly detected in the last energy bin ($234\lt \hat{E}\lt 316$ TeV; TS = 0.39) so a 90% confidence level upper limit is computed. In all of the other sub-bins, the TS is >16 (>4σ).

This spectral hardening feature is presently not significant. The last two flux points in Figure 7 that are significantly detected are only 1.5σ and 1.8σ away from the best-fit log-parabola line. Adding these two points in quadrature, the total significance of this potential feature is 2.3σ. Improvements to the HAWC reconstruction algorithms or studies using future detectors will decrease the uncertainties and make a more definitive statement. It is worth noting that this feature is remains visible by eye if detector response functions designed to probe systematic uncertainties related to the mis-modeling of the HAWC detector are used, so it is unlikely to be an instrumental effect

If this feature is shown to be significant in the future, it would be strong evidence for hadronic processes associated with this source. See the Theoretical Modeling section (Section 4) for a further discussion of this.

If the >100 TeV gamma rays originate from hadronic processes, this requires proton energies corresponding to the knee of the cosmic-ray spectrum. It is unclear what the origin of these cosmic rays would be. One scenario is that they could be from the SNR interacting with molecular clouds in the region. Another scenario is that hadronic acceleration mechanisms occur in PWN. While this does not occur in conventional models, this possibility has been explored in the literature (Bednarek 2003; Di Palma et al. 2017).

A purely hadronic model is disfavored due to a lack of sufficient energy to power the hadronic emission. To show this, first we calculate the target proton density in the region. Molecular clouds in the region can affect this density. There are several known clouds near 3HWC J1908+063. Here, we use data from a CO survey published in Dame et al. (2001). Data from the 35–50 km s⁻¹ range is used, as this corresponds to the distance to PSR J1907+0602. Atomic gas is also considered, using the HI4PI survey (HI4PI Collaboration et al. 2016). Assuming a circle of radius 1° centered on the HAWC source, the average combined column density from these two surveys totals 9.9 × 10²¹ cm⁻². A rough estimate of the mass can then be obtained. The radius of the 1° circle is 41.36 pc, assuming that the gas is 2.37 kpc away (the distance to PSR J1907+0602). Then,

$$\begin{eqnarray}&&m=\pi {r}^{2}{{Nm}}_{{\rm{p}}},\end{eqnarray} \tag{ 10 }$$

where m is the total mass in the region, r is the radius, N is the column density, and m_p is the mass of the proton. The total mass is 8.5 × 10³⁵ kg, which leads us to a proton density (n) of ∼60 protons cm⁻³. This is considerably higher than the typical ISM value of 1 cm⁻³. This calculation assumes that the gas is spherical.

Assuming that the source is purely hadronic and that the population of protons is trapped within the cloud, the total energy in nonthermal hadrons is

$$\begin{eqnarray}&&{W}_{\mathrm{pp}}={L}_{\gamma }\frac{{t}_{\mathrm{pp}}}{{\eta }_{{\rm{n}}}}\end{eqnarray} \tag{ 11 }$$

where L_γ is the gamma-ray luminosity:

$$\begin{eqnarray}&&{L}_{\gamma }=4\pi {d}^{2}{\int }_{{E}_{\min }}^{{E}_{\max }}\displaystyle \frac{{dN}}{{dE}}{EdE},\end{eqnarray} \tag{ 12 }$$

t_pp, the proton energy loss timescale, is equal to 2 × 10¹⁵ n⁻¹ sec cm⁻³ and η_n is equal to 1.5. The parameter η_n accounts for the fact that there are nuclei heavier than hydrogen in both cosmic rays and interstellar matter (HESS Collaboration et al. 2016). In Equation (12), $\tfrac{{dN}}{{dE}}$ is the gamma-ray flux per energy. It is not related to N (the column density) from Equation (10).

From the fitted HAWC spectrum (see Section 3.1), E_min is 470 GeV and E_max is 213 TeV. Using the best-fit values for the spectrum to calculate L_γ and 60 protons cm⁻³ for the number density, n, we obtain a value of 1.7 × 10⁴⁸ erg for W_pp over the energy range of HAWC. This is the total energy available based on the gas in the region. However, in order to explain the HAWC data points, a simple hadronic model requires a total energy on the order of ∼10⁵⁰ erg, assuming that n is ∼60 cm⁻³ and the parent proton spectrum can be modeled by a single power law. This is two orders of magnitude higher in energy than the available energy calculated above, so the hadronic model is difficult to explain.

We cannot completely reject a one-population model because the lower energies are unconstrained by the HAWC data. The energy requirement of 10⁵⁰ erg assumes that the protons can be modeled by a power law all the way down to the lowest energies. However, if the spectrum of the protons is instead curved, it is conceivable that the total energy required could be lower, as calculating the energy budget involves integrating over the spectrum.

An inverse Compton (IC) origin of the emission is attractive given the absence of a clear correlation with target material.

In a single zone/population model and in the case that synchrotron losses are dominant at all energies, the equilibrium high-energy electron spectrum is steep and Klein-Nishina (KN) suppression results in a very steep IC spectrum that is not consistent with the HAWC data. Alternatives to this picture are the existence of more than one emission component (see Section 4.3) or that the magnetic field energy density lies below that of radiation fields such that IC cooling effects become important. This case was considered in detail by Breuhaus et al. (2021a, 2021b) and here we adopt the same approach.

Electrons following an exponential cutoff power-law equation

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}{\left(\displaystyle \frac{E}{1\mathrm{erg}}\right)}^{-\alpha }\exp \left(-\displaystyle \frac{E}{{E}_{\mathrm{cut}}}\right)\end{eqnarray} \tag{ 13 }$$

are injected into a region with constant, isotropic, and homogeneous radiation fields and a fixed magnetic field of 3 μG. The radiation field adopted is that of the large-scale Galactic radiation field model of Popescu et al. (2017) at the location of the possibly associated pulsar PSR 1907+0602 (plus the CMB). The injection at a constant rate takes place for 19.5 kyr, the characteristic age of the pulsar. A fit to the data with N₀, the injection index α and the cutoff energy E_cut as free parameters gives N₀ = (1.3 ± 0.2) × 10³⁶ erg s⁻¹ and α = 2.68 ± 0.04. The cutoff energy E_cut can only be constrained to be larger than 610 TeV with 95% confidence and therefore a value of E_cut = 10 PeV is used. The quoted values correspond to a mean electron/positron luminosity of 50% of the current spin-down power of the pulsar injected above 1 TeV. For the numerical calculations, we made use of the open-source code GAMERA (Hahn 2015).

To account for absorption of the γ-ray emission due to pair production in the interstellar photon fields on the way to Earth, the same radiation model together with the CMB was used. The model is shown in Figure 8.

Zoom In Zoom Out Reset image size

Assume instead that there is a region near the PWN where the magnetic field is high (>10 μG) and there is a TeV halo driven by the pulsar, so that u_B is very high. Electrons have been continuously injected with constant power over the lifetime of the system. This leads to a complete dominance of synchrotron losses over Inverse Compton losses. Under these conditions, one population of leptonic particles cannot explain the TeV emission. The inverse Compton component suffers from considerable Klein-Nishina effects at the highest energies, which suppresses the high-energy gamma-ray flux.

In this section, we explore the possibility that there are two populations responsible for the TeV particles present in the emission region: a primary population of nonthermal electrons that is responsible for the bulk of the emission, and a secondary nonthermal particle population that dominates the emission above 50 TeV. The secondary component may be either of leptonic or hadronic origin.

Due to the unconstrained parameter space, the number of free parameters in the models described in this section are larger than the number of data points. Therefore, both models presented here are merely possible combinations of parameters that describe the data well. Both models assume that the source is located 2.37 kpc from the Earth and the ISM magnetic field strength is assumed to be 3 μG. Only contributions from the CMB photon field are considered. The spin-down power of the pulsar is assumed to be 2.8 × 10³⁶ erg s⁻¹.

In both two-population models presented here, the first component is a nonthermal electron distribution with a broken-power-law shape. The electron spectrum is given by:

$$\begin{eqnarray}\frac{{{dn}}_{{\rm{e}}}}{d{\gamma }_{{\rm{e}}}}=\left\{\begin{array}{cc}{n}_{{\rm{e}}0}{\left(\frac{{\gamma }_{{\rm{e}}}}{{\gamma }_{{\rm{e}}0,b}}\right)}^{-{p}_{{\rm{e}}0}}{e}^{-{\gamma }_{{\rm{e}}}/{\gamma }_{{\rm{e}}0,\mathrm{cut}}} & {\gamma }_{{\rm{e}}}\lt {\gamma }_{{\rm{e}}0,b}\\ {n}_{e0}{\left(\frac{{\gamma }_{{\rm{e}}}}{{\gamma }_{{\rm{e}}0,b}}\right)}^{-({p}_{{\rm{e}}0}+1)}{e}^{-{\gamma }_{{\rm{e}}}/{\gamma }_{{\rm{e}}0,\mathrm{cut}}} & {\gamma }_{{\rm{e}}}\gt {\gamma }_{{\rm{e}}0,b}\end{array}\right.,\end{eqnarray} \tag{ 14 }$$

where γ_e is the electron Lorentz factor, n_e0 is the spectral normalization factor, γ_e0,b is the Lorentz factor at the cooling break, γ_e0,cut is the electron Lorentz factor at the spectral cutoff, and p_e0 is the spectral index. The spectral index changes by one unit after the break, as is expected from radiative cooling. The spectrum begins at ${\gamma }_{{\rm{e}}0,\min }$ = 10³ (∼500 Mev).

The nonthermal electron energy density in the emission region is given by

$$\begin{eqnarray}&&{u}_{{\rm{e}},n\mathrm{th}}={\int }_{1}^{\infty }\frac{{{dn}}_{{\rm{e}}}}{d{\gamma }_{{\rm{e}}}}{\gamma }_{{\rm{e}}}{m}_{{\rm{e}}}{c}^{2}d{\gamma }_{{\rm{e}}}.\end{eqnarray} \tag{ 15 }$$

Values of γ_e0,b = 5 × 10⁷ (26 TeV), γ_e0,cut = 1.2 × 10⁸ (61 TeV), and p_e0 = 2 provide a good description of the data. Assuming that the emission region is homogeneous, the total nonthermal electron energy is 1.6 × 10⁴⁸ erg. The location of the cooling break (γ_e0,b ) corresponds to a cooling time of τ_cool = 2.57 × 10⁴ yr, which is roughly the same as the characteristic age of the pulsar. 93% of the spin-down power of the pulsar is contained in this component.

4.3.1. Two-population Purely Leptonic Model

To describe the hard spectrum at ≳50 TeV, we introduce a second nonthermal electron population. The second population is assumed to be from a more recent active phase of the source and has the form:

$$\begin{eqnarray}&&\frac{{{dn}}_{{\rm{e}}}}{d{\gamma }_{{\rm{e}}}}={n}_{{\rm{e}}1}{\gamma }_{{\rm{e}}}^{-{p}_{{\rm{e}}1}}{e}^{-{\gamma }_{{\rm{e}}}/{\gamma }_{{\rm{e}}1,\mathrm{cut}}},\end{eqnarray} \tag{ 16 }$$

where the spectrum extends from ${\gamma }_{{\rm{e}}1,\min }={10}^{3}$ (∼500 MeV) to γ_e1,cut = 10⁹ (∼500 TeV) and p_e1 = 2. The total energy in this component is 4 × 10⁴⁶ erg. This is only a fraction of the primary leptonic component. We assume that the electrons in this second component were injected in the last ∼2000 yr. This component contains 2% of the spin-down power of the pulsar.

Note that the second pulsar in the region, PSR J1907+0631, is unlikely to be able to accelerate particles to the energies discussed here, so this second component is still likely associated with PSR J1907+0602.

This model is shown as the thick blue line in Figure 9. Note that this model violates the X-ray upper limit from XMM-Newton. However, we note that the XMM-Newton upper limit is extracted from a region that is smaller (a 075 by 075 square) than the HAWC extent. The observation was centered on the pulsar, which is ∼03 away from the centroid of the HAWC source. It is difficult for X-ray satellites to observe regions that are 1 degree across, as this source is, due to their relatively small fields of view. Nevertheless, the second electron population results in extra flux in the X-ray band, which may be examined by future X-ray telescopes.

Zoom In Zoom Out Reset image size

4.3.2. Two-population Hybrid Lepto-hadronic Model

A hybrid lepto-hadronic model can also explain the HAWC data well. In this scenario, the extra nonthermal particle distribution consists of protons. The TeV gamma-ray emission is created via proton-proton collision. The hadronic population is modeled as follows:

$$\begin{eqnarray}&&\frac{{{dn}}_{{\rm{p}}}}{d{\gamma }_{{\rm{p}}}}={n}_{{\rm{p}}}{\gamma }_{{\rm{p}}}^{-{p}_{p}}{e}^{-{\gamma }_{{\rm{p}}}/{\gamma }_{{\rm{p}},\mathrm{cut}}}\end{eqnarray} \tag{ 17 }$$

where the proton spectrum extends from ${\gamma }_{{\rm{p}},\min }=1$ (∼900 MeV) to γ_p,cut = 10⁷ (∼9 PeV), and the spectral index p_p is 2. Assuming these model parameters, the energy in protons is 1.5 × 10⁴⁸ ergs. Note that this is close to the energy in the electrons. The solid red line in Figure 9 shows the predicted flux for the lepto-hadronic model. There is much less tension with the X-ray upper limit than in the two-population leptonic model. This is expected, because the pair synchrotron emission, a byproduct of proton-proton collision, gives trivial contribution to the emission due to the very low ISM magnetic field (3 μG). Note that this model also overshoots the X-ray upper limit, although only slightly.

A major issue is how to accelerate nonthermal protons in the source. There are two possible scenarios. The first scenario is that the protons were accelerated when the supernova exploded, and then diffused out of the shock. The other scenario is that the protons were accelerated at the termination shock by the pulsar wind. In either scenario, the nearby molecular clouds provide a target for the gamma-ray emission.

4.3.3. Comparison of Two-population Model to Recent Fermi-LAT Results

4.3.4. Comparing the different Two-population Models