© ESO 2020

1. Introduction

Theoretical models for the origin of cosmic rays (CRs) consider young supernova remnants (SNRs) as Pevatrons, that is, as hadron accelerators for energies up to the “knee” of the cosmic-ray spectrum measured at around 3 PeV (Baade & Zwicky 1934; Hillas 2005). The detection of a neutral pion decay signature in γ-ray data from SNR observations confirmed that SNRs accelerate hadronic particles (see e.g., Ackermann et al. 2013; Jogler & Funk 2016). However, the maximal energy to which SNRs can accelerate hadrons is not determined and, in particular, whether or not SNRs are sources of PeV CRs is unclear.

The association of a γ-ray source with a Pevatron can be tested with a measurement of the γ-ray spectrum above TeV energies. An exponential cutoff E_cut,  γ in the γ-ray spectrum at a few hundred TeV is expected when CRs are accelerated up to PeV energies and produce γ-rays in interactions with ambient material (Gabici & Aharonian 2018). The Pevatron nature of a γ-ray source is in turn constrained when an exponential cutoff at an energy well below 100 TeV is detected in the observed γ-ray spectrum of the source. For example, the recent detection of an exponential cutoff at E_cut,  γ = 3.5 TeV in the γ-ray spectrum of Cassiopeia A rules out simple Pevatron models for this young Galactic SNR (Ahnen et al. 2017).

The most constraining lower limits on the energy cutoff E_cut,  h of Galactic hadron accelerators are currently derived from the diffuse γ-ray emission in the vicinity of the radio source Sagittarius A* (SgrA*) and for the γ-ray source HESS J1641−463. Assuming that hadronic particles generate the γ-ray emission from these sources, E_cut,  h >  600 TeV and E_cut,  h >  100 TeV are inferred at 90% and 99% confidence level (CL) respectively (Abramowski et al. 2016, 2014). More than 100 isolated sources of TeV γ-rays are currently known. However, no systematic search for Pevatron candidates among the known γ-ray sources has yet been presented.

This paper is organized as follows. Section 2 outlines the general Pevatron search method followed in this work. The later analysis relies in large parts on the public catalog of the High Energy Stereoscopic System (HESS) Galactic plane survey (HGPS). This catalog is introduced in Sect. 3. A method to constrain energy cutoffs in measurements of power-law energy spectra is discussed in Sect. 4. The method is argued to be applicable to the HGPS catalog introduced above. A search for Pevatron candidates in the HGPS catalog is summarized in Sect. 5. For one Pevatron candidate identified in the previous analysis of HGPS data, additional public data are available from the HESS, Very Energetic Radiation Imaging Telescope Array System (VERITAS), and Milagro experiments. A combined analysis of these data is presented in Sect. 6. Finally, Sects. 7 and 8 discuss the results and summarize the conclusions.

2. Pevatron candidate search method

The detection of TeV γ-rays from an isolated region of the sky can indicate cosmic sites where particles are accelerated to energies of greater than 1 TeV. However, the identification of the primary particle type is nontrivial when only TeV γ-ray data are available. Both leptonic and hadronic scenarios are typically conceivable. Electrons are accelerated in the leptonic scenario and generate γ-rays via the inverse Compton up-scattering of low-energy photons. This scenario is not relevant in the following search for hadronic Pevatrons as origins of Galactic CRs. However, the unambiguous exclusion of the leptonic scenario for a given γ-ray source is very challenging. For example, as recently discussed, even the production of γ-rays with an energy exceeding 100 TeV is possible in a leptonic scenario (Amenomori et al. 2019). The nondetection of a γ-ray energy cutoff at energies below few 100 TeV is therefore not considered as a sufficient criterium for the identification of a hadronic Pevatron. Instead of an unambiguous identification of hadronic Pevatrons, the following search aims to find Pevatron candidates, that is, γ-ray sources where the hadronic Pevatron scenario cannot be excluded. The resulting list of Pevatron candidates allows dedicated multi-messenger analyses with reduced statistical trial factors, involving, for example, searches for high-energy neutrinos as fingerprints of hadronic interactions.

Many known Galactic γ-ray sources are identified as pulsar wind nebulae (PWNe) for which the γ-ray emission is modeled in a leptonic scenario (Gaensler & Slane 2006). Another large class of γ-ray sources has yet to be identified, with sources of radiation in other energy ranges. The search for new Pevatron candidates presented here is restricted to these unidentified γ-ray sources in the Galactic disk. The motivation for the restriction to unidentified γ-ray sources is at least twofold. The exclusion of PWNe removes a large class of γ-ray sources for which the leptonic scenario is very likely. Additionally, Gabici & Aharonian (2007) discuss that the time period where SNRs act as Pevatrons might be short, that is, less than 1000 years. Easier than the direct detection of Pevatrons might therefore be the indirect identification via the detection of delayed γ-ray emission from molecular clouds in the vicinity of SNRs. Gabici & Aharonian (2007) discuss this scenario and suggest that unidentified γ-ray sources can be associated with molecular clouds illuminated by nearby SNRs.

Following Kelner et al. (2006), the spectrum of γ-rays generated in the vicinity of a Pevatron in interactions of accelerated CRs with ambient material is modeled by

$$\begin{array}{c}\hfill \varphi (E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}exp(-\sqrt{16{\lambda }_{\mathrm{h}}E}).\end{array}$$(1)

Here, ${\varphi }_{0}exp(-\sqrt{16{\lambda }_{\mathrm{h}}{E}_0})$ is the flux normalization at energy E₀ and the parent CR spectrum has a cutoff energy of E_cut,  h = 1/λ_h (Kelner et al. 2006). Models for diffuse shock acceleration predict a power-law index Γ ≈ 2 (see e.g., Hinton & Hofmann 2009). Following this model, TeV γ-ray sources that are selected as Pevatron candidates must have a spectrum at TeV energies which is compatible with a power-law with index Γ = 2. Additionally, the inferred lower limit on the energy cutoff E_cut,  h in the hadronic model given by Eq. (1) must be at least O(100 TeV), that is, of the same order of magnitude as the current best constraints derived for HESS J1641−463 and the γ-ray emission in the vicinity of SgrA*.

To further constrain a leptonic scenario for the generation of γ-rays, the presence of a break in the TeV γ-ray spectrum where the power-law index changes can be excluded. For example, in a leptonic scenario in which inverse Compton losses dominate over synchrotron cooling, a Klein–Nishina (KN) break is expected in the γ-ray spectrum at an energy E_KN which depends on the target photon field. For cosmic microwave background (CMB) target photons, the break is expected at a few hundred TeV. For infrared target photons or optical starlight, the energy break is expected around 10 TeV and 30 GeV respectively (Hinton & Hofmann 2009). The γ-ray spectrum at energies far away from E_KN can be described by a broken power-law with E_break = E_KN

$$\begin{array}{c}\hfill \varphi (E)=\{\begin{array}{cc}{\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}\hfill & E<E_{\mathrm{break}}\hfill \\ {\varphi }_0{\left(\frac{E_{\mathrm{break}}}{E_0}\right)}^{\mathrm{\Delta }\mathrm{\Gamma }}{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }-\mathrm{\Delta }\mathrm{\Gamma }}\hfill & E>E_{\mathrm{break}}.\hfill \end{array}\end{array}$$(2)

It is expected that ΔΓ = Γ for the index change at the KN break (see e.g., Hinton & Hofmann 2009). The selection of sources with a power-law index around 2 at TeV energies excludes leptonic scenarios where the inverse Compton target photons are of a much higher energy than CMB photons, for example optical starlight, and the KN break must occur below TeV energies. However, leptonic scenarios with CMB target photons cannot typically be constrained with TeV γ-ray data.

The selection of a γ-ray source as a Pevatron candidate in the following search is not sufficient to identify the source with a Pevatron. Furthermore, leptonic emission scenarios cannot be ruled out and the lower limits on E_cut,  h are not sufficient to argue for the acceleration of hadrons to PeV energies. Additionally, the selection of a γ-ray source as a Pevatron is not necessary for the association with a Pevatron. For example, the γ-ray source Cassiopeia A would not pass the selection due to the measured TeV cutoff discussed in Sect. 1. However, Cassiopeia A is still being discussed as a Pevatron candidate in multi-zone models (Zhang & Liu 2019).

Despite the Pevatron candidate selection being neither necessary nor sufficient for the association of a γ-ray source with a Pevatron, selected sources can be a particularly interesting target for further observations. For example, targeted signal searches with instruments sensitive to high-energy neutrinos towards a few selected Pevatron candidates can be performed.

3. HESS Galactic plane survey data and quality selection

The currently largest catalog of Galactic TeV γ-ray sources is based on the HGPS. This dataset, as far as is relevant for the data analysis below, and additional quality selection criteria are discussed in this section.

3.1. HESS Galactic plane survey data

The HGPS comprises data acquired from 2004 to 2013 with the HESS I (Aharonian et al. 2006b) array of Cherenkov telescopes in Namibia. Large parts of the Milky Way plane were observed to detect γ-ray sources in a nominal energy range of 0.2−100 TeV with an angular resolution better than 0.1°. Observations of the Milky Way plane were performed in a scanning mode with individual observation runs of typically 28 min towards different sky directions. In total, 78 γ-ray sources were detected in the plane survey. The majority (47) of the detected sources are unidentified, that is, no firm identification is currently possible with objects detected in other energy ranges.

The HGPS catalog is based on a general-purpose analysis without special adoption for each individual source. Systematic effects of this analysis cannot be ruled out, in particular for sky regions with multiple γ-ray sources.

An adaptive ring background algorithm is used in the analysis of data from the HGPS. This algorithm estimates the number of signal (N_ON) and background (N_OFF) events and enables a calculation of the excess events N_γ = N_ON − αN_OFF. The acceptance normalization factor α is obtained from independent observations of different sky regions void of γ-ray sources. Energy-dependent information on N_γ, N_ON, N_OFF and α is not part of the public HGPS catalog. However, spectral flux data for each γ-ray source detected in the HGPS are typically available in six logarithmic energy bins. The bin centers, E_i, i = 0, …, N − 1, of spectral points are related via E_i = ξⁱE_min with a binning factor ξ where E_min is the center of the lowest energy bin. Statistical errors on the γ-ray flux are given as asymmetric 68% confidence level (CL) intervals and originate from the Poisson counting error on N_ON and N_OFF. Other error sources, for example an error on the acceptance normalization, are not included in the statistical error on the γ-ray flux. A complete description of the HGPS dataset can be found in Abdalla et al. (2018).

3.2. Data selection

A special quality selection is applied in this analysis of the public HGPS data. The energy bin-width for spectral points is required to be much larger than the instrumental energy resolution to avoid a statistical correlation between spectral points. Following Aharonian et al. (2006b), the energy resolution of the HGPS analysis is assumed to be 15% or better. Only sources in the HGPS with a spectral energy binning factor of ξ >  1.5 are selected.

In the HGPS catalog, data on the error of spectral points are given as a lower (σ_low) and an upper (σ_high) error. The average error σ = 1/2 (σ_low + σ_high) is used in the analysis. Spectral points with an error asymmetry η = max(|σ − σ_low|, |σ − σ_high|)/σ >  10% are discarded from the analysis. Also discarded are spectral points ϕ which are detected at low significance level S = ϕ/σ <  1.5. Sources in the HGPS catalog are not analyzed when less than five spectral points remain after the quality selection. A total of 25 out of 78 sources were selected for further analysis. A summary of the selected HGPS data is given in Table 1.

4. Constraints on the energy cutoff of a power-law spectrum

A statistical hypothesis test can be used to search for an exponential cutoff in a power-law spectrum. A confidence interval for the energy cutoff can be constructed via the inversion of the test. However, special care must be taken in the analysis of public data from the HGPS because the number of available spectral points per γ-ray source is typically small. The small sample size questions the application of frequently applied asymptotic results for the distribution of a test statistic when the null hypothesis is true. Additionally, systematic effects on the acceptance normalization cannot be ruled out and the catalog data are based on a general-purpose analysis. Differences between the estimated and the true acceptance normalization due to instrumental effects could, for some observation runs, lead to outliers in the number of excess events, which are not modeled with the published statistical error. Additional systematic effects can originate from the confusion of γ-ray sources in regions with complex morphology.

This section contains a description of statistical methods to constrain energy cutoffs in power-law spectra. The application of the likelihood ratio (LR) test and the F-test for the presence of cutoffs in power-law spectra are discussed in Sect. 4.1. In Sect. 4.2, it is argued that the F-test is more robust against potential violations of the normal error model of the HGPS catalog data than the LR-test. A method to derive a lower limit on the energy cutoff of a power-law γ-ray energy spectrum is discussed in Sect. 4.3.

4.1. Models and hypotheses tests

Consider a power-law model for the γ-ray flux ϕ(E) at energy E with an exponential cutoff at E_cut,  γ = 1/λ_γ. The model can be parameterized by

$$\begin{array}{c}\hfill \varphi (E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}exp(-{({\lambda }_{\gamma }E)}^{\beta }),\end{array}$$(3)

where again Γ is the spectral index and ϕ₀ the flux normalization at energy E₀. The parameter E₀ is assumed to be fixed to 1 TeV in the following. Also fixed is the parameter β which controls the shape of the exponential cutoff. The case β = 1, for example, is frequently discussed in the literature (see e.g., Ahnen et al. 2017). The hadronic model given by Eq. (1) reduces to Eq. (3) with β = 0.5 and λ_γ = 16λ_h. The cutoff parameter λ_γ is constrained to λ_γ >  0 because a γ-ray flux suppression for energies above E_cut,  γ is expected for Pevatrons. The model given by Eq. (3) has three free parameters θ₁ = (ϕ₀,  Γ,  λ_γ). The question remains as to whether the power-law model with exponential cutoff (Eq. (3)) gives a better description of a dataset than a power-law model parameterized by

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{PL}}(E)={\varphi }_0{\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }}.\end{array}$$(4)

This model has two free parameters θ₀ = (ϕ₀,  Γ) while E₀ is again fixed to 1 TeV. To select the best-fitting model, the null hypothesis H₀, that is, the absence of an exponential cutoff (λ_γ = 0), can be tested against the physically constrained alternative hypothesis H₁, λ_γ >  0, by means of a hypothesis test at a given confidence level (CL).

In the following ϕ_i, i = 0, …, N − 1 is used to denote binned γ-ray flux measurements with bin-centers at energies E_i. All N measurements of ϕ_i are assumed to be independent. In practice, the independence of the ϕ_i can be expected when the energy bin-width for the calculation of ϕ_i is much larger than the instrumental energy resolution. Further, σ_i is used to represent the measurement errors corresponding to the γ-ray fluxes ϕ_i. The measurements ϕ_i are, in a first step, assumed to be normally distributed around a hypothetical true value ${\widehat{\varphi }}_i$. A model ϕ(E_i|θ) with parameters θ predicts the true flux ${\widehat{\varphi }}_i$. Two frequently used tests for a power-law hypothesis (Eq. (4)) against a model with exponential cutoff (Eq. (3)) are presented in the following. Subsequently, the robustness of the tests against deviations from the assumed normal error model is discussed.

4.1.1. Likelihood ratio test

Because the flux measurements are assumed to be independent, the likelihood function for the parameters θ of a given model factorizes to $L(\mathit{\theta })=\prod _{i=0}^{N-1}\frac{1}{\sqrt{2\pi {\sigma }_i^2}}exp(-\frac{{({\varphi }_i-\varphi (E_i|\mathit{\theta }))}^2}{2{\sigma }_i^2})$. Let ${\widehat{\mathit{\theta }}}_k={\mathrm{argmax}}_{{\mathrm{H}}_k}L(\mathit{\theta })$ be the maximum likelihood estimates for the model parameters under the hypotheses H_k with k ∈ {0, 1}. Consider the LR-test with the statistic $-2ln\mathrm{\Lambda }=-2ln(L({\widehat{\mathit{\theta }}}_0)/L({\widehat{\mathit{\theta }}}_1))=\mathrm{RSS}({\widehat{\mathit{\theta }}}_0)-\mathrm{RSS}({\widehat{\mathit{\theta }}}_1)$ where

$$\begin{array}{c}\hfill \mathrm{RSS}({\widehat{\mathit{\theta }}}_k)={\displaystyle \sum _{i=0}^{N-1}}{\left(\frac{{\varphi }_i-\varphi (E_i|{\widehat{\mathit{\theta }}}_k)}{{\sigma }_i}\right)}^2\end{array}$$(5)

for the hypotheses H_k. The null hypothesis is discarded at confidence level CL if −2lnΛ >  Λ_crit(CL). The critical value Λ_crit can be calculated using the results from Chernoff (1954) because the two models given by Eqs. (4) and (3) are nested, that is, Eq. (3) reduces to Eq. (4) when λ_γ → 0. When $F_{{\chi }_{\text{1}}^{\text{2}}}^{-1}$ denotes the inverse cumulative density function of a χ²-distributed random variable with one degree of freedom, the critical value is given by ${\mathrm{\Lambda }}_{\text{crit}}=F_{{\chi }_{\text{1}}^{\text{2}}}^{-1}(2\text{CL}-1)$; see Appendix A. The test will have the expected false positive error rate 1 − CL when the normal error model holds and the sample size N is large such that $\varphi (E_i|{\widehat{\mathit{\theta }}}_k)\to {\widehat{\varphi }}_i$ when H₀ is true.

4.1.2. F-test

The F-test is, in the case of N independent measurements and two nested models with two and three parameters, based on the test statistic

$$\begin{array}{c}\hfill F=(N-3)\frac{\mathrm{RSS}({\widehat{\mathit{\theta }}}_0)-\mathrm{RSS}({\widehat{\mathit{\theta }}}_1)}{\mathrm{RSS}({\widehat{\mathit{\theta }}}_1)}·\end{array}$$(6)

The null hypothesis λ_γ = 0 is discarded in favor of the alternative hypothesis λ_γ >  0 at confidence level CL when F >  F_crit(CL). Similar to the LR-test, the critical value is given by

$$\begin{array}{c}\hfill F_{\mathrm{crit}}=F_{1,N-3}^{-1}(2\mathrm{CL}-1),\end{array}$$(7)

where $F_{1,N-3}^{-1}$ is the inverse cumulative density function of a F-distributed random variable with 1 (numerator) and N − 3 (denominator) degrees of freedom. The F-test is exact for small samples N when linear models are being compared. For nonlinear models, such as in the case considered here, the false-positive error rate can only be expected to be asymptotically equal to 1 − CL (Gallant 1975). Like the special LR-test constructed above, the F-test relies on the assumption of a normal error model.

4.2. Robustness of the LR- and F-tests

The upper plot in Fig. 1 shows spectral data points from the HGPS for the γ-ray source HESS J1800−240. A fit of the data to a power-law model (Eq. (4)) results in the best fit parameters θ_s = (ϕ₀,  Γ) = ((4.3 ± 0.1) × 10⁻¹³ TeV⁻¹ cm⁻² s⁻¹, 2.44 ± 0.02). The robustness of the LR- and F-tests constructed above is tested in a Monte Carlo simulation of true null models where random spectra are created based on the best-fit model ϕ(E|θ_s). For each simulated spectrum, all five spectral points are sampled from a Student t-distribution t₅ (location,  scale) with location parameter ϕ(E_i|θ_s), scale parameter σ_i, and 5 degrees of freedom. For comparison, the same simulation is repeated by sampling from a normal distribution with mean ϕ(E_i|θ_s) and standard deviation σ_i. The Student error model predicts more flux outliers than expected in the normal error model on which the constructed LR- and F-tests rely.

Fig. 1. Spectral fits to the data from the HGPS source for which the most (upper panel) and least (lower panel) constraining lower limit on Ecut,  h is derived. Shown in green are the best fits to the data with a power-law spectrum. Red lines are best fits to the data with a γ-ray power-law spectrum with exponential cutoff at the lower limit Ecut,  γ derived from the HGPS data. Yellow lines are the corresponding best fits to the data in a hadronic model with cutoff at the lower limit Ecut,  h. Limits for the lower plot are at 90% CL. In the upper panel, best fits with lower cutoff limit at 90% CL are hardly distinguishable from the best power-law fit. Therefore, 99% CL best fits are shown in this case.

Figure 2 shows the distribution of LR- and F-test p-values for Monte Carlo simulations with the Student (upper panel) and normal (lower panel) error models. The shown p-values refer to tests of a power-law model against a model with exponential cutoff in the γ-ray spectrum (Eq. (3) with β = 1). Similar simulations with unchanged conclusions were performed when a power-law model is tested against a model where the γ-ray emission is modeled in a hadronic scenario (Eq. (1)). The lower panel of Fig. 2 shows that the distribution of p-values for a true null hypothesis is compatible with the uniform expectation when the normal error model is used. However, severe discrepancies between the uniform expectation and the LR-test are observed in the upper panel of Fig. 2. In contrast, the F-test shows no indication for a deviation from the uniform expectation when the normal error model is not true. The bias towards low p-values of the LR-test in the Student error model is equivalent to a bias towards larger significances. Similar Monte Carlo simulations for all other γ-ray sources that were selected for the analysis (see Sect. 3) confirm the observation that the LR-test is not robust against violations of the normal error model in the HGPS catalog.

Fig. 2. Distribution of p-values obtained in LR- and F-tests in a Monte Carlo simulation of true null hypotheses. The expectation for the p-value distribution is uniform in the interval 0 ≤ p ≤ 0.5 where 50% of the p-values are predicted. The other half of the p-values are expected at p = 1, which is due to the constraint on positive energy cutoffs (see also Appendix A). Upper panel: resulting p-value distribution in the range 0 ≤ p ≤ 0.5 when the simulated spectral points disperse around the true data according to a Student t-distribution with 5 degrees of freedom. It is clear that the LR-test (red) deviates from the uniform expectation for small p-values. This means that the LR-test severely overestimates the significance in this case. The distribution of p-values for the F-test (blue) is compatible with the uniform expectation. The distribution of p-values in the lower panel is generated when the simulated spectral points are normally distributed around their true value. In this case, the p-value distribution for both tests is compatible with the uniform expectation. The parameters of the Monte Carlo simulation are discussed in Sect. 4.2.

The improved robustness of the F-test compared to the LR-test comes at the price of a reduced test power. The test power can be estimated with a toy Monte Carlo simulation of true alternative hypotheses. For this, five spectral points are calculated at energies E_i = ξⁱE_min where ξ = 2.4 and E_min = 0.4 TeV are the median values of the binning factor and the minimal energy bin for the selected HGPS data (see Table 1). At each energy, the flux is sampled from ϕ₀(E/E₀)^−Γ exp(−E/E_cut,  γ) where ϕ₀ = 3.51 × 10⁻¹² cm⁻² s⁻¹ TeV⁻¹ and Γ = 2.11 are the median flux normalization and power-law index of the selected HGPS data (see Table 1). The flux normalization energy E₀ = 1 TeV is fixed and the cutoff energy E_cut,  γ is varied. The error of each sampled flux point is set to 20% of the flux in this toy simulation. Figure 3 compares the energy cutoff detection power between the LR- and F-tests. It is clearly seen that, compared to the LR-test, the F-test has less or at most equal power for all considered true energy cutoffs. For example, for a true energy cutoff of 100 TeV, the LR-test has 10% more power than the F-test.

Fig. 3. Simulated coverage of the interval I given by Eq. (8) at 90% CL (blue). The coverage is compatible with the nominal confidence level for true energy cutoffs below 50 TeV where the F-test power (magenta) is large. The interval I overcovers for large true energy cutoffs where the F-test power is small. The LR-test power is shown in yellow for comparison.

Despite the reduced power of the F-test when compared to the LR-test, the possible influence of systematic effects in HGPS data motivates the use of the F-test in the following discussion.

4.3. Lower limit on the energy cutoff

In practice, a hypothesis test for the presence of a cutoff is performed at high confidence level, for example at a significance level of 5σ. If the power-law hypothesis is not discarded, a lower limit on the energy cutoff E_cut is required at reduced confidence level (e.g., 90%) to constrain theoretical models. Depending on the alternative hypothesis, E_cut can refer to E_cut,  γ (Eq. (3)) or to E_cut,  h (Eq. (1)). A lower limit $E_{\text{cut}}^{-}$ on the energy cutoff E_cut corresponds to a one-sided confidence interval

$$\begin{array}{c}\hfill I=(E_{\mathrm{cut}}^{-},\infty ).\end{array}$$(8)

$E_{\text{cut}}^{\text{MLE}}$ is in the following used to refer to the maximum likelihood estimate of the energy cutoff as determined in a least-square fit. The lower limit $E_{\text{cut}}^{-}$ is the solution to the equation

$$\begin{array}{c}\hfill F(E)=F_{\mathrm{crit}},\end{array}$$(9)

in the interval $0<E<E_{\text{cut}}^{\text{MLE}}$ where

$$\begin{array}{c}\hfill F(E)=(N-3)\frac{\mathrm{RSS}({\widehat{\mathit{\theta }}}_0^{\ast }|E)-\mathrm{RSS}({\widehat{\mathit{\theta }}}_1)}{\mathrm{RSS}({\widehat{\mathit{\theta }}}_1)},\end{array}$$(10)

and RSS is given by Eq. (5). The parameters ${\widehat{\mathit{\theta }}}_0^{\ast }$ are the maximum likelihood estimators of ϕ₀ and Γ for a power-law model with fixed energy cutoff at E. Similarly, ${\widehat{\mathit{\theta }}}_1$ is the maximum likelihood estimator for the parameters of a power-law model with variable exponential energy cutoff, again with the constraint that E_cut >  0. The critical value F_crit = F_crit(CL) is given by Eq. (7).

Equation (10) is the result of an F-test of the hypothesis H₀ :   E_cut = E against the alternative hypothesis H₁ :   E_cut ≠ E and E_cut >  0. The existence of $E_{\text{cut}}^{-}$, defined in Eq. (8), is shown in Appendix B. In the following, it is assumed that $E_{\text{cut}}^{-}$ is the unique solution of Eq. (10) in the energy interval $0<E<E_{\text{cut}}^{\text{MLE}}$ (see also the discussion in Appendix B).

It is shown in Appendix C that the interval I is expected to overcover, that is, the probability to contain the true value is larger than the desired confidence level, when the true energy cutoff is large. However, excellent frequentist coverage properties are expected when the true energy cutoff is small such that the F-test power at the given CL is large. Figure 3 shows the result of a Monte Carlo simulation where γ-ray energy spectra with energy cutoff shape β = 1 were simulated following the respective discussion in Sect. 4.2. The result of this simulation shows that the frequentist coverage of the interval I agrees with the nominal 90% confidence level when the true energy cutoff is smaller than 50 TeV. For energy cutoffs larger than 50 TeV, the test power at the confidence level of 90% decreases fast and the interval overcovers.

It is concluded that the lower limit $E_{\text{cut}}^{-}$ defined in Eq. (8) has very good coverage properties when the true cutoff of the γ-ray spectrum is small. Otherwise it is conservative. For typical source and instrumental parameters representative of the HESS Galactic plane survey, good coverage can be expected when the true cutoff of the γ-ray energy spectrum is below 50 TeV.

4.4. Systematic errors

The lower limit on the energy cutoff is the lower end of the confidence interval I given by Eq. (8). This is a confidence interval for the maximum likelihood estimator $E_{\text{cut}}^{\text{MLE}}$ of the energy cutoff. Systematic effects on the energy cutoff become dominant when the systematic variation of $E_{\text{cut}}^{\text{MLE}}$ is not covered by the confidence interval I. For data from Cherenkov telescope arrays like HESS and VERITAS, possible reasons for systematic variations of $E_{\text{cut}}^{\text{MLE}}$ can be atmospheric and analysis effects or instrumental problems like broken camera elements. In case of the HGPS, the systematic errors on the fit parameters of a power-law model are estimated to be 30% on the flux normalization and 0.2 on the power-law index (Abdalla et al. 2018). These results were obtained from the comparison of the fit results obtained with different analysis methods. However, no public information is available on the systematic variation of $E_{\text{cut}}^{\text{MLE}}$ in the HGPS.

5. Analysis of HGPS data

The analysis of the HGPS data, after the quality selection discussed in Sect. 3.2, starts with a consistency check. It is ensured that the HGPS spectral points can be fit with results that are in reasonable agreement with the HGPS analysis. Afterwards, lower limits on the spectral cutoffs are derived for the selected γ-ray sources.

The HGPS catalog contains best-fitting power-law parameters for each source. The fit in the HGPS analysis is performed with an unbinned forward folding method (Piron et al. 2001). In contrast, binned least-square fits of the public spectral HGPS data are used in this analysis because unbinned data are not publicly available. However, the best least-square fit parameters are in reasonable agreement with the cataloged best-fit parameters. The maximum absolute difference between the power-law index obtained in this analysis and the power-law index given in the HGPS catalog is 1.8σ for HESS J1708−443 where σ denotes the HGPS catalog index error. The median absolute difference between the fitted power-law indices is 0.6σ. The corresponding largest absolute difference between the flux normalizations is 1.2σ for HESS J1026−582. The median absolute flux difference for the dataset is 0.1σ.

Table 2 shows the results of a search for spectral cutoffs for the 25 selected γ-ray sources in the HGPS. The best-fit exponential cutoff models (Eq. (3) with β = 1 and Eq. (1)) have a p-value larger than 1% in a χ² test for all considered spectra. All calculated solutions to Eq. (10) are unique. No energy cutoff is detected with a statistical significance exceeding 5σ for any of the analyzed sources. The lower limit on the energy cutoff in a hadronic scenario is larger than 100 TeV for five sources. For these five sources, the best fit power-law index is in the range between 2.0 and 2.5. Figure 1 shows the spectral data together with the best-fit power-law and the excluded cutoffs for the sources HESS J1800−240 and HESS J1746−285. For these sources, the most and least constraining lower limit on the energy cutoff is derived in this analysis.

No systematic error on the lower limit on the energy cutoff can be derived based on HGPS data alone since the systematic variation of $E_{\text{cut}}^{\text{MLE}}$ is not public. However, in the case of the source HESS J1908+063, an at least partially independent dataset is available. This dataset is discussed in the following section.

6. Analysis of MGRO J1908+06 data

The γ-ray source MGRO J1908+06 is associated with HESS J1908+063. This γ-ray source is one of the unidentified HGPS sources for which the lower limit on the energy cutoff in a hadronic scenario is found to be larger than 100 TeV in Sect. 5. Compared to the HGPS catalog, additional public data for this γ-ray source are available. These data, acquired with different experiments, are described and analyzed in this section.

6.1. Data for MGRO J1908+06

Milagro detected the γ-ray source at a median energy of 20 TeV with a flux of (8.8 ± 2.4) × 10⁻¹⁵ TeV⁻¹ cm⁻² s⁻¹ (Abdo et al. 2007). HESS spectral data are published in Aharonian et al. (2009b), based on the analysis of 27 h of data. This dataset will be denoted by HESSd1 in the following while HESSd2 will be used to identify the spectral HGPS data on MGRO J1908+06. The 27 h HESSd1 dataset must be at least partially independent from the 12 h HESSd2 because of the greater livetime. VERITAS also observed MGRO J1908+063 and collected 62 h of data between the years 2007 and 2012. Spectral datapoints derived from the VERITAS dataset are published in Aliu et al. (2014). The spectral properties of the HESSd1, VERITAS, and Milagro datasets are compatible with each other (Aliu et al. 2014). Two lower limits on the energy cutoff E_cut,  γ of the γ-ray spectrum were derived before. A 90% CL lower limit on E_cut,  γ is derived from a combination of the VERITAS and Milagro data at 17.7 TeV (Aliu et al. 2014). Also at 90% CL, the combination of the HESSd1 and Milagro data leads to a lower limit of E_cut,  γ = 19.1 TeV (Aharonian et al. 2009b).

The additional Cherenkov telescope data for MGRO J1908+063 from HESS and VERITAS must pass the same quality selection as discussed for the HGPS spectral data in Sect. 3.2. The HESSd1 dataset passes this selection. However, the public spectral data for MGRO J1908+063 from VERITAS is binned with a factor ξ = 1.4 while the energy resolution is similar to that of HESS (Holder 2011). Only every second spectral point from the public VERITAS data is therefore used in this analysis. The following analysis results are checked with the full VERITAS dataset to ensure that no bias towards better results is introduced with this spectral point selection.

6.2. Data analysis

The analysis of the selected spectral data from HESSd1, VERITAS, and Milagro results in a best-fit power-law spectrum with ϕ₀ = (4.2 ± 0.2) × 10⁻¹² TeV⁻¹ cm⁻² s⁻¹ and Γ = 2.17 ± 0.04. The power-law index is compatible with the index Γ = 2.26 ± 0.06 cataloged in the HGPS while the flux normalization in the HGPS, (1.05 ± 0.09) × 10⁻¹¹ TeV⁻¹ cm⁻² s⁻¹, is incompatible. The fit to a power-law model has a p-value of 0.91 in a χ² goodness-of-fit test. The lower 90% CL limits on the spectral energy cutoffs derived from this dataset are E_cut,  γ = 29.5 TeV and E_cut,  h = 141.5 TeV. These limits confirm the limits derived from the HGPS data. Overall, the agreement between this analysis and the HGPS analysis result is very good with regards to the derived lower limits on spectral energy cutoffs. Data and best-fitting models are shown in Fig. 4.

Fig. 4. Spectral γ-ray data for the source MGRO J1908+06. Shown as green circles are the HESS data from Aharonian et al. (2009b). Shown in cyan open circles are the spectral data from the HGPS. Blue spectral points are from VERITAS (Aliu et al. 2014). VERITAS data marked with a square are not selected for the analysis to avoid correlations between spectral points. The Milagro data point is from Abdo et al. (2007). Shown in black, yellow, and red are best fits to all datapoints indicated by filled circles. The black line is a fit to a power-law model. Yellow and red lines are best fits to a power-law model with cutoff at the 90% CL lower limits Ecut,  γ = 29.5 TeV (yellow) and Ecut,  h = 141.5 TeV (red).

Compared to the HESSd2 dataset with 6 spectral points, the combined data from HESSd1, VERITAS, and Milagro with 15 spectral points considered in this section contain many more flux measurements. This enables the search for a spectral break in the TeV energy range. Multiple F-tests are performed where the fit to the data of a power-law model (H₀ :   ΔΓ = 0) and a broken power-law (H₁ :   ΔΓ >  0, see Eq. (2)) are compared. The comparison cannot be performed in a single test because the break energy E_break is not defined under the null hypothesis. Using E_i again to denote the energies of the N available spectral flux data points, the energy range E₂ = 0.75 TeV to E_N − 2 = 14.8 TeV is scanned for a break energy E_break with a logarithmic binning factor of 1.05. The F-test for the comparison of the fit quality of the power-law model and the broken power-law model, assuming the respective energy break energy, is performed in each scanning step. The F-test is inverted to derive an upper limit on the index change ΔΓ as a function of the assumed break energy E_break. The result is shown in Fig. 5. An index change ΔΓ >  0.5 is ruled out at 90% CL for energies between 1 TeV and 10 TeV.

Fig. 5. Upper limit on a break in the power-law spectrum by ΔΓ as a function of the break energy. Shown as red line is the 90% CL upper limit on ΔΓ derived from the combined data discussed in Sect. 6.

The largest local significance obtained in the 61 tests of the energy break scan is p_local = 0.14 with an F-test statistic c = 1.23. The local significance must be transformed into a global significance p_global considering the number of performed tests. Based on results for the global p-value in a multiple testing scenario (see Gross & Vitells 2010, and references therein), Algeri & van Dyk (2017) derive a global p-value correction for an F-process. In the case considered here, this correction reads

$$\begin{array}{c}\hfill p_{\mathrm{global}}\le p_{\mathrm{local}}+{\left(\frac{N-3+c}{N-3+c_0}\right)}^{-N/2+2}E[N_{c_0}],\end{array}$$(11)

where E[N_c0] is the expected number of upcrossings over the level c₀ of the process of F-test statistics during the energy cutoff scan when the null hypothesis is true. Equation (11) allows us to extrapolate the global p-value correction from a low test statistic level c₀, where a small number of simulations is required for the estimation of E[N_c0], to a larger test statistic value c. Following Algeri et al. (2016b), E[N_c0] is estimated in a parametric bootstrap simulation at c₀ = 0.3 to be E[N_c0] = 0.52 ± 0.02. With this correction, the global p-value of the F-test for a spectral break in the selected MGRO J1908+06 data from HESS, VERITAS, and Milagro can be estimated to be p_global = 0.5. The power-law hypothesis is therefore not discarded in favor of a broken power-law.

7. Discussion

The analysis of HGPS data resulted in five γ-ray sources for which the energy spectrum is compatible with a power-law, and in a hadronic scenario the lower limit on the energy cutoff of the accelerated particle population is larger than 100 TeV. These five γ-ray sources are discussed in detail in this section. An emphasis of the discussion is on the plausibility of the hadronic scenario for these sources.

7.1. HESS J1800−240

The region around HESS J1800−240 is found to have a complex γ-ray morphology in Aharonian et al. (2008) based on a HESS dataset of 42 h livetime. Four γ-ray sources are detected in this region. The SNR W28 is associated with the source HESS J1801−233. The region around HESS J1800−240 itself is subdivided into three hotspots that are spatially coincident with molecular clouds.

Based on a smaller dataset of 10 h, the spectral HGPS analysis detects only one γ-ray source in the region. The lower limit on the energy cutoff for accelerated particles derived above is E_cut,  h = 510 TeV at 90% CL. This constraint is comparable to the lower limit of E_cut,  h = 600 TeV at 90% CL on the energy cutoff of particles in the vicinity of the Galactic center in a hadronic scenario (Abramowski et al. 2016).

A physics case for the association of this γ-ray source with a Pevatron is the model discussed in Gabici & Aharonian (2007) where delayed CRs from a SNR illuminate a molecular cloud. The plausible physics case and the interesting constraint on the energy cutoff in a hadronic scenario may motivate a systematic investigation of differences between the HGPS analysis and the analysis in Aharonian et al. (2008) as well as a deeper exposure of the region with current Cherenkov telescopes.

7.2. HESS J1641−463

The γ-ray source is discussed in Abramowski et al. (2014), where the TeV spectrum extracted from 72 h of HESS data is found to be compatible with a power-law model. Additionally, a lower limit on the energy cutoff at E_cut,  h = 100 TeV is derived at 99% CL in a hadronic model. The lower limit of E_cut,  h = 185 TeV at 90% CL on the energy cutoff derived in this work from 27 h of HGPS data is compatible with the result in Abramowski et al. (2014).

A detailed hadronic model for the γ-ray emission of HESS J1651−463 is presented in Tang et al. (2015). Within this model, runaway CR particles accelerated inside the young and nearby SNR G338.3−0.0 interact with a molecular cloud where γ-rays are produced via the decay of neutral pions. The SNR G338.3−0.0 itself is associated with the γ-ray source HESS J1640−465.

Given the presence of a nearby SNR and a molecular cloud, a hadronic scenario for this γ-ray source is plausible. Observations with instruments which have an improved sensitivity, at a few tens of TeV to a few hundred TeV, or a neutrino detection are necessary to confirm a Pevatron scenario.

7.3. MGRO J1908+06

The γ-ray source MGRO J1908+06 was discovered by Milagro (Abdo et al. 2007) and later confirmed by HESS (Aharonian et al. 2009b). Abdo et al. (2009) find a spatial coincidence between MGRO J1908+06 and the PWN of the pulsar PSR 1907.5+0602 detected by Fermi/LAT. However, three significant emission regions are detected in the field with data from VERITAS (Aliu et al. 2014). This raises the question as to whether the entire emission can originate from the PWN or an additional γ-ray source is present. The SNR G40.5−0.5 is discussed in Aliu et al. (2014) as a candidate for an additional γ-ray source.

The HESS spectrum is confirmed by the measurement with VERITAS (Aliu et al. 2014). However, the spectrum measured by the Astrophysical Radiation with Ground-based Observatory (ARGO-YBJ) is incompatible with the spectra measured by HESS and VERITAS (Bartoli et al. 2012). The High-Altitude Water Cherenkov Gamma-Ray Observatory (HAWC) detects the γ-ray source at an energy of 7 TeV and derives a flux which varies within a factor of 2.5 depending on the assumed source extension (Abeysekara et al. 2017). The γ-ray source is discussed as a possible source of high-energy neutrinos (Halzen et al. 2017). Aartsen et al. (2017b) find a pre-trial p-value of 2.5% when testing a background-only hypothesis against the emission of astrophysical high-energy neutrinos from MGRO J1908+06 based on data from IceCube.

It is concluded that the spectral properties of the γ-ray emission of MGRO J1908+06 are currently not unambiguously measured. Moreover, the source identification is not finally resolved, although an at least partial association with the PWN of PSR 1907.05+0602 is likely.

The lower limit on the exponential cutoff of the γ-ray spectrum derived above, E_cut,  γ ≈ 30 TeV, is a factor of almost two larger than previous results discussed in Aharonian et al. (2009b) and Aliu et al. (2014). The lower limit E_cut,  γ ≈ 30 TeV is derived from HGPS data (E_cut,  γ = 31.1 TeV at 90% CL) and confirmed in an analysis of combined data from HESS, VERITAS, and Milagro (E_cut,  γ = 29.5 TeV at 90% CL). The compatibility of the lower limits on the energy cutoff derived from two different datasets shows that the systematic error on the energy cutoff is not likely to be dominant in the case of this source.

A sharp energy break is ruled out at 90% CL in the energy range between 1 TeV and 10 TeV (see Fig 5). A search for such an energy break could, for example, be motivated by the sharp energy break observed around 1 TeV in the electron spectrum measured on earth (see e.g., Aharonian et al. 2009a. The power-law index of the electron spectrum measured on earth changes within measured errors by ΔΓ_e = 1 at E_break = 1 TeV. This energy break can be interpreted as a cooling break (Recchia et al. 2019). When electrons generate the γ-ray signal detected towards MGRO J1908+06 via inverse Compton scattering in the Thomson regime, a break in the spectrum of electrons by ΔΓ_e = 1 would translate into a change of the γ-ray spectrum index by ΔΓ = 0.5 (see e.g., Hinton & Hofmann 2009). An energy break of this kind is ruled out at 90% CL in the energy range between 1 TeV and 10 TeV.

Although the derived lower limits on the γ-ray energy cutoff are more constraining than previous lower limits and a cooling break can be ruled out in the energy range between 1 and 10 TeV, the results cannot rule out leptonic scenarios for the γ-ray emission. Deeper exposures and observations with more sensitive instruments are necessary to reveal the nature of this source.

7.4. HESS J1852−000

It is discussed in Abdalla et al. (2018) that the spectral HGPS data for this source must be treated with caution due to deviations between the main HGPS result and an independent cross check data analysis. HESS data for this source have been analyzed before (Kosack et al. 2011), but an energy spectrum has not been derived. No other spectral data in the TeV γ-ray energy range are publicly available. The analysis of the spectral cutoff for this source is therefore not conclusive.

Different scenarios for the production of γ-rays are discussed in Kosack et al. (2011) and Bamba et al. (2016). Among them is an association with the SNR Kes 78 and a nearby molecular cloud.

It is concluded that a hadronic source scenario in which runaway CRs from the SNR Kes 78 illuminate a molecular cloud can currently not be ruled out. The inconclusive spectral data may motivate a more detailed analysis and possibly further observations with current-generation Cherenkov telescopes.

7.5. HESS J1634−472

HESS J1634−472 is cataloged as an unidentified source in the HGPS without further discussion (Abdalla et al. 2018). The source was also detected in a previous survey of the inner Galaxy with HESS (Aharonian et al. 2006a). Here, the SNR G337.2+0.1 and a source of X-rays are found to be within 0.2° of HESS J1634−472 but no association is claimed. Acero et al. (2013) discuss a detection of the source with Fermi/LAT and find that there is no counterpart pulsar in this region which is energetic enough to power the γ-ray source HESS J1634−472.

No further public TeV γ-ray data are available for this source. The interesting constraint on E_cut,  h >  108 TeV at 90% CL derived from the analysis of HGPS data above may motivate a dedicated re-analysis of the HESS data and possibly a deeper exposure of the source.

8. Conclusion

Five γ-ray sources are found in the HGPS catalog for which the maximal energy of accelerated particles is at least 100 TeV in a hadronic model. For at least three of these sources, a hadronic scenario for the γ-ray emission is plausible, that is, as a result of the interaction of runaway CRs with a molecular cloud. One of the Pevatron candidates found in the HGPS is MGRO J1908+06. The γ-ray spectrum of this source extends to at least E_cut,  γ = 30 TeV without indication of a cutoff. The lower limit on E_cut,  γ for this source is a factor of almost two larger than previous results. A break ΔΓ >  0.5 for this γ-ray source can be ruled out at 90% CL in the energy range between 1 and 10 TeV.

The search presented in this work can only find Pevatron candidates. A conclusive identification must be performed with data from neutrino telescopes (Aartsen et al. 2017a) or γ-ray data around and above 100 TeV. The extended air shower array HAWC (DeYoung 2012), now operating for more than 3 years, can add spectral γ-ray data at a few tens of TeV. An improved sensitivity and higher energies will be accessible with the upcoming γ-ray detectors, for example the Cherenkov Telescope Array (CTA; The CTA Consortium 2011) and the Large High Altitude Air Shower Observatory (LHAASO; Sciascio 2016). Finally, the development of new experimental techniques, such as for the Tunka Advanced Instrument for Gamma-ray and Cosmic ray Astrophysics (TAIGA; Tluczykont et al. 2014), will help to identify Galactic Pevatrons, thereby aiding the quest for the origin of cosmic rays.

Acknowledgments

The author acknowledges financial support by the German Ministry for Education and Research (BMBF).

References

Appendix A: Critical values

Consider, as in Sect. 4.1, the two nested hypotheses H₀ :   λ = 0 and H₁ :   λ >  0 for the single parameter λ. The parameter set for H₀ is not in the interior of the parameter set of H₁. The asymptotic result for the LR test statistic −2lnΛ under H₀ discussed in Wilks (1935) is therefore not applicable. However, a modification is discussed in Chernoff (1954) where it is predicted that −2lnΛ is asymptotically distributed like an equal mixture of two random variables when H₀ is true. In the considered case, one of the random variables is χ²-distributed with one degree of freedom. The other random variable has a Dirac delta function probability density function δ(0) (see also Algeri et al. 2016a). With this asymptotic distribution f_H0 of the test statistic −2lnΛ under H₀, the critical value Λ_crit of the LR test statistic can be inferred from ${\int }_{-\infty }^{{\mathrm{\Lambda }}_{\text{crit}}}\text{d}x{f}_{\text{H0}}(x)=\text{CL}$. Let $F_{{\chi }_{\text{1}}^{\text{2}}}$ denote the cumulative density functions of a χ²-distribution with one degree of freedom. Subsequently, ${\int }_{-\infty }^{{\mathrm{\Lambda }}_{\text{crit}}}\text{d}x{f}_{\text{H0}}(x)=1/2F_{{\chi }_{\text{1}}^{\text{2}}}({\mathrm{\Lambda }}_{\text{crit}})+1/2{\int }_{-\infty }^{{\Lambda }_{\text{crit}}}\text{d}x\delta (0)=\text{CL}$. This leads to the equality ${\mathrm{\Lambda }}_{\text{crit}}=F_{{\chi }_{\text{1}}^{\text{2}}}^{-1}(2\text{CL}-1)$ which can be used to calculate the critical value Λ_crit of the LR test statistic at a given CL.

A similar result follows for the critical value of the asymptotic F statistic in Sect. 4.1.2: the F statistic (Eq. (6)) and the LR statistic −2lnΛ are related via $F=(N-3)(-2ln\mathrm{\Lambda })/\mathrm{RSS}({\widehat{\mathit{\theta }}}_1)$ (see Sect. 4.1.1). Let now F_{1, N − 3} be the cumulative density function for a F-distributed random variable with 1 and N − 3 degrees of freedom. Because the LR statistic is asymptotically expected to be zero in half of the tests when the null hypothesis is true, the critical value F_crit of the F-test can again be inferred from $1/2F_{1,N-3}({\mathrm{\Lambda }}_{\text{crit}})+1/2{\int }_{-\infty }^{{\mathrm{\Lambda }}_{\text{crit}}}\text{d}x\delta (0)=\text{CL}$ to be $F_{\mathrm{crit}}=F_{1,N-3}^{-1}(2\mathrm{CL}-1)$.

Appendix B: Existence of the lower limit

Let F(E) and F_crit be given by Eqs. (7) and (10) and. The lower limit on the energy cutoff is defined as a solution to the equation

$$\begin{array}{c}\hfill F(E)=F_{\mathrm{crit}},\end{array}$$(B.1)

in the interval $(0,E_{\text{cut}}^{\text{MLE}})$. Here, $E_{\text{cut}}^{\text{MLE}}$ is the maximum likelihood estimator for the energy cutoff (see Sect. 4.3). It holds that F_crit >  0 and F(E) is a continuous function. A solution to Eq. (B.1) does always exist in the given interval because $F(E_{\text{cut}}^{\text{MLE}})=0$ and F(E) → ∞ for E → 0 by definition of F(E).

However, the solution to Eq. (B.1) in the interval $(0,E_{\text{cut}}^{\text{MLE}})$ must not necessarily be unique when general data is considered. Multiple solutions can, for example, occur in case of a mismatch between the model with an exponential cutoff and the true model where the likelihood function has multiple local maxima. In this case, the meaning of a lower limit on the energy cutoff is questionable.

In any case, special care must be taken when multiple solutions to Eq. (B.1) are found in the interval $(0,E_{\text{cut}}^{\text{MLE}})$. In the data analyzed in this work, all lower limits were unique when using a confidence level of 90%.

Appendix C: Frequentist coverage of the lower limit

Consider the interval I given by Eq. (8). It is argued that I has excellent frequentist coverage properties when the true energy cutoff is small such that the test power is large. When the true energy cutoff is large such that the test power is small, the interval I overcovers and the lower limit on the energy cutoff is conservative.

In the following, F, F(E), and F_crit are given by Eqs. (6), (7), and (10). It is assumed that only one solution $E_{\text{cut}}^{-}$ for Eq. (B.1) in the interval $(0,E_{\text{cut}}^{\text{MLE}})$ exists. Similarly, it is also assumed that at most one solution to Eq. (B.1) exists in the interval $(E_{\text{cut}}^{\text{MLE}},\infty )$. Moreover, $E_{\text{cut}}^{+}$ is the solution to Eq. (B.1) in $(E_{\text{cut}}^{\text{MLE}},\infty )$, if it exists, and $E_{\text{cut}}^{+}=\infty$ if no solution exists. Now consider the interval $(E_{\text{cut}}^{-},E_{\text{cut}}^{+})$. This interval is the acceptance interval A = {E|F(E) ≤ F_crit} of an F-test, and as such defines a confidence interval for $E_{\text{cut}}^{\text{MLE}}$. However, the confidence level of the F-test to which A is the acceptance interval depends on E. This is a result of the fit being constrained to λ >  0.

When the true energy cutoff is large, F(E) will need to be evaluated at large energies E to obtain A. By definition of F(E), it holds that F(E)→F when E → ∞ and F(E) is the F-test statistic for a test of a power-law hypothesis against a power-law model with exponential cutoff. This means that A is, in the limit E → ∞, the acceptance interval for an F-test constructed at confidence level CL. In other words, it holds that P(E_true ∈ A) = CL. Additionally, by definition of the interval I, it holds that P(E_true ∈ I) ≥ P(E_true ∈ A). Together it follows that P(E_true ∈ I) ≥ CL, which is the overcoverage of I as a frequentist confidence interval when the true energy cutoff is large.

Now, let the true energy cutoff be small but positive such that the power-law hypothesis is discarded at the given CL. The fit constraint λ >  0 is in this case irrelevant because the best fit λ is always large. By definition of F_crit, it holds that P(E_true ∈ A) = 2CL − 1. This is because without the constraint λ >  0, A is the acceptance interval of an F-test that is constructed at confidence level 2CL − 1. Close to $E_{\text{cut}}^{\text{MLE}}$, F(E) can be quadratically approximated; see also Fig. C.2. Using the nomenclature of Sect. 4.1, it holds that $F(E)\approx (N-3){(E-E_{\mathrm{cut}}^{\mathrm{MLE}})}^2/(\mathrm{RSS}({\widehat{\mathit{\theta }}}_1){\sigma }^2)$ where σ² is the inverse Fisher information of L(θ₁) evaluated at ${\widehat{\mathit{\theta }}}_1$. This means that the confidence interval A will be symmetric around $E_{\text{cut}}^{\text{MLE}}$ when E is small such that the solution to Eq. (B.1) is found in the range where the quadratic approximation holds. It follows that $P(E_{\text{true}}\in (0,E_{\text{cut}}^{-}))=P(E_{\text{true}}\in (E_{\text{cut}}^{+},\infty ))=(1-(2\text{CL}-1))/2=1-\text{CL}$. For the coverage probability, it follows that $P(E_{\text{true}}\in I)=P(E_{\text{true}}\in A)+P(E_{\text{true}}\in (E_{\text{cut}}^{+},\infty ))=\text{CL}$. This states that very good coverage properties of I as a frequentist confidence interval are expected when the F-test power is large at the given CL.

Fig. C.1. Example for the dependence of F(E) given by Eq. (10) on the cutoff energy E. The example is based on the HGPS data for the source HESS J1858+020. The maximum likelihood estimate $E_{\text{cut},\gamma }^{\text{MLE}}$ for an exponential cutoff in the γ-ray spectrum (Eq. (3) with β = 1) is indicated with a green dashed line. The blue and the magenta lines indicate the critical value Fcrit and the inferred lower limit on Ecut,  γ at 90% CL. The magenta region is the excluded range of cutoff energies.

Fig. C.2. Same as in Fig. C.1 but based on the HGPS data for the source HESS J1507−622. The power-law hypotheses is discarded at 90% CL and the acceptance interval A = {E|F(E) ≤ Fcrit} is now a finite interval roughly symmetric around $E_{\text{cut},\gamma }^{\text{MLE}}$. Shown in yellow is the quadratic approximation of F(E) around $E_{\text{cut},\gamma }^{\text{MLE}}$.

All Tables

All Figures

Fig. 1. Spectral fits to the data from the HGPS source for which the most (upper panel) and least (lower panel) constraining lower limit on Ecut,  h is derived. Shown in green are the best fits to the data with a power-law spectrum. Red lines are best fits to the data with a γ-ray power-law spectrum with exponential cutoff at the lower limit Ecut,  γ derived from the HGPS data. Yellow lines are the corresponding best fits to the data in a hadronic model with cutoff at the lower limit Ecut,  h. Limits for the lower plot are at 90% CL. In the upper panel, best fits with lower cutoff limit at 90% CL are hardly distinguishable from the best power-law fit. Therefore, 99% CL best fits are shown in this case.

In the text

Fig. 2. Distribution of p-values obtained in LR- and F-tests in a Monte Carlo simulation of true null hypotheses. The expectation for the p-value distribution is uniform in the interval 0 ≤ p ≤ 0.5 where 50% of the p-values are predicted. The other half of the p-values are expected at p = 1, which is due to the constraint on positive energy cutoffs (see also Appendix A). Upper panel: resulting p-value distribution in the range 0 ≤ p ≤ 0.5 when the simulated spectral points disperse around the true data according to a Student t-distribution with 5 degrees of freedom. It is clear that the LR-test (red) deviates from the uniform expectation for small p-values. This means that the LR-test severely overestimates the significance in this case. The distribution of p-values for the F-test (blue) is compatible with the uniform expectation. The distribution of p-values in the lower panel is generated when the simulated spectral points are normally distributed around their true value. In this case, the p-value distribution for both tests is compatible with the uniform expectation. The parameters of the Monte Carlo simulation are discussed in Sect. 4.2.

In the text

	Fig. 3. Simulated coverage of the interval I given by Eq. (8) at 90% CL (blue). The coverage is compatible with the nominal confidence level for true energy cutoffs below 50 TeV where the F-test power (magenta) is large. The interval I overcovers for large true energy cutoffs where the F-test power is small. The LR-test power is shown in yellow for comparison.
In the text

Fig. 4. Spectral γ-ray data for the source MGRO J1908+06. Shown as green circles are the HESS data from Aharonian et al. (2009b). Shown in cyan open circles are the spectral data from the HGPS. Blue spectral points are from VERITAS (Aliu et al. 2014). VERITAS data marked with a square are not selected for the analysis to avoid correlations between spectral points. The Milagro data point is from Abdo et al. (2007). Shown in black, yellow, and red are best fits to all datapoints indicated by filled circles. The black line is a fit to a power-law model. Yellow and red lines are best fits to a power-law model with cutoff at the 90% CL lower limits Ecut,  γ = 29.5 TeV (yellow) and Ecut,  h = 141.5 TeV (red).

In the text

	Fig. 5. Upper limit on a break in the power-law spectrum by ΔΓ as a function of the break energy. Shown as red line is the 90% CL upper limit on ΔΓ derived from the combined data discussed in Sect. 6.
In the text

Fig. C.1. Example for the dependence of F(E) given by Eq. (10) on the cutoff energy E. The example is based on the HGPS data for the source HESS J1858+020. The maximum likelihood estimate $E_{\text{cut},\gamma }^{\text{MLE}}$ for an exponential cutoff in the γ-ray spectrum (Eq. (3) with β = 1) is indicated with a green dashed line. The blue and the magenta lines indicate the critical value Fcrit and the inferred lower limit on Ecut,  γ at 90% CL. The magenta region is the excluded range of cutoff energies.

In the text

	Fig. C.2. Same as in Fig. C.1 but based on the HGPS data for the source HESS J1507−622. The power-law hypotheses is discarded at 90% CL and the acceptance interval A = {E|F(E) ≤ Fcrit} is now a finite interval roughly symmetric around $E_{\text{cut},\gamma }^{\text{MLE}}$. Shown in yellow is the quadratic approximation of F(E) around $E_{\text{cut},\gamma }^{\text{MLE}}$.
In the text