11institutetext: Institut für Physik, Humboldt-Universität zu Berlin,
Newtonstr. 15, 12489 Berlin, Germany
11email<EMAIL_ADDRESS>
# Search for galactic Pevatron candidates in a population of unidentified
$\gamma$-ray sources

Gerrit Spengler

(Received September 4, 2019; accepted December 10, 2019)

###### Abstract

Aims. A list of Pevatron candidates is presented to enable deeper observations
and dedicated analyses.

Methods. Lower limits on the energy cutoff for unidentified $\gamma$-ray
sources detected in the HESS galactic plane survey are derived. Additional
public data from the VERITAS, HESS and Milagro experiments are used for MGRO
J1908+06 to confirm the limit derived from the HESS galactic plane survey data
and to enable further conclusions on the presence of spectral breaks.

Results. Five Pevatron candidates are identified in the HESS galactic plane
survey. The cutoff of the $\gamma$-ray spectrum for these sources is larger
than 20 TeV at 90% confidence level. The $\gamma$-ray sources MGRO J1908+06
and HESS J1641-463, found to be Pevatron candidates in the analysis of the
HESS galactic plane survey catalog, had been discussed as Pevatron candidates
before. For MGRO J1908+06, the lower limit on the $\gamma$-ray energy cutoff
is $30$ TeV at $90\%$ confidence level. This is a factor of almost two larger
than previous results. Additionally, a break in the $\gamma$-ray spectrum at
energies between $1$ TeV and $10$ TeV with an index change $\Delta\Gamma>0.5$
can be excluded at $90\%$ confidence level for MGRO J1908+06.
The energy cutoff of accelerated particles is larger than $100$ TeV at $90\%$
confidence level in a hadronic scenario for all five Pevatron candidates. A
hadronic scenario is plausible for at least three of the Pevatron candidates,
based on the presence of nearby molecular clouds and supernova remnants.

###### Key Words.:

astroparticle physics – gamma rays: general – methods: statistical

##  1 Introduction

Theoretical models for the origin of cosmic rays (CRs) consider young
supernova remnants (SNRs) as Pevatrons, i.e. 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 $\gamma$-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 it is in particular unclear whether SNRs are sources of PeV
CRs.
The association of a $\gamma$-ray source with a Pevatron can be tested with a
measurement of the $\gamma$-ray spectrum above TeV energies. An exponential
cutoff $E_{\mathrm{cut,\;\gamma}}$ in the $\gamma$-ray spectrum at a few $100$
TeV is expected when CRs are accelerated up to PeV energies and produce
$\gamma$-rays in interactions with ambient material (Gabici & Aharonian,
2018). The Pevatron nature of a $\gamma$-ray source is in turn constrained
when an exponential cutoff at an energy well below $100$ TeV is detected in
the observed $\gamma$-ray spectrum of the source. For example, the recent
detection of an exponential cutoff at $E_{\mathrm{cut,\;\gamma}}=3.5$ TeV in
the $\gamma$-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_{\mathrm{cut,\;h}}$
of galactic hadron accelerators are currently derived from the diffuse
$\gamma$-ray emission in the vicinity of the radio source SgrA* and for the
$\gamma$-ray source HESS J1641-463. Assuming that hadronic particles generate
the $\gamma$-ray emission from these sources, $E_{\mathrm{cut,\;h}}>600$ TeV
and $E_{\mathrm{cut,\;h}}>100$ TeV is inferred at $90\%$ and $99\%$ confidence
level (CL) respectively (Abramowski et al., 2016; Abramowski et al., 2014).
More than $100$ isolated sources of TeV $\gamma$-rays are currently known.
However, no systematic search for Pevatron candidates among the known
$\gamma$-ray sources has been presented yet.
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 HESS galactic plane survey (HGPS). This catalog
is introduced in Sec. 3. A method to constrain energy cutoffs in measurements
of powerlaw energy spectra is discussed in Sec. 4. The method is argued to be
applicable to the HGPS catalog introduced before. A search for Pevatron
candidates in the HGPS catalog is summarized in Sec. 5. For one Pevatron
candidate, identified in the previous analysis of HGPS data, additional public
data is available from the HESS, VERITAS and Milagro experiments. A combined
analysis of these data is presented in Sec. 6. Finally, sections 7 and 8
discuss the results and summarize the conclusions.

##  2 Pevatron candidate search method

The detection of TeV $\gamma$ rays from an isolated region of the sky can
indicate cosmic sites where particles are accelerated to energies larger than
$1$ TeV. However, the identification of the primary particle type is non-
trivial when only TeV $\gamma$-ray data are available. A leptonic and a
hadronic scenario are typically conceivable. Electrons are accelerated in the
leptonic scenario and generate $\gamma$-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. The
unambiguous exclusion of the leptonic scenario for a given $\gamma$-ray source
is, however, very challenging. For example, as recently discussed, even the
production of $\gamma$-rays with an energy exceeding $100$ TeV is possible in
a leptonic scenario (Amenomori et al., 2019). The non-detection of a
$\gamma$-ray energy cutoff at energies below few $100$ TeV is therefore not
considered as a sufficient criteria for the identification of a hadronic
Pevatron. Instead of an unambiguous identification of hadronic Pevatrons, the
following search aims to find Pevatron candidates, i.e. $\gamma$-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 e.g. searches for high energy neutrinos
as fingerprints of hadronic interactions.
Many known galactic $\gamma$-ray sources are identified as pulsar wind nebulae
(PWNe) for which the $\gamma$-ray emission is modeled in a leptonic scenario
(Gaensler & Slane, 2006). Another large class of $\gamma$-ray sources is not
yet identified with sources of radiation in other energy ranges. The search
for new Pevatron candidates will be restricted to these unidentified
$\gamma$-ray sources in the galactic disc. The motivation for the restriction
to unidentified $\gamma$-ray sources is at least twofold. The exclusion of
PWNe removes a large class of $\gamma$-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, i.e. less than
$1000$ years. Easier than the direct detection of Pevatrons might therefore be
the indirect identification via the detection of delayed $\gamma$-ray emission
from molecular clouds in the vicinity of SNRs. Gabici & Aharonian (2007)
discuss this scenario and motivate that unidentified $\gamma$-ray sources can
be associated with molecular clouds illuminated by nearby SNRs.
Following Kelner et al. (2006), the spectrum of $\gamma$ rays generated in the
vicinity of a Pevatron in interactions of accelerated CRs with ambient
material is in the following modeled by

$\phi(E)=\phi_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma}\exp(-\sqrt{16\lambda_{\mathrm{h}}E})\>\mathrm{.}$
(1)

Here, $\phi_{0}\exp(-\sqrt{16\lambda_{\mathrm{h}}E_{0}})$ is the flux
normalization at energy $E_{0}$ and the parent CR spectrum has a cutoff energy
of $E_{\mathrm{cut,\;h}}=1/\lambda_{h}$ (Kelner et al., 2006). Models for
diffuse shock acceleration predict a powerlaw index $\Gamma\approx 2$ (see
e.g. Hinton & Hofmann (2009)). Following this model, TeV $\gamma$-ray sources
that are selected as Pevatron candidates must have a spectrum at TeV energies
which is compatible with a powerlaw with index $\Gamma=2$. Additionally, the
inferred lower limit on the energy cutoff $E_{\mathrm{cut,\;h}}$ in the
hadronic model given by Eq. 1 must be at least
$\mathrm{O}(100\;\mathrm{TeV})$, i.e. in the order of magnitude of the current
best constraints derived for HESS J1641-463 and the $\gamma$-ray emission in
the vicinity of SgrA*.
To further constrain a leptonic scenario for the generation of $\gamma$ rays,
the presence of a break in the TeV $\gamma$-ray spectrum where the powerlaw
index changes can be excluded. Consider, for example, a leptonic scenario in
which inverse Compton losses dominate over synchrotron cooling. A Klein-
Nishina (KN) break is expected in the $\gamma$-ray spectrum at an energy
$E_{\mathrm{KN}}$ which depends on the target photon field. For cosmic
microwave background (CMB) target photons, the break is expected at 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 $\gamma$-ray spectrum at energies far away from $E_{\mathrm{KN}}$
can be described by a broken powerlaw with
$E_{\mathrm{break}}=E_{\mathrm{KN}}$

$\phi(E)=\cases{\phi}_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma}&E<E_{\mathrm{break}}\\\
\phi_{0}\left(\frac{E_{\mathrm{break}}}{E_{0}}\right)^{\Delta\Gamma}\left(\frac{E}{E_{0}}\right)^{-\Gamma-\Delta\Gamma}E>E_{\mathrm{break}}\;\mathrm{.}{}$
(2)

It is expected that $\Delta\Gamma=\Gamma$ for the index change at the KN break
(see e.g. Hinton & Hofmann (2009)). The selection of sources with a powerlaw
index around $2$ at TeV energies excludes leptonic scenarios where the inverse
Compton target photons have a much larger energy than CMB photons, e.g.
optical starlight, and the KN break must occur below TeV energies. However,
leptonic scenarios with CMB target photons can typically not be constrained
with TeV $\gamma$-ray data.
The selection of a $\gamma$-ray source as a Pevatron candidate in the
following search is not sufficient to identify the source with a Pevatron.
Neither can leptonic emission scenarios be ruled out nor are the lower limits
on $E_{\mathrm{cut,\;h}}$ sufficient to argue for the acceleration of hadrons
to PeV energies. Additionally, the selection of a $\gamma$-ray source as a
Pevatron is not necessary for the association with a Pevatron. For example,
the $\gamma$-ray source Cassiopeia A would not pass the selection due to the
measured TeV cutoff discussed in Sec. 1. However, Cassiopeia A is still being
discussed as a Pevatron candidate in multi-zone models (Zhang & Liu, 2019).
Despite of the Pevatron candidate selection being neither necessary nor
sufficient for the association of a $\gamma$-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 few selected Pevatron candidates can be
performed.

##  3 HGPS data and quality selection

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

###  3.1 HGPS data

The HGPS comprises data acquired from 2004 to 2013 with the HESS I (Aharonian
et al., 2006) array of Cherenkov telescopes in Namibia. Large parts of the
Milky Way plane were observed to detect $\gamma$-ray sources in a nominal
energy range of $0.2-100$ TeV with an angular resolution better than
$0.1^{\circ}$. 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$ $\gamma$-ray sources were detected in
the plane survey. The majority ($47$) of the detected sources are
unidentified, i.e. 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 $\gamma$-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_{\mathrm{ON}}$) and
background ($N_{\mathrm{OFF}}$) events and enables a calculation of the excess
events $N_{\gamma}=N_{\mathrm{ON}}-\alpha N_{\mathrm{OFF}}$. The acceptance
normalization factor $\alpha$ is obtained from independent observations of
different sky regions without $\gamma$-ray source. Energy dependent
information on $N_{\gamma}$, $N_{\mathrm{ON}}$, $N_{\mathrm{OFF}}$ and
$\alpha$ is not part of the public HGPS catalog. However, spectral flux data
for each $\gamma$-ray source detected in the HGPS is typically available in
$6$ logarithmic energy bins. The bin centers $E_{i}$, $i=0,\dots,N-1$, of
spectral points are related via $E_{i}=\xi^{i}E_{\mathrm{min}}$ with a binning
factor $\xi$ where $E_{\mathrm{min}}$ is the center of the lowest energy bin.
Statistical errors on the $\gamma$-ray flux are given as asymmetric $68\%$
confidence level (CL) intervals and originate from the Poisson counting error
on $N_{\mathrm{ON}}$ and $N_{\mathrm{OFF}}$. Other error sources, e.g. an
error on the acceptance normalization, are not included in the statistical
error on the $\gamma$-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. (2006), the energy
resolution of the HGPS analysis is assumed to be $15\%$ or better. Only
sources in the HGPS with spectral energy binning factor $\xi>1.5$ are
selected.
In the HGPS catalog, data on the error of spectral points is given as lower
($\sigma_{\mathrm{low}}$) and upper ($\sigma_{\mathrm{high}}$) error. The
average error $\sigma=1/2\;(\sigma_{\mathrm{low}}+\sigma_{\mathrm{high}})$ is
used in the analysis. Spectral points with an error asymmetry
$\eta=\mathrm{max}(|\sigma-\sigma_{\mathrm{low}}|,|\sigma-\sigma_{\mathrm{high}}|)/\sigma>10\%$
are discarded from the analysis. Also discarded are spectral points $\phi$
which are detected at low significance level $S=\phi/\sigma<1.5$. Sources in
the HGPS catalog are not analyzed when less than $5$ spectral points remain
after the quality selection. In total $25$ out of $78$ sources are selected
for further analysis. A summary of the selected HGPS data is given in Table 1.

Table 1: Characterization of the HGPS spectral data for the $25$ selected $\gamma$-ray sources. The columns show the minimum, maximum and median of the distribution of the quantity in the corresponding row for the selected sources. $E_{\mathrm{min}}$ and $E_{\mathrm{max}}$ are the minimal and maximal bin-center energy of selected spectral points. The livetime is the dead-time corrected observation time as obtained from the HGPS dataset. $\Gamma$ and $\phi_{0}$ are the best fit index and flux normalization for a powerlaw model. The binning factor $\xi$, error asymmetry $\eta$ and spectral point significance $S$ are defined in Sec. 3.2. Quantity | Minimum | Maximum | Median
---|---|---|---
$E_{\mathrm{min}}$ (TeV) | 0.3 | 1.47 | 0.4
$E_{\mathrm{max}}$ (TeV) | 6.32 | 31.62 | 12.71
$\xi$ | 1.6 | 2.5 | 2.4
Livetime (h) | 4.0 | 259.0 | 17.0
$\Gamma$ | 1.81 | 2.69 | 2.31
$\phi_{0}$ ($\mathrm{cm^{-2}s^{-1}TeV^{-1}}$) | 1.80e-13 | 1.56e-11 | 3.51e-12
$\eta$ (%) | 0.0 | 8.8 | 1.7
S | 1.5 | 22.8 | 4.3

##  4 Constraints on the energy cutoff of a powerlaw spectrum

Figure 1: Spectral fits to the data from the HGPS source for which the most
(upper panel) and least (lower panel) constraining lower limit on
$E_{\mathrm{cut,\;h}}$ is derived. Shown in green are the best fits to the
data with a powerlaw spectrum. Red lines are best fits to the data with a
$\gamma$-ray powerlaw spectrum with exponential cutoff at the lower limit
$E_{\mathrm{cut,\;\gamma}}$ 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 $E_{\mathrm{cut,\;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 powerlaw fit. $99\%$ CL best fits are
therefore shown in this case.

A statistical hypothesis test can be used to search for an exponential cutoff
in a powerlaw spectrum. Via the inversion of the test, a confidence interval
for the energy cutoff can be constructed. However, special care must be taken
in the analysis of public data from the HGPS because the number of available
spectral points per $\gamma$-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 is 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
$\gamma$-ray sources in regions with complex morphology.
This section contains a description of statistical methods to constrain energy
cutoffs in powerlaw spectra. The application of the likelihood ratio (LR) test
and the F-test for the presence of cutoffs in powerlaw spectra are discussed
in section 4.1. In section 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 powerlaw $\gamma$-ray energy spectrum is discussed in section 4.3.

###  4.1 Models and hypotheses tests

Figure 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\leq p\leq 0.5$ where $50\%$ of the
p-values are predicted. The other half of the p-values is, due to the
constraint to positive energy cutoffs, expected at $p=1$ (see also appendix
A). The upper panel shows the resulting p-value distribution in the range
$0\leq p\leq 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
clearly visible 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 normal
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 Sec. 4.2.

Consider a powerlaw model for the $\gamma$-ray flux $\phi(E)$ at energy $E$
with an exponential cutoff at $E_{\mathrm{cut,\;\gamma}}=1/\lambda_{\gamma}$.
The model can be parameterized by

$\phi(E)=\phi_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma}\exp(-(\lambda_{\gamma}E)^{\beta})$
(3)

where again $\Gamma$ is the spectral index and $\phi_{0}$ the flux
normalization at energy $E_{0}$. The parameter $E_{0}$ is assumed to be fixed
to $1$ TeV in the following. Also fixed is the parameter $\beta$ which
controls the shape of the exponential cutoff. Frequently discussed is, for
example, the case $\beta=1$ (see e.g. Ahnen et al. (2017)). The hadronic model
given by Eq. 1 reduces to Eq. 3 with $\beta=0.5$ and
$\lambda_{\gamma}=16\lambda_{h}$. The cutoff parameter $\lambda_{\gamma}$ is
constrained to $\lambda_{\gamma}>0$ because a $\gamma$-ray flux suppression
for energies above $E_{\mathrm{cut,\;\gamma}}$ is expected for Pevatrons. The
model given by Eq. 3 has three free parameters
$\boldsymbol{\theta_{1}}=(\phi_{0},\;\Gamma,\;\lambda_{\gamma})$. The question
is whether the powerlaw model with exponential cutoff (Eq. 3) gives a better
description of a dataset than a powerlaw model parameterized by

$\phi_{\mathrm{PL}}(E)=\phi_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma}\;\mathrm{.}$
(4)

This model has two free parameters
$\boldsymbol{\theta_{0}}=(\phi_{0},\;\Gamma)$ while $E_{0}$ is again fixed to
$1$ TeV. To select the best fitting model, the null hypothesis
$\mathrm{H_{0}}$, i.e. the absence of an exponential cutoff
($\lambda_{\gamma}=0$), can be tested against the physically constrained
alternative hypothesis $\mathrm{H_{1}}$, $\lambda_{\gamma}>0$, by means of a
hypothesis test at a given confidence level $\mathrm{CL}$.
Let in the following $\phi_{i}$, $i=0,\dots,N-1$, denote binned $\gamma$-ray
flux measurements with bin-centers at energies $E_{i}$. All $N$ measurements
of $\phi_{i}$ are assumed to be independent. In practice, the independence of
the $\phi_{i}$ can be expected when the energy bin-width for the calculation
of $\phi_{i}$ is much larger than the instrumental energy resolution. Let
further $\sigma_{i}$ be the measurement errors corresponding to the
$\gamma$-ray fluxes $\phi_{i}$. The measurements $\phi_{i}$ are, in a first
step, assumed to be normal distributed around a hypothetical true value
$\hat{\phi}_{i}$. A model $\phi(E_{i}|\boldsymbol{\theta})$ with parameters
$\boldsymbol{\theta}$ predicts the true flux $\hat{\phi}_{i}$. Two frequently
used tests for a powerlaw hypothesis (Eq. 4) against a model with exponential
cutoff (Eq. 3) are presented in the following. Afterwards, 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 $\boldsymbol{\theta}$ of a given model factorizes
to
$L(\boldsymbol{\theta})=\prod_{i=0}^{N-1}\frac{1}{\sqrt{2\pi\sigma_{i}^{2}}}\exp\left(-\frac{(\phi_{i}-\phi(E_{i}|\boldsymbol{\theta}))^{2}}{2\sigma_{i}^{2}}\right)$.
Let
$\boldsymbol{\hat{\theta}_{k}}=\mathrm{argmax}_{H_{k}}L(\boldsymbol{\theta})$
be the maximum likelihood estimates for the model parameters under the
hypotheses $\mathrm{H_{k}}$ with $k\in\\{0,1\\}$. Consider the LR-test with
the statistic
$-2\ln\Lambda=-2\ln(L(\boldsymbol{\hat{\theta}_{0}})/L(\boldsymbol{\hat{\theta}_{1}}))=\mathrm{RSS}(\boldsymbol{\hat{\theta}_{0}})-\mathrm{RSS}(\boldsymbol{\hat{\theta}_{1}})$
where

$\mathrm{RSS}(\boldsymbol{\hat{\theta}_{k}})=\sum_{i=0}^{N-1}\left(\frac{\phi_{i}-\phi(E_{i}|\boldsymbol{\hat{\theta}_{k}})}{\sigma_{i}}\right)^{2}$
(5)

for the hypotheses $\mathrm{H_{k}}$. The null hypothesis is discarded at
confidence level $\mathrm{CL}$ if
$-2\ln\Lambda>\Lambda_{\mathrm{crit}}(\mathrm{CL})$. The critical value
$\Lambda_{\mathrm{crit}}$ can be calculated using the results from Chernoff
(1954) because the two models given by Eqs. 4 and 3 are nested, i.e. Eq. 3
reduces to Eq. 4 when $\lambda_{\gamma}\to 0$. When
$F_{\mathrm{\chi^{2}_{1}}}^{-1}$ denotes the inverse cumulative density
function of a $\chi^{2}$-distributed random variable with one degree of
freedom, the critical value is given by
$\Lambda_{\mathrm{crit}}=F_{\mathrm{\chi^{2}_{1}}}^{-1}(2\mathrm{CL}-1)$, see
appendix A. The test will have the expected false positive error rate
$1-\mathrm{CL}$ when the normal error model holds and the sample size $N$ is
large such that $\phi(E_{i}|\boldsymbol{\hat{\theta}_{k}})\to\hat{\phi}_{i}$
when $\mathrm{H_{0}}$ is true.

####  4.1.2 F-test

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

$F=(N-3)\frac{\mathrm{RSS}(\boldsymbol{\hat{\theta}_{0}})-\mathrm{RSS}(\boldsymbol{\hat{\theta}_{1}})}{\mathrm{RSS}(\boldsymbol{\hat{\theta}_{1}})}\;\mathrm{.}$
(6)

The null hypothesis $\lambda_{\gamma}=0$ is discarded in favor of the
alternative hypothesis $\lambda_{\gamma}>0$ at confidence level CL when
$F>F_{\mathrm{crit}}(\mathrm{CL})$. Similar to the LR-test, the critical value
is given by

$F_{\mathrm{crit}}=F_{1,N-3}^{-1}(2\mathrm{CL}-1)$ (7)

where $F^{-1}_{1,N-3}$ 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 non-linear models, such as in the considered
case, the false positive error rate can only be expected to be asymptotically
(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
$\gamma$-ray source HESS J1800-240. A fit of the data to a powerlaw model (Eq.
4) results in the best fit parameters
$\boldsymbol{\theta_{\mathrm{s}}}=(\phi_{0},\;\Gamma)=((4.3\pm 0.1)\cdot
10^{-13}\;\mathrm{TeV^{-1}cm^{-2}s^{-1}},2.44\pm 0.02)$. The robustness of the
LR- and F-test constructed above is tested in a Monte Carlo simulation of true
null models where random spectra are created based on the best fit model
$\phi(E|\boldsymbol{\theta_{\mathrm{s}}})$. For each simulated spectrum, all
five spectral points are sampled from a Student t-distribution
$t_{5}(\mathrm{location},\;\mathrm{scale})$ with location parameter
$\phi(E_{i}|\boldsymbol{\theta_{\mathrm{s}}})$, scale parameter $\sigma_{i}$
and $5$ degrees of freedom. For comparison, the same simulation is repeated by
sampling from a Normal distribution with mean
$\phi(E_{i}|\boldsymbol{\theta_{\mathrm{s}}})$ and standard deviation
$\sigma_{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.
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
model. The shown p-values refer to tests of a powerlaw model against a model
with exponential cutoff in the $\gamma$-ray spectrum (Eq. 3 with $\beta=1$).
Similar simulations with unchanged conclusions were performed when a powerlaw
model is tested against a model where the $\gamma$-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 $\gamma$-ray sources that were selected for the
analysis (see Sec. 3) confirm the observation that the LR-test is not robust
against violations of the normal error model in the HGPS catalog.
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}=\xi^{i}E_{\mathrm{min}}$ where
$\xi=2.4$ and $E_{\mathrm{min}}=0.4$ TeV are the median values of the binning
factor and the minimal energy bin for the selected HGPS data (see Tab. 1). At
each energy, the flux is sampled from
$\phi_{0}(E/E_{0})^{-\Gamma}\exp(-E/E_{\mathrm{cut,\;\gamma}})$ where
$\phi_{0}=3.51\cdot 10^{-12}\;\mathrm{cm^{-2}s^{-1}TeV^{-1}}$ and
$\Gamma=2.11$ are the median flux normalization and powerlaw index of the
selected HGPS data (see Tab. 1). The flux normalization energy $E_{0}=1$ TeV
is fixed and the cutoff energy $E_{\mathrm{cut,\;\gamma}}$ 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 the F-test. It is clearly seen that there is no true energy cutoff
where the F-test has more power than the LR-test. However, for example for a
true energy cutoff of $100$ TeV, the LR-test has $10\%$ more power than the
F-test.
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 usage 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, e.g. at a significance level of $5\sigma$. If the
powerlaw hypothesis is not discarded, a lower limit on the energy cutoff
$E_{\mathrm{cut}}$ is requested at reduced confidence level (e.g. $90\%$) to
constrain theoretical models. Depending on the alternative hypothesis,
$E_{\mathrm{cut}}$ can refer to $E_{\mathrm{cut,\;\gamma}}$ (Eq. 3) or to
$E_{\mathrm{cut,\;h}}$ (Eq. 1). A lower limit on the energy cutoff corresponds
to a one-sided confidence interval. Let $E_{\mathrm{cut}}^{\mathrm{MLE}}$ be
the maximum likelihood estimate of the energy cutoff as determined in a least-
square fit. Consider the interval

$I=(E_{\mathrm{cut}}^{-},\infty)$ (8)

where $E_{\mathrm{cut}}^{-}$ is the solution to the equation

$F(E)=F_{\mathrm{crit}}$ (9)

in the interval $0<E<E_{\mathrm{cut}}^{\mathrm{MLE}}$. Here,

$F(E)=(N-3)\frac{\mathrm{RSS}(\boldsymbol{\hat{\theta}_{0}^{*}}|E)-\mathrm{RSS}(\boldsymbol{\hat{\theta}_{1}})}{\mathrm{RSS}(\boldsymbol{\hat{\theta}_{1}})}$
(10)

and $\mathrm{RSS}$ is given by Eq. 5. The parameters
$\boldsymbol{\hat{\theta}_{0}^{*}}$ are the maximum likelihood estimators of
$\phi_{0}$ and $\Gamma$ for a powerlaw model with fixed energy cutoff at $E$.
Similarly, $\boldsymbol{\hat{\theta}_{1}}$ is the maximum likelihood estimator
for the parameters of a powerlaw model with variable exponential energy
cutoff, again with the constraint that $E_{\mathrm{cut}}>0$. The critical
value $F_{\mathrm{crit}}=F_{\mathrm{crit}}(\mathrm{CL})$ is given by Eq. 7.
Equation 10 is the result of an F-test of the hypothesis
$\mathrm{H_{0}}:\;E_{\mathrm{cut}}=E$ against the alternative hypothesis
$\mathrm{H_{1}}:\;E_{\mathrm{cut}}\neq E\;\mathrm{and}\;E_{\mathrm{cut}}>0$.
The existence of $E_{\mathrm{cut}}^{-}$, defined in Eq. 8, is shown in
appendix B. In the following, it is assumed that $E_{\mathrm{cut}}^{-}$ is the
unique solution of Eq. 10 in the energy interval
$0<E<E_{\mathrm{cut}}^{\mathrm{MLE}}$ (see also the discussion in appendix B).
It is shown in appendix C that the interval $I$ is expected to overcover 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 $\gamma$-ray energy spectra with energy cutoff shape
$\beta=1$ were simulated following the respective discussion in Sec. 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_{\mathrm{cut}}^{-}$ defined in Eq. 8
has very good coverage properties when the true cutoff of the $\gamma$-ray
spectrum is small. Otherwise it is conservative. For typical source and
instrumental parameters representative for the HESS galactic plane survey,
good coverage can be expected when the true cutoff of the $\gamma$-ray energy
spectrum is below $50$ TeV.

Figure 3: Shown in blue is the simulated coverage of the interval $I$ given by
Eq. 8 at $90\%$ CL. 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, for comparison, shown in yellow.

###  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_{\mathrm{cut}}^{\mathrm{MLE}}$ of the energy cutoff.
Systematic effects on the energy cutoff become dominant when the systematic
variation of $E_{\mathrm{cut}}^{\mathrm{MLE}}$ is not covered by the
confidence interval $I$. For data from Cherenkov telescope arrays like HESS
and VERITAS, reasons for systematic variations of
$E_{\mathrm{cut}}^{\mathrm{MLE}}$ can be, for example, 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 powerlaw model
are estimated to be $30\%$ on the flux normalization and $0.2$ on the powerlaw
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_{\mathrm{cut}}^{\mathrm{MLE}}$ in the HGPS.

##  5 Analysis of HGPS data

The analysis of the HGPS data, after the quality selection discussed in Sec.
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 $\gamma$-ray sources.
The HGPS catalog contains best fitting powerlaw 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. The best least square fit parameters are, however, in
reasonable agreement with the cataloged best fit parameters. The maximum
absolute difference between the powerlaw index obtained in this analysis and
the powerlaw index given the HGPS catalog is $1.8\sigma$ for HESS J1708-443
where $\sigma$ denotes the HGPS catalog index error. The median absolute
difference between the fitted powerlaw indices is $0.6\sigma$. The
corresponding largest absolute difference between the flux normalizations is
$1.2\sigma$ for HESS J1026-582. The median absolute flux difference for the
dataset is $0.1\sigma$.
Table LABEL:analysis_result shows the results of a search for spectral cutoffs
for the $25$ selected $\gamma$-ray sources in the HGPS. The best fit
exponential cutoff models (Eq. 3 with $\beta=1$ and Eq. 1) have a p-value
larger than $1\%$ in a $\chi^{2}$ test for all considered spectra. All
calculated solutions to Eq. 10 are unique. No energy cutoff is detected with a
statistical significance exceeding $5\sigma$ 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 powerlaw
index is in the range between $2.0$ and $2.5$. Figure 1 shows the spectral
data together with the best fit powerlaw 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_{\mathrm{cut}}^{\mathrm{MLE}}$ is not public. However, in case of the
source HESS J1908+063, an at least partially independent dataset is available.
This dataset is discussed in the next section.

Table 2: Results of the energy cutoff analysis for $25$ unidentified $\gamma$-ray sources in the HGPS catalog: The source name refers to the HGPS (Abdalla et al., 2018). $E_{\mathrm{cut,\;\gamma}}$ is the lower limit on the exponential energy cutoff of the $\gamma$-ray spectrum (Eq. 3 with $\beta=1$). $\mathrm{S_{\gamma}}$ is the significance of an F-test for a powerlaw model against an exponential cutoff in the $\gamma$-ray spectrum (again modeled by Eq. 3 with $\beta=1$). $E_{\mathrm{cut,\;h}}$ is the upper limit on the energy cutoff of a hadron population in a purely hadronic model (Eq. 1). $E_{\mathrm{max}}$ is the maximum bin-center energy of the spectral points for the source used in the analysis. $\Gamma$ is the powerlaw index of the best least square fit to the public HGPS data. Errors are statistical only. N is the number of spectral points for the source used in the analysis. The p-values refers to a $\chi^{2}$ goodness of fit test for the powerlaw model. The HESS livetime used for the spectral analysis of the source in the HGPS is indicated by T. All upper limits are at $90\%$ CL. The results are ordered descending in $E_{\mathrm{cut,\;h}}$. Source |  $E_{\mathrm{cut,\;\gamma}}$ (TeV) | $\mathrm{S_{\gamma}}$ |  $E_{\mathrm{cut,\;h}}$ (TeV) |  $E_{\mathrm{max}}$ (TeV) | $\Gamma$ | N | p-value | T (h)
---|---|---|---|---|---|---|---|---
HESS J1800-240 | 34.7 | 0.0 | 510.1 | 12 | 2.44$\pm$0.03 | 5 | 9.6e-01 | 10
HESS J1641-463 | 21.3 | 0.1 | 184.9 | 13 | 2.38$\pm$0.03 | 5 | 9.8e-01 | 27
HESS J1908+063 | 31.1 | 1.0 | 154.0 | 30 | 2.22$\pm$0.04 | 6 | 8.2e-01 | 12
HESS J1852-000 | 29.2 | 0.0 | 119.7 | 30 | 2.12$\pm$0.06 | 6 | 8.0e-01 | 23
HESS J1634-472 | 20.3 | 0.5 | 107.8 | 30 | 2.30$\pm$0.05 | 6 | 2.8e-01 | 14
HESS J1828-099 | 14.5 | 0.0 | 64.6 | 12 | 2.19$\pm$0.07 | 5 | 8.3e-01 | 20
HESS J1023-575 | 13.5 | 2.0 | 35.3 | 32 | 2.38$\pm$0.07 | 6 | 2.1e-01 | 23
HESS J1858+020 | 10.6 | 0.3 | 33.5 | 13 | 2.30$\pm$0.04 | 5 | 9.5e-01 | 22
HESS J1841-055 | 8.4 | 2.7 | 29.9 | 12 | 2.42$\pm$0.07 | 5 | 9.8e-03 | 14
HESS J1503-582 | 8.3 | 0.0 | 28.4 | 14 | 2.7$\pm$0.1 | 5 | 3.4e-01 | 18
HESS J1507-622 | 10.6 | 2.4 | 26.3 | 15 | 2.10$\pm$0.07 | 5 | 4.5e-01 | 13
HESS J1646-458 | 7.4 | 0.0 | 22.7 | 13 | 2.6$\pm$0.1 | 5 | 7.1e-01 | 10
HESS J1457-593 | 7.7 | 0.0 | 18.1 | 15 | 2.6$\pm$0.1 | 5 | 7.3e-01 | 4
HESS J1843-033 | 9.9 | 2.8 | 16.6 | 30 | 2.2$\pm$0.1 | 6 | 1.7e-03 | 17
HESS J1018-589 B | 8.6 | 2.3 | 14.7 | 15 | 2.1$\pm$0.1 | 5 | 6.3e-01 | 21
HESS J1632-478 | 5.3 | 0.1 | 10.0 | 13 | 2.46$\pm$0.06 | 5 | 4.0e-01 | 15
HESS J1741-302 | 4.5 | 0.0 | 9.2 | 7 | 2.2$\pm$0.2 | 6 | 5.3e-01 | 145
HESS J1702-420 | 5.2 | 0.2 | 9.1 | 13 | 2.1$\pm$0.1 | 5 | 1.7e-01 | 4
HESS J1809-193 | 3.4 | 1.0 | 5.0 | 12 | 2.3$\pm$0.1 | 5 | 1.6e-01 | 10
HESS J1616-508 | 3.2 | 1.1 | 3.8 | 13 | 2.3$\pm$0.1 | 5 | 1.1e-02 | 4
HESS J1808-204 | 2.1 | 0.0 | 2.7 | 12 | 2.1$\pm$0.2 | 5 | 1.3e-01 | 37
HESS J1026-582 | 5.0 | 2.0 | 2.3 | 32 | 2.0$\pm$0.2 | 5 | 7.0e-03 | 22
HESS J1804-216 | 1.3 | 1.1 | 1.8 | 12 | 2.7$\pm$0.1 | 5 | 2.6e-04 | 19
HESS J1708-443 | 2.2 | 1.2 | 1.3 | 13 | 2.0$\pm$0.2 | 5 | 4.2e-02 | 8
HESS J1746-285 | 1.2 | 0.6 | 0.7 | 6 | 2.0$\pm$0.2 | 7 | 4.0e-01 | 259
Figure 4: Spectral $\gamma$-ray data for the source MGRO J1908+06. Shown as
green circles is the HESS data from Aharonian et al. (2009). Shown in cyan
open circles is the spectral data from the HGPS. Blue spectral points are from
VERITAS (Aliu et al., 2014). VERITAS data marked with a square is 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 powerlaw model. Yellow and red lines are best fits to a powerlaw model with
cutoff at the $90\%$ CL lower limits $E_{\mathrm{cut,\;\gamma}}=29.5$ TeV
(yellow) and $E_{\mathrm{cut,\;h}}=141.5$ TeV (red).

##  6 Analysis of MGRO J1908+06 data

The $\gamma$-ray source MGRO J1908+06 is associated with HESS J1908+063. This
$\gamma$-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 section 5. Compared to the HGPS catalog, additional public
data for this $\gamma$-ray source is available. This data, acquired with
different experiments, is described and analyzed in this section.

###  6.1 Data for MGRO J1908+06

Milagro detected the $\gamma$-ray source at a median energy of $20$ TeV with a
flux of $(8.8\pm 2.4)\cdot 10^{-15}\;\mathrm{TeV^{-1}cm^{-2}s^{-1}}$ (Abdo et
al., 2007). Based on an analysis of $27$h data acquired with the HESS
experiment, spectral data is published in Aharonian et al. (2009). 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 the
livetime is larger. VERITAS also observed MGRO J1908+063 and collected $62$h
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 dataset are compatible with each other
(Aliu et al., 2014). Two lower limits on the energy cutoff
$E_{\mathrm{cut,\;\gamma}}$ of the $\gamma$-ray spectrum were derived before.
A $90\%$ CL lower limit on $E_{\mathrm{cut,\;\gamma}}$ 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_{\mathrm{cut,\;\gamma}}=19.1$ TeV (Aharonian et al., 2009).
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 section 3.2. The HESSd1 dataset passes this selection.
However, the public spectral data for MGRO J1908+063 from VERITAS is binned
with a factor $\xi=1.4$ while the energy resolution is similar to 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 powerlaw spectrum with $\phi_{0}=(4.2\pm 0.2)\cdot
10^{-12}\;\mathrm{cm^{-2}s^{-1}TeV^{-1}}$ and $\Gamma=2.17\pm 0.04$. The
powerlaw index is compatible with the index $\Gamma=2.26\pm 0.06$ cataloged in
the HGPS while the flux normalization in the HGPS, $(1.05\pm 0.09)\cdot
10^{-11}\;\mathrm{cm^{-2}s^{-1}TeV^{-1}}$, is incompatible. The fit to a
powerlaw model has a p-value of $0.91$ in a $\chi^{2}$ goodness of fit test.
The lower $90\%$ CL limits on the spectral energy cutoffs derived from this
dataset are $E_{\mathrm{cut,\;\gamma}}=29.5$ TeV and
$E_{\mathrm{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 in regard to the derived lower limits on spectral
energy cutoffs. Data and best fitting models are shown in Fig. 4.
Compared to the HESSd2 dataset with $6$ spectral points, the combined data
from HESSd1, VERITAS and Milagro considered in this section contains with $15$
spectral points many more flux measurements. This allows to search for a
spectral break in the TeV energy range. Multiple F-tests are performed where
the fit to the data of a powerlaw model ($H_{0}:\;\Delta\Gamma=0$) and a
broken powerlaw ($H_{1}:\;\Delta\Gamma>0$, see Eq. 2) are compared. The
comparison cannot be performed in a single test because the break energy
$E_{\mathrm{break}}$ is not defined under the null hypothesis. Let $E_{i}$
again be the energies of the $N$ available spectral flux data points. The
energy range $E_{2}=0.75$ TeV to $E_{N-2}=14.8$ TeV is scanned for a break
energy $E_{\mathrm{break}}$ with a logarithmic binning factor of $1.05$. The
F-test for the comparison of the fit quality of the powerlaw model and the
broken powerlaw 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 $\Delta\Gamma$ as a function of the assumed break
energy $E_{\mathrm{break}}$. The result is shown in Fig. 5. An index change
$\Delta\Gamma>0.5$ is ruled out at $90\%$ CL for energies between $1$ TeV and
$10$ TeV.
The largest local significance obtained in the $61$ tests of the energy break
scan is $p_{\mathrm{local}}=0.14$ with a F-test statistic $c=1.23$. The local
significance must be transformed into a global significance
$p_{\mathrm{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. This correction reads in the case
considered here

$p_{\mathrm{global}}\leq
p_{\mathrm{local}}+\left(\frac{N-3+c}{N-3+c_{0}}\right)^{-N/2+2}E[N_{c_{0}}]\;\mathrm{.}$
(11)

$E[N_{c_{0}}]$ is the expected number of upcrossings over the level $c_{0}$ of
the process of F-test statistics during the energy cutoff scan when the null
hypothesis is true. Equation 11 allows to extrapolate the global p-value
correction from a low test statistic level $c_{0}$, where a small number of
simulations is required for the estimation of $E[N_{c_{0}}]$, to a larger
value $c$ of the test statistic. Following Algeri et al. (2016b),
$E[N_{c_{0}}]$ is estimated in a parametric bootstrap simulation at
$c_{0}=0.3$ to be $E[N_{c_{0}}]=0.52\pm 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_{\mathrm{global}}=0.5$. The powerlaw hypothesis is therefore not discarded
in favor of a broken powerlaw.

Figure 5: Upper limit on a break in the powerlaw spectrum by $\Delta\Gamma$ as
a function of the break energy. Shown as red line is the $90\%$ CL upper limit
on $\Delta\Gamma$ derived from the combined data discussed in Sec. 6.

##  7 Discussion

The analysis of HGPS data resulted in five $\gamma$-ray sources for which the
energy spectrum is compatible with a powerlaw and, in a hadronic scenario, the
lower limit on the energy cutoff of the accelerated particle population is
larger than $100$ TeV. These five $\gamma$-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 $\gamma$-ray
morphology in Aharonian et al. (2008), based on a HESS dataset of $42$h
livetime. Four $\gamma$-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, spatially coincident with molecular
clouds.
Based on a smaller dataset of $10$h, the spectral HGPS analysis detects only
one $\gamma$-ray source in the region. The lower limit on the energy cutoff
for accelerated particles derived above is $E_{\mathrm{cut,\;h}}=510$ TeV at
$90\%$ CL. This constraint is comparable to the lower limit of
$E_{\mathrm{cut,\;h}}=600$ TeV at $90\%$ CL on the energy cutoff of particles
in the vicinity of the Galactic centre in a hadronic scenario (Abramowski et
al., 2016).
A physics case for the association of this $\gamma$-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 $\gamma$-ray source is discussed in Abramowski et al. (2014) where the TeV
spectrum, extracted from $72$h HESS data, is found to be compatible with a
powerlaw model. Additionally, a lower limit on the energy cutoff at
$E_{\mathrm{cut,\;h}}=100$ TeV is derived at $99\%$ CL in a hadronic model.
The lower limit of $E_{\mathrm{cut,\;h}}=185$ TeV at $90\%$ CL on the energy
cutoff derived in this work from $27$h HGPS data is compatible with the result
in Abramowski et al. (2014) but not more constraining.
A detailed hadronic model for the $\gamma$-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 $\gamma$-rays are produced via the decay of neutral
pions. The SNR G338.3-0.0 itself is associated with the $\gamma$-ray source
HESS J1640-465.
Given the presence of a nearby SNR and a molecular cloud, a hadronic scenario
for this $\gamma$-ray source is plausible. Observations with instruments which
have an improved sensitivity at few $10$ TeV to few $100$ TeV or a neutrino
detection are necessary to confirm a Pevatron scenario.

###  7.3 MGRO J1908+06

The $\gamma$-ray source MGRO J1908+06 was discovered by Milagro (Abdo et al.,
2007) and later confirmed by HESS (Aharonian et al., 2009). 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 whether the entire emission can originate from the PWN or
an additional $\gamma$-ray source is present. As candidate for an additional
$\gamma$-ray source, the SNR G40.5-0.5 is discussed in Aliu et al. (2014).
The HESS spectrum is confirmed by the measurement with VERITAS (Aliu et al.,
2014). However, the spectrum measured by ARGO-YBJ is incompatible with the
spectra measured by HESS and VERITAS (Bartoli et al., 2012). HAWC detects the
$\gamma$-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 $\gamma$-ray source is discussed as a possible
source of high energy neutrinos (Halzen et al., 2016). Aartsen et al. (2017)
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 $\gamma$-ray emission of
MGRO J1908+06 are currently not unambiguously measured. Also 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 gamma-ray spectrum derived
above, $E_{\mathrm{cut,\;\gamma}}\approx 30$ TeV, is a factor of almost $2$
more constraining than previous results discussed in Aharonian et al. (2009)
and Aliu et al. (2014). The lower limit $E_{\mathrm{cut,\;\gamma}}\approx 30$
TeV is derived from HGPS data ($E_{\mathrm{cut,\;\gamma}}=31.1$ TeV at $90\%$
CL) and confirmed in an analysis of combined data from HESS, VERITAS and
Milagro ($E_{\mathrm{cut,\;\gamma}}=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 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. (2009). The
powerlaw index of the electron spectrum measured on earth changes within
measured errors by $\Delta\Gamma_{e}=1$ at $E_{\mathrm{break}}=1$ TeV. This
energy break can be interpreted as cooling break (Recchia et al., 2019). When
electrons generate the $\gamma$-ray signal detected towards MGRO J1908+06 via
inverse Compton scattering in the Thomson regime, a break in the spectrum of
electrons by $\Delta\Gamma_{e}=1$ would translate into a change of the
$\gamma$-ray spectrum index by $\Delta\Gamma=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 $\gamma$-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 $\gamma$-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
has been analyzed before (Kosack et al., 2011), however, without derivation of
an energy spectrum. No other spectral data in the TeV $\gamma$-ray energy
range is publicly available. The analysis of the spectral cutoff for this
source is therefore not conclusive.
Different scenarios for the production of $\gamma$-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 unidentified source in the HGPS without further
discussion (Abdalla et al., 2018). The source is also detected in a previous
survey of the inner Galaxy with HESS (Aharonian et al., 2006). Here, the SNR
G337.2+0.1 and a source of X-rays are found to be within $0.2^{\circ}$ 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 $\gamma$-ray
source HESS J1634-472.
No further public TeV $\gamma$-ray data is available for this source. The
interesting constraint on $E_{\mathrm{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 $\gamma$-ray sources are found in the HGPS catalog for which the maximal
energy of accelerated particles is, in a hadronic model, at least $100$ TeV.
For at least $3$ of these sources, a hadronic scenario for the $\gamma$-ray
emission, as result from the interaction of runaway CRs with a molecular
cloud, is plausible. One of the Pevatron candidates found in the HGPS is MGRO
J1908+06. The $\gamma$-ray spectrum of this source extends to at least
$E_{\mathrm{cut,\;\gamma}}=30$ TeV without indication of a cutoff. The lower
limit on $E_{\mathrm{cut,\;\gamma}}$ for this source is a factor of almost two
more constraining than previous results. A break $\Delta\Gamma>0.5$ for this
$\gamma$-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., 2017) or $\gamma$-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 $\gamma$-ray data at few $10$ TeV. An improved
sensitivity and higher energies will be accessible with the upcoming
$\gamma$-ray detectors CTA (The CTA Consortium et al., 2011) and LHAASO
(Sciascio, 2016). Finally, the development of new experimental techniques,
such as for TAIGA/HiSCORE (Tluczykont et al., 2014) will help to identify
galactic Pevatrons and solve the quest of the origin of cosmic rays.

###### Acknowledgements.

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

## References

  * Aartsen et al. (2017) Aartsen, M., Ackermann, M., Adams, J., et al. 2017, Journal of Instrumentation, 12, P03012
  * Aartsen et al. (2017) Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017, ApJ, 835, 151
  * Abdalla et al. (2018) Abdalla et al. 2018, A&A, 612
  * Abdo et al. (2009) Abdo, A. A., Allen, B. T., Aune, T., et al. 2009, ApJ, 700
  * Abdo et al. (2007) Abdo et al. 2007, ApJ letters, 664
  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, ApJ, 843
  * Abramowski et al. (2014) Abramowski, A. et al. 2014, ApJ, 794
  * Abramowski et al. (2016) Abramowski et al. 2016, Nature, 531
  * Acero et al. (2013) Acero, F., Ackermann, M., Ajello, M., et al. 2013, ApJ, 773
  * Ackermann et al. (2013) Ackermann, M. et al. 2013, Science, 339
  * Aharonian et al. (2009) Aharonian, F., Akhperjanian, A., Anton, G., et al. 2009, A&A, 508
  * Aharonian et al. (2006) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, ApJ, 636
  * Aharonian et al. (2006) Aharonian, F. et al. 2006, A&A, 457
  * Aharonian et al. (2008) Aharonian, F. et al. 2008, A&A, 481
  * Aharonian et al. (2009) Aharonian, F. et al. 2009, A&A, 499
  * Ahnen et al. (2017) Ahnen et al. 2017, MNRAS, 472
  * Algeri et al. (2016a) Algeri, S., Conrad, J., & van Dyk, D. A. 2016a, MNRAS: Letters, 458, L84
  * Algeri et al. (2016b) Algeri, S., van Dyk, D., Conrad, J., & Anderson, B. 2016b, Journal of Instrumentation, 11
  * Algeri & van Dyk (2017) Algeri, S. & van Dyk, D. A. 2017, arXiv preprint, 1701.06820v4
  * Aliu et al. (2014) Aliu, E. et al. 2014, ApJ, 787
  * Amenomori et al. (2019) Amenomori et al. 2019, Phys. Rev. Lett., 123
  * Baade & Zwicky (1934) Baade, W. & Zwicky, F. 1934, Proc. of the National Academy of Sciences, 20
  * Bamba et al. (2016) Bamba, A., Terada, Y., Hewitt, J., et al. 2016, ApJ, 818, 63
  * Bartoli et al. (2012) Bartoli et al. 2012, ApJ, 760
  * Chernoff (1954) Chernoff, H. 1954, Ann. Math. Statist., 25
  * DeYoung (2012) DeYoung, T. 2012, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 692, 72
  * Gabici & Aharonian (2018) Gabici, S. & Aharonian, F. 2018, A&A, 612
  * Gabici & Aharonian (2007) Gabici, S. & Aharonian, F. A. 2007, ApJ, 665
  * Gaensler & Slane (2006) Gaensler, B. M. & Slane, P. O. 2006, Annual review of Astronomy and Astrophysics, 44, 17
  * Gallant (1975) Gallant, A. R. 1975, The American Statistician, 29
  * Gross & Vitells (2010) Gross, E. & Vitells, O. 2010, The European Physical Journal C, 70
  * Halzen et al. (2016) Halzen, F., Kheirandish, A., & Niro, V. 2016, arXiv preprint, 1609.03072v1
  * Hillas (2005) Hillas, A. M. 2005, J. Phys. G: Nucl. Part. Phys., 31
  * Hinton & Hofmann (2009) Hinton, J. A. & Hofmann, W. 2009, Annual Review of Astronomy and Astrophysics, 47
  * Holder (2011) Holder, J. 2011, International Cosmic Ray Conference, 12
  * Jogler & Funk (2016) Jogler, T. & Funk, S. 2016, ApJ, 816, 100
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74
  * Kosack et al. (2011) Kosack, K., Chaves, R. C. G., & Acero, F. 2011, International Cosmic Ray Conference, 7
  * Piron et al. (2001) Piron, F. et al. 2001, A&A, 374
  * Recchia et al. (2019) Recchia, S., Gabici, S., Aharonian, F. A., & Vink, J. 2019, Phys. Rev. D, 99
  * Sciascio (2016) Sciascio, G. D. 2016, Nuclear and Particle Physics Proceedings, 279-281, 166 , proceedings of the 9th Cosmic Ray International Seminar
  * Tang et al. (2015) Tang, Y., Yang, C., Zhang, L., & Wang, J. 2015, ApJ, 812
  * The CTA Consortium et al. (2011) The CTA Consortium, Actis, M., Agnetta, G., et al. 2011, Experimental Astronomy, 32
  * Tluczykont et al. (2014) Tluczykont, M., Hampf, D., Horns, D., et al. 2014, Astroparticle Physics, 56, 42
  * Wilks (1935) Wilks, S. S. 1935, Ann. Math. Statist., 9
  * Zhang & Liu (2019) Zhang, X. & Liu, S. 2019, ApJ, 874

##  Appendix A Critical values

Consider, as in Sec. 4.1, the two nested hypotheses $H_{0}:\lambda=0$ and
$H_{1}:\lambda>0$ for the single parameter $\lambda$. The parameter set for
$\mathrm{H_{0}}$ is not in the interior of the parameter set of
$\mathrm{H_{1}}$. The asymptotic result for the LR test statistic
$-2\ln\Lambda$ under $H_{0}$ discussed in Wilks (1935) is therefore not
applicable. A modification is, however, discussed in Chernoff (1954) where it
is predicted that $-2\ln\Lambda$ is asymptotically distributed like an equal
mixture of two random variables when $\mathrm{H_{0}}$ is true. In the
considered case, one of the random variables is $\chi^{2}$-distributed with
one degree of freedom. The other random variable has a Dirac delta function
probability density function $\delta(0)$ (see also Algeri et al. (2016a)).
With this asymptotic distribution $f_{\mathrm{H0}}$ of the test statistic
$-2\ln\Lambda$ under $\mathrm{H_{0}}$, the critical value
$\Lambda_{\mathrm{crit}}$ of the LR test statistic can be inferred from
$\int_{-\infty}^{\Lambda_{\mathrm{crit}}}dx\;f_{\mathrm{H0}}(x)=\mathrm{CL}$.
Let $F_{\mathrm{\chi^{2}_{1}}}$ denote the cumulative density functions of a
$\chi^{2}$-distribution with one degree of freedom. Then
$\int_{-\infty}^{\Lambda_{\mathrm{crit}}}dx\;f_{\mathrm{H0}}(x)=1/2F_{\mathrm{\chi^{2}_{1}}}(\Lambda_{\mathrm{crit}})+1/2\int_{-\infty}^{\Lambda_{\mathrm{crit}}}dx\;\delta(0)=\mathrm{CL}$.
This leads to the equality
$\Lambda_{\mathrm{crit}}=F_{\mathrm{\chi^{2}_{1}}}^{-1}(2\mathrm{CL}-1)$ which
can be used to calculate the critical value $\Lambda_{\mathrm{crit}}$ of the
LR test statistic at a given CL.
A similar result follows for the critical value of the asymptotic F statistic
in Sec. 4.1.2. Note that the F statistic (Eq. 6) and the LR statistic
$-2\ln\Lambda$ are related via
$F=(N-3)(-2\ln\Lambda)/\mathrm{RSS(\boldsymbol{\hat{\theta}_{1}})}$ (see Sec.
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_{\mathrm{crit}}$ of
the F test can again be inferred from
$1/2F_{1,N-3}(\Lambda_{\mathrm{crit}})+1/2\int_{-\infty}^{\Lambda_{\mathrm{crit}}}dx\;\delta(0)=\mathrm{CL}$
to be $F_{\mathrm{crit}}=F^{-1}_{1,N-3}(2\mathrm{CL}-1)$.

##  Appendix B Existence of the lower limit

Let $F(E)$ and $F_{\mathrm{crit}}$ be given by Eqs. 10 and 7. The lower limit
on the energy cutoff is defined as solution to the equation

$F(E)=F_{\mathrm{crit}}$ (12)

in the interval $(0,\;E_{\mathrm{cut}}^{\mathrm{MLE}})$. Here,
$E_{\mathrm{cut}}^{\mathrm{MLE}}$ is the maximum likelihood estimator for the
energy cutoff (see Sec. 4.3). It holds that $F_{\mathrm{crit}}>0$ and $F(E)$
is a continuous function. A solution to Eq. 12 does always exist in the given
interval because $F(E_{\mathrm{cut}}^{\mathrm{MLE}})=0$ and $F(E)\to\infty$
for $E\to 0$ by definition of $F(E)$.
However, the solution to Eq. 12 in the interval
$(0,\;E_{\mathrm{cut}}^{\mathrm{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. 12 are
found in the interval $(0,\;E_{\mathrm{cut}}^{\mathrm{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.
Let in the following $F$, $F(E)$ and $F_{\mathrm{crit}}$ be given by Eqs. 6,
10 and 7. It is assumed that only one solution $E_{\mathrm{cut}}^{-}$ for Eq.
12 in the interval $(0,E_{\mathrm{cut}}^{\mathrm{MLE}})$ exists. Similarly, it
is also assumed that at most one solution to Eq. 12 exists in the interval
$(E_{\mathrm{cut}}^{\mathrm{MLE}},\infty)$. Let $E_{\mathrm{cut}}^{+}$ be the
solution to Eq. 12 in $(E_{\mathrm{cut}}^{\mathrm{MLE}},\infty)$, if it
exists, and $E_{\mathrm{cut}}^{+}=\infty$ if no solution exists. Now consider
the interval $(E_{\mathrm{cut}}^{-},E_{\mathrm{cut}}^{+})$. This interval is
the acceptance interval $A=\\{E|F(E)\leq F_{\mathrm{crit}}\\}$ of an F-test
and, as such, defines a confidence interval for
$E_{\mathrm{cut}}^{\mathrm{MLE}}$. The confidence level of the F-test to which
$A$ is the acceptance interval depends, however, on $E$. This is a result of
the fit being constrained to $\lambda>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)\to F$ when $E\to\infty$ and $F(E)$ is the F-test statistic for a test of
a powerlaw hypothesis against a powerlaw model with exponential cutoff. This
means that $A$ is, in the limit $E\to\infty$, the acceptance interval for an
F-test constructed at confidence level CL. In other words, it holds that
$P(E_{\mathrm{true}}\in A)=\mathrm{CL}$. Additionally, by definition of the
interval $I$, it holds that $P(E_{\mathrm{true}}\in I)\geq
P(E_{\mathrm{true}}\in A)$. Together it follows $P(E_{\mathrm{true}}\in
I)\geq\mathrm{CL}$, i.e. the overcoverage of $I$ as a frequentist confidence
interval when the true energy cutoff is large.
Let the true energy cutoff now be small but positive such that the powerlaw
hypothesis is discarded at the given CL. The fit constraint $\lambda>0$ is in
this case irrelevant because the best fit $\lambda$ is always large. By
definition of $F_{\mathrm{crit}}$, it holds that $P(E_{\mathrm{true}}\in
A)=2\mathrm{CL}-1$. This is because without the constraint $\lambda>0$, $A$ is
the acceptance interval of an F-test that is constructed at confidence level
$2\mathrm{CL}-1$. Close to $E_{\mathrm{cut}}^{\mathrm{MLE}}$, $F(E)$ can be
quadratically approximated, see also Fig. 7. Using the nomenclature of Sec.
4.1, it holds that
$F(E)\approx(N-3)(E-E_{\mathrm{cut}}^{\mathrm{MLE}})^{2}/(\mathrm{RSS}(\boldsymbol{\hat{\theta}_{1}})\sigma^{2})$
where $\sigma^{2}$ is the inverse Fisher information of
$L(\boldsymbol{\theta_{1}})$ evaluated at $\boldsymbol{\hat{\theta}_{1}}$.
This means that the confidence interval $A$ will be symmetric around
$E_{\mathrm{cut}}^{\mathrm{MLE}}$ when $E$ is small such that the solution to
Eq. 12 is found in the range where the quadratic approximation holds. It
follows that
$P(E_{\mathrm{true}}\in(0,E_{\mathrm{cut}}^{-}))=P(E_{\mathrm{true}}\in(E_{\mathrm{cut}}^{+},\infty))=(1-(2\mathrm{CL}-1))/2=1-\mathrm{CL}$.
For the coverage probability, it follows that $P(E_{\mathrm{true}}\in
I)=P(E_{\mathrm{true}}\in
A)+P(E_{\mathrm{true}}\in(E_{\mathrm{cut}}^{+},\infty))=\mathrm{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.

Figure 6: 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_{\mathrm{cut,\;\gamma}}^{\mathrm{MLE}}$ for an exponential cutoff in the
$\gamma$-ray spectrum (Eq. 3 with $\beta=1$) is indicated with a green dashed
line. The blue and the magenta lines indicate the critical value
$F_{\mathrm{crit}}$ and the inferred lower limit on
$E_{\mathrm{cut,\;\gamma}}$ at $90\%$ CL. The magenta region is the excluded
range of cutoff energies. Figure 7: Same as Fig. 6 but based on the HGPS data
for the source HESS J1507-622. The powerlaw hypotheses is discarded at $90\%$
CL and the acceptance interval $A=\\{E|F(E)\leq F_{\mathrm{crit}}\\}$ is now a
finite interval roughly symmetric around
$E_{\mathrm{cut,\;\gamma}}^{\mathrm{MLE}}$. Shown in yellow is the quadratic
approximation of $F(E)$ around $E_{\mathrm{cut,\;\gamma}}^{\mathrm{MLE}}$.


Current affiliation: ]Quantum Matter Institute, University of British
Columbia, Vancouver, BC V6T 1Z4, Canada

# Direct Photon Scattering by Plasmons in BiTeI

A. C. Lee<EMAIL_ADDRESS>Department of Physics & Astronomy,
Rutgers University, Piscataway, NJ 08854, USA  S. Sarkar  Department of
Physics, Concordia University, 7141 Sherbrooke St. W, Montreal, Quebec  K. Du
Department of Physics & Astronomy, Rutgers University, Piscataway, NJ 08854,
USA  Rutgers Center for Emergent Materials, Rutgers University, Piscataway, NJ
08854, USA  H.-H. Kung  [  Department of Physics & Astronomy, Rutgers
University, Piscataway, NJ 08854, USA  C. J. Won  Laboratory for Pohang
Emergent Materials and Max Planck POSTECH Center for Complex Phase Materials,
Department of Physics, Pohang University of Science and Technology, Pohang
37673, Korea  K. Wang  Department of Physics & Astronomy, Rutgers University,
Piscataway, NJ 08854, USA  S.-W. Cheong  Department of Physics & Astronomy,
Rutgers University, Piscataway, NJ 08854, USA  Rutgers Center for Emergent
Materials, Rutgers University, Piscataway, NJ 08854, USA  Laboratory for
Pohang Emergent Materials and Max Planck POSTECH Center for Complex Phase
Materials, Department of Physics, Pohang University of Science and Technology,
Pohang 37673, Korea  S. Maiti<EMAIL_ADDRESS>Department of
Physics, Concordia University, 7141 Sherbrooke St. W, Montreal, Quebec  G.
Blumberg<EMAIL_ADDRESS>Department of Physics & Astronomy, Rutgers
University, Piscataway, NJ 08854, USA  National Institute of Chemical Physics
and Biophysics, 12618 Tallinn, Estonia

###### Abstract

We use polarization resolved Raman spectroscopy to show that for 3D giant
Rashba system the bulk plasmon collective mode directly couples to the Raman
response even in the long wavelength $\mathbf{q}\rightarrow 0$ limit although
the standard theory predicts that the plasmon spectral weight should scale as
the square of its quasi-momentum and hence be negligibly weak in the Raman
spectra. Such plasmon coupling to the Raman response at $\mathbf{q}\rightarrow
0$ arises for a polar system with spin-orbit coupling when the incoming photon
excitation is turned to a resonance with Rashba-split intermediates states
involved in the resonant Raman process. As an example, we identify special
features of BiTeI’s chiral band structure that enable the appearance of
plasmon mode in the Raman spectrum.

## Introduction:

Plasmons, the coherent oscillations of the electronic plasma relative to the
nuclei lattice, are longitudinal waves whose energy,
$\Omega_{pl}(\mathbf{q})$, is finite in 3D systems even in the long wavelength
limit where the quasi-momentum $\mathbf{q}\rightarrow 0$ [1, 2]. Although
Raman spectroscopy can probe optical modes, the spectral weight of scattering
from bulk plasmons scales as the square of its quasi-momentum and is therefore
vanishingly weak [3, 4].

This suppression of the spectral weight for Raman scattering from conventional
metals can be understood from the charge conservation perspective. In the
long-wavelength limit, there are no fluctuations in the total charge, without
which no dipole or quadrupole fields that could scatter photons are produced.
At finite $\mathbf{q}$, a longitudinal dipole field is produced and Raman
scattering being a two photon process couples to it as its square leading to
quadrupolar coupling. The selection rule that applies to the plasmon-photon
coupling dictates that the response would be picked up in the irreducible
representations (irreps) that contain the x,y,z fields with strengths
$q_{x}^{2}$, $q_{y}^{2}$ and $q_{z}^{2}$. Under conventional circumstances,
even if the scattering signal were boosted by a resonant Raman process, this
$q^{2}$ suppression would still persist [5]. This makes it challenging, if
not, impossible to directly observe scattering from or absorption by a
plasmon.

The numerous experimental Raman scattering studies [6, 7, 8, 5] that were able
to couple to the plasmon involved indirect help from extrinsic sources. E.g.,
in early studies, a technique that couples the plasmon to a boson of similar
energy at finite momentum to produce a polariton was commonly applied to
observe a plasmon in Raman spectra, particularly for low-carrier concentration
semiconductors [9, 10, 11, 12]. The application of a high magnetic field has
also been shown to produce observable by Raman collective charge-density
fluctuations [13, 14] such as in the two-dimensional fractional quantum Hall
systems [15, 16, 17, 18, 19].

In systems with inversion, there is no intrinsic vector that could produce a
dipole field and hence a finite $q$ is required for Raman coupling to light.
In prior experiments the visibility of the plasmon is enhanced either by a
grating-induced $\mathbf{q}$ that breaks the translational symmetry and
inversion or by borrowing the spectral weight from other bosons with which the
plasmon hybridises [9, 10, 11]. In polar systems with broken inversion,
however, one can expect non-$q$ sources of creating this dipolar field. It
then opens up the question of if the photon scattering from plasmons could
have a different weight from $q^{2}$ in such systems leading to a potential
direct observation of plasmon in Raman scattering.

Figure 1:  Basic properties of BiTeI. (a) The crystal structure. (b) Spin
texture of Rashba-spin split bands. The inner and outer Rashba bands have
opposite helicities for the same in-plane quasi-momentum direction. The
orbital Rashba effect results in slight z-canting of the spins. (c) Energy-
band dispersion relative to the Fermi energy obtained by DFT calculations
[20].

## Choice of material:

The observation above is the prime motivation to study the 3D giant Rashba
polar material, BiTeI [20]. The unit cell has Bi3+ atoms that are sandwiched
between layers of I- and Te2- atoms forming a tri-layer of atoms [space group
P3m1 (C3v)] stacked along the c-axis, Fig. 1(a), which gives it a non-
centrosymmetric nature and is the origin of the orbital Rashba effect which
intrinsically couples the charge and spin degrees of freedom without the need
of external gating or surface-induced symmetry breaking. The different layers
are van der Waals-bonded strong enough to result in 3D band structure. This
results in spin split helical electronic bands [21] with the spins also
showing a slight out-of-plane canting, Fig. 1(b). Although expected to be an
insulator, the I- vacancies that occur during the crystal growth cooling
process raise the chemical potential into the lower subbands [22] [labelled as
CB0 in Fig. 1(c)], resulting in a metallic character with low carrier
concentrations, $n_{c}$, such that the out-of-plane plasma frequency is
expected to be in 60-80 meV range [23]. The self-doping is such that chemical
potential remains below the Dirac-like point. This aspect introduces a gap
below the inter-sub-band excitations [24, 25] which will play a central role
in the data interpretation.

## Main finding:

In this Letter, we report the observation of a new spin mediated photon-
plasmon scattering channel in BiTeI which leads to a prominent mode around 74
meV in the spectrum which is below the inter-sub-band continuum of excitations
and is shown to be attributed to scattering in fully symmetric (A1) channel,
see Fig. 2. By studying the resonant Raman excitation profile (RREP), Fig. 3,
we infer that the scattering turns on by resonant excitation between spin-
split bands CB0 and CB1, see Fig. 1(c), but not between bands where there is
no mutual spin splitting like between CB0 and CB2. By studying different
samples from the same batch with varying carrier density we show that the
square of the mode’s energy scales linearly with the doped carriers
concentration $n_{c}$, thereby confirming the plasmonic nature of the mode. We
also track the spectral line-shape of the plasmon as a function of temperature
and show that its relaxation rate is consistent with the behavior of an in-gap
collective mode.

## Experimental details:

Single crystals of BiTeI were grown using the vertical Bridgman technique [26,
27]. Four-probe Hall effect measurements were performed using a Physical
Property Measurement System (PPMS) by Quantum Design to determine the $n_{c}$
of the samples; see 111See Supplementary Material at [URL will be inserted by
publisher] for the details of crystal growth and characterization, measurement
methods, analysis of incident angle dependent Raman data, Raman selection
rules, analysis of Raman data temperature dependence, which includes Refs.
[26, 27, 36, 37, 38, 39, 40, 41, 41, 42, 20] for full experimental details.
The samples used in this work are denoted as A, B, C, D, and E. The crystals
were exfoliated in a nitrogen-rich environment and then transferred into a
continuous He-flow optical cryostat. We use polarization-resolved Raman
scattering measurements in the quasi-back-scattering geometry. Scattering
polarization geometries, defined by the incident (collected) light
polarizations, ei(s), and light propagation directions relative to the crystal
structure, $\hat{k}_{i(s)}$, will be denoted as
$\hat{k}_{i}$(eies)$\hat{k}_{s}$, where $\hat{k}_{i}=-\hat{k}_{s}$. The linear
light polarization directions, X and Y, lie within the xy-plane and are
orthogonal to one another; due to the symmetries of $C_{3v}$, the choice of
in-plane linear light polarization directions does not affect the results. R
and L denote right and left circularly polarized light, respectively, such
that R(L) = X $\pm i$Y [28]. Finally, the energy of the incident (collected)
photons will be denoted as $\omega_{i(s)}$, respectively with
$\omega=\omega_{i}-\omega_{s}$ being the Raman shift. The Raman susceptibility
will be denoted as $\chi^{{}^{\prime\prime}}(\omega)$.

Figure 2: Secondary emission and Raman response of BiTeI. (a) Secondary
emission spectra taken from sample A in the z(RR)$\bar{z}$ (blue),
z(RL)$\bar{z}$ (green), z(XX)$\bar{z}$ (red), and z(YX)$\bar{z}$ (black)
scattering geometries using low-resolution settings. The vertical dotted lines
represent the bounds of the Rashba continuum, $\Omega_{\pm}$ [25]. (b)
Symmetry decomposed Raman spectra into $A_{1}$, $A_{2}$, and $E$ channels
[28]. No noticeable $E$-symmetry excitations where detected within this
frequency range. All measurements were taken at 15 K using $\omega_{i}$ = 1.83
eV excitation energy. The upper scale denotes absolute energy of emission
$\omega_{s}$ and the lower scale denotes Raman shift energy
$\omega=\omega_{i}-\omega_{s}$.

## Data and interpretation:

In Fig. 2(a), we show spectra of secondary emission from the sample A measured
in four in-plane polarization geometries. The broad centered at
$\omega_{s}=1.48$ eV energy feature (shaded in violet) is photo-luminescence
(PL) which appears in all polarization channels [25]. By subtracting the PL
and decomposing the remaining signal into symmetrized Raman response channels,
see [28] for details, several excitation features become evident, Fig. 2(b): A
sharp in-gap collective mode appears at 74 meV in the $A_{1}$ symmetry channel
that we denote as $A_{1}^{*}$ and two electronic continua of excitations; the
$A_{1}$ continuum that extends from just above the $A_{1}^{*}$ mode to 500 meV
and the $A_{2}$ Rashba continuum, which was studied in Ref. [25], that is an
order of magnitude stronger than the $A_{1}$ continuum. No appreciable signal
was collected in the E symmetry channel.

Figure 3: Resonance Raman excitation profile of the $A_{1}^{*}$ mode. (a) The
Raman response in fully symmetric symmetry channel from the sample A taken in
z(XX)$\bar{z}$ scattering polarization for $\omega_{i}=$1.55–2.60 eV
excitation energies. (b) The resonance Raman excitation profile (RREP) -
spectral intensity of the $A_{1}^{*}$ mode integrated over Raman shift and
plotted against $\omega_{i}$. The real part of the optical conductivity,
$\sigma_{1}(\omega_{i})$, is overlaid for reference.

We collected Raman spectra of the $A_{1}^{*}$ mode as a function of excitation
energy $\omega_{i}$ to track its spectral weight. In Fig. 3(a), we show the
Raman response of BiTeI using a set of laser excitation energies between 1.55
eV and 2.60 eV. In Fig. 3(b) we plot the resonance Raman excitation profile
(RREP) of the $A_{1}^{*}$ mode, $I(\omega_{i})$, where the real part of the
optical conductivity, $\sigma_{1}(\omega_{i})$, is also overlaid for
reference. The RREP of the $A_{1}^{*}$ mode peaks at $\omega_{i}=1.83$ eV and
is nearly undetectable away from this resonance excitation. Based on the RREP
and calculations of the electronic band structure of BiTeI, the resonance
enhancement of the mode is a result of intermediate state transitions between
the CB0 and CB1, see Fig. 1(c). The energy gap between the two bands can be
discerned from the broad maxima in the $\sigma_{1}(\omega_{i})$ spectra.
Although there is an additional broad maximum in $\sigma_{1}(\omega_{i})$ at
$\omega_{i}\sim 2.7$ eV, the $A_{1}^{*}$ mode’s RREP does not have a second
peak at this excitation energy. Since the CB1 bands are also Rashba spin-split
[20], this suggests that the physical mechanism that allows for the
observation of the mode involves the Rashba spin splitting of the intermediate
states as well.

Figure 4: Raman response of BiTeI as a function of temperature. The Raman
response for the sample A taken in z(XX)$\bar{z}$ scattering geometry at
temperatures between 10 to 300 K using $\omega_{i}=1.83$ eV resonant
excitation energy. (inset) The energy and HWHM values derived from Lorentzian
fits to the $A_{1}^{*}$ mode lineshape accounting for spectral resolution.

We also investigated the temperature dependence of the $A_{1}^{*}$ mode near
its excitation resonance. In Fig. 4, we show the measured Raman response for
10–300 K temperatures range. In the inset of Fig. 4, we plot the energy,
$\omega_{p}(T)$, and the half-width-at-half-maxima (HWHM), $\gamma(T)$, of the
mode as a function of temperature. The spectral lineshape monotonically red-
shifts and sharpens upon cooling. These results are consistent with the
picture of a $T$-dependent damped oscillator. However, the lineshape of the
$A_{1}^{*}$ mode is slightly asymmetrical, suggesting that a more
sophisticated model is required to understand damping of the in-gap collective
mode.

## Analysis and justification:

The selective presence of the mode in $A_{1}$ channel establishes the charge
nature of the excitation and thus rules out possible chiral spin modes
expected in Rashba metals, as they would appear in the $A_{2}$ symmetry
channel [29]. We repeated the same experimental procedure for samples B, C, D,
and E which had different concentration of I- vacancies. We found the carrier
densities, $n_{c}$, for all the samples using Hall measurements, Fig. 5(b),
and tracked the energy of the $A_{1}^{*}$ mode, Fig. 5(a). Upon plotting the
square of the mode’s energy as function of the carrier concentration $n_{c}$,
Fig. 5(c), we observe a straight line confirming its plasmonic nature [30].
The solid black line is the best fit to the data points with the constraint
$\Omega_{pl}^{2}(0)=0$. Within parabolic electron dispersion model with out-
of-plane effective mass $m_{c}^{*}$, the slope of this line is given by $4\pi
e^{2}\hbar^{2}/m_{c}^{*}\epsilon_{\infty}=1.2\times 10^{-16}$ [meV2-cm3]. The
derived value of $m_{c}^{*}=1.02m_{e}$, where $m_{e}$ is the mass of a free
electron, well agrees within the value acquired by the optical reflectivity
measurements [23]. We have assumed that for the c-axis direction
$\epsilon_{\infty}=12.2$, which was found from grazing-angle reflectivity
measurements [23].

Figure 5: The relationship between the square of the $A_{1}^{*}$ mode energy
and $n_{c}$. (a) Secondary emission spectra of the $A_{1}^{*}$ mode for
samples B, C, D, and E taken at room temperature. The collective mode is
fitted to a Lorentzian-like function [31] after background subtraction and
shifted vertically by a constant proportional to $n_{c}$ of the sample to aid
the reader’s eye. (b) Transverse resistivity data for the samples listed in
(a). The data has been normalized to account for sample dimensions. (c)
$\Omega_{pl}^{2}$ plotted against $n_{c}$. The solid black line is the line-
of-best-fit for the data points fixing the y-intercept to be zero.

There is another important aspect of the observed plasmon mode that requires
discussion. One usually thinks of the plasmon to exist above the continuum of
particle-hole excitations. From Fig. 2(b) we observe that the plasmon is below
the continuum of excitations. This apparent discrepancy is resolved once we
note that the continuum here is the chiral continuum of spin-flip excitations
[25]. The out-of-plane plasmon is still above the charge sector particle hole
continuum that is not visible in Raman spectrum due to $\mathbf{q}\rightarrow
0$. It is known that in 2D Rashba metals at finite $\mathbf{q}$ the plasmon
can enter and be damped by the spin-flip continuum [32]. It is also known that
that in a 3D Rashba metal, the out-of-plane plasmon does not couple to this
spin-flip continuum [33]. The novel feature of this plasmon is that it couples
to and is damped by the proximity to the chiral spin-flip continuum [31]. This
is evident from the weak scattering signal in the continuum region of the
$A_{1}$ channel, and also from the temperature evolution of the scattering
rate of the $A_{1}^{*}$ mode. What distinguishes BiTeI from a 3D Rashba metal
[which only models the in-plane spin-orbit coupling (SOC)] is the presence
$z$-canting of spins which is about 10% of the in-plane spin [34, 35]. We
attribute the sensitivity of the $A_{1}$ response to the spin-flip continuum,
and hence the coupling of the plasmon to the spin-flip continuum, to the
presence of the $z$-canting in the chiral band structure. The details of this
coupling are present in a concomitant work that explores the observations of
plasmons in SOC systems on general grounds [31]. Because the $z$-canting is
small, we would expect the plasmon to still be well resolved even if it
entered the spin-flip continuum.

From a group theoretical perspective, the in-plane plasmons couple to the $E$
representation and the out-of-plane plasmon to the $A_{1}$ representation of
the C3v point group, but both still couple with spectral weight $q^{2}$
without SOC. Since BiTeI is a polar system where the out-of-plane dipole
($\bf{z}$) and the in-plane quadrupole ($x^{2}+y^{2}$) excitation transform as
the same fully-symmetric A1 representation, they mix allowing the plasmon to
be detectable in $A_{1}$. The finite spectral weight (non-$\mathbf{q}^{2}$)
for the coupling is main discovery of this work and is attributed to the
$z$-canting of the spin states [31].

## Conclusions:

In this work, we have demonstrated that the bulk plasmon appears in the Raman
spectra of the 3D giant Rashba metal, BiTeI. The tell-tale signatures of
$\sqrt{n_{c}}$-dependence, $A_{1}$ symmetry, and temperature dependence of the
plasma frequency and linewidth are present. These results are supported by the
fact that the plasmon energy is below the Rashba continuum of single-particle
excitations, thus the mode is not overdamped. The excitation profile of the
plasmon demonstrates that observation of the collective mode is possible only
when the laser excitation energy is resonantly tuned to the CB0-CB1 transition
energy. This discovery sheds light on the important role of inversion breaking
in the manifestation of special effects in optical spectroscopy. Two important
questions arise from these results: (1) why do the initial and intermediate
states need to be Rashba spin-split and (2) why is resonance needed to observe
the plasmon? In a concomitant study [31], we delineate how the coupling
between spin-split Rashba bands in the resonant Raman effect could dominate
over the conventional $q^{2}$-coupling in Rashba systems.

We thank Bo Peng, Zuzhang Lin and Bartmeu Monserrat for discussions. The
spectroscopic work at Rutgers (A.C.L., H.-H.K. and G.B.) was supported by the
NSF under Grant No. DMR-2105001. The characterization work (K.D. and S.-W.C)
was supported by the DOE under Grant No. DOE: DE-FG02-07ER46382. The crystal
growth at Postech was supported by the National Research Foundation of Korea
(NRF) funded by the Ministry of Science and ICT, grant no. 2020M3H4A2084417.
S.S. and S.M. acknowledge support from the Natural Sciences and Engineering
Research Council of Canada (NSERC) Grant No. RGPIN-2019-05486. The work at
NICPB (G.B.) was supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme Grant
Agreement No. 885413.

## References

  * Pines [2018] D. Pines,  _Elementary Excitations in Solids Lectures on Protons, Electrons, and Plasmons_ (CRC Press, 2018).
  * Mahan [1981] G. Mahan,  _Many-Particle Physics_ (Plenum Publishers, 1981).
  * Pines [1966] D. Pines,  _The Theory of Quantum Liquids_ (Benjamin, New York, 1966).
  * Hayes and Loudon [2012] W. Hayes and R. Loudon,  _Scattering of Light by Crystals_ (Courier Corporation, 2012).
  * Das Sarma and Wang [1999] S. Das Sarma and D.-W. Wang, Resonant Raman Scattering by Elementary Electronic Excitations in Semiconductor Structures, Phys. Rev. Lett. 83, 816 (1999).
  * Pinczuk  _et al._ [1971] A. Pinczuk, L. Brillson, E. Burstein, and E. Anastassakis, Resonant light scattering by single-particle electronic excitations in $n-\mathrm{GaAs}$, Phys. Rev. Lett. 27, 317 (1971).
  * Pinczuk  _et al._ [1989] A. Pinczuk, S. Schmitt-Rink, G. Danan, J. P. Valladares, L. N. Pfeiffer, and K. W. West, Large exchange interactions in the electron gas of GaAs quantum wells, Phys. Rev. Lett. 63, 1633 (1989).
  * Goñi  _et al._ [1991] A. R. Goñi, A. Pinczuk, J. S. Weiner, J. M. Calleja, B. S. Dennis, L. N. Pfeiffer, and K. W. West, One-dimensional plasmon dispersion and dispersionless intersubband excitations in GaAs quantum wires, Phys. Rev. Lett. 67, 3298 (1991).
  * Varga [1965] B. B. Varga, Coupling of Plasmons to Polar Phonons in Degenerate Semiconductors, Phys. Rev. 137, A1896 (1965).
  * Mooradian and Wright [1966] A. Mooradian and G. B. Wright, Observation of the Interaction of Plasmons with Longitudinal Optical Phonons in GaAs, Phys. Rev. Lett. 16, 999 (1966).
  * Mooradian and McWhorter [1967] A. Mooradian and A. L. McWhorter, Polarization and Intensity of Raman Scattering from Plasmons and Phonons in Gallium Arsenide, Phys. Rev. Lett. 19, 849 (1967).
  * Abstreiter  _et al._ [1979] G. Abstreiter, R. Trommer, M. Cardona, and A. Pinczuk, Coupled plasmon-LO phonon modes and Lindhard-Mermin dielectric function of n-GaAs, Solid State Communications 30, 703 (1979).
  * Slusher  _et al._ [1967] R. E. Slusher, C. K. N. Patel, and P. A. Fleury, Inelastic light scattering from landau-level electrons in semiconductors, Phys. Rev. Lett. 18, 77 (1967).
  * Patel and Slusher [1968] C. K. N. Patel and R. E. Slusher, Light scattering from electron plasmas in a magnetic field, Phys. Rev. Lett. 21, 1563 (1968).
  * Girvin  _et al._ [1985] S. M. Girvin, A. H. MacDonald, and P. M. Platzman, Collective-Excitation Gap in the Fractional Quantum Hall Effect, Phys. Rev. Lett. 54, 581 (1985).
  * Girvin  _et al._ [1986] S. M. Girvin, A. H. MacDonald, and P. M. Platzman, Magneto-roton theory of collective excitations in the fractional quantum Hall effect, Phys. Rev. B 33, 2481 (1986).
  * Pinczuk  _et al._ [1993] A. Pinczuk, B. S. Dennis, L. N. Pfeiffer, and K. West, Observation of collective excitations in the fractional quantum Hall effect, Phys. Rev. Lett. 70, 3983 (1993).
  * Pinczuk  _et al._ [1994] A. Pinczuk, B. Dennis, L. Pfeiffer, K. West, and E. Burstein, Inelastic light scattering by the two-dimensional electron gas: fractional quantum Hall regime and beyond, Philosophical Magazine B 70, 429 (1994).
  * Eriksson  _et al._ [1999] M. A. Eriksson, A. Pinczuk, B. S. Dennis, S. H. Simon, L. N. Pfeiffer, and K. W. West, Collective Excitations in the Dilute 2D Electron System, Phys. Rev. Lett. 82, 2163 (1999).
  * Ishizaka  _et al._ [2011] K. Ishizaka, M. Bahramy, H. Murakawa, M. Sakano, T. Shimojima, T. Sonobe, K. Koizumi, S. Shin, H. Miyahara, A. Kimura, _et al._ , Giant Rashba-type spin splitting in bulk BiTeI, Nature materials 10, 521 (2011).
  * Park  _et al._ [2011] S. R. Park, C. H. Kim, J. Yu, J. H. Han, and C. Kim, Orbital-Angular-Momentum Based Origin of Rashba-Type Surface Band Splitting, Phys. Rev. Lett. 107, 156803 (2011).
  * Tomokiyo  _et al._ [1977] A. Tomokiyo, T. Okada, and S. Kawano, Phase Diagram of System (Bi2Te3)-(BiI3)and Crystal Structure of BiTeI, Japanese Journal of Applied Physics 16, 291 (1977).
  * Lee  _et al._ [2011] J. S. Lee, G. A. H. Schober, M. S. Bahramy, H. Murakawa, Y. Onose, R. Arita, N. Nagaosa, and Y. Tokura, Optical Response of Relativistic Electrons in the Polar BiTeI Semiconductor, Phys. Rev. Lett. 107, 117401 (2011).
  * Ye  _et al._ [2015] L. Ye, J. G. Checkelsky, F. Kagawa, and Y. Tokura, Transport signatures of Fermi surface topology change in BiTeI, Phys. Rev. B 91, 201104 (2015).
  * Lee  _et al._ [2022] A. C. Lee, B. Peng, K. Du, H.-H. Kung, B. Monserrat, S. W. Cheong, C. J. Won, and G. Blumberg, Chiral electronic excitations in the quasi-two-dimensional Rashba system BiTeI, Phys. Rev. B 105, L161105 (2022).
  * Sankar  _et al._ [2014] R. Sankar, I. Panneer Muthuselvam, C. J. Butler, S.-C. Liou, B. H. Chen, M.-W. Chu, W. L. Lee, M.-T. Lin, R. Jayavel, and F. C. Chou, Room temperature agglomeration for the growth of BiTeI single crystals with a giant Rashba effect, CrystEngComm 16, 8678 (2014).
  * Kanou and Sasagawa [2013] M. Kanou and T. Sasagawa, Crystal growth and electronic properties of a 3D Rashba material, BiTeI, with adjusted carrier concentrations, Journal of Physics: Condensed Matter 25, 135801 (2013).
  * Note [1] See Supplementary Material at [URL will be inserted by publisher] for the details of crystal growth and characterization, measurement methods, analysis of incident angle dependent Raman data, Raman selection rules, analysis of Raman data temperature dependence, which includes Refs. [26, 27, 36, 37, 38, 39, 40, 41, 41, 42, 20].
  * Kung  _et al._ [2017] H.-H. Kung, S. Maiti, X. Wang, S.-W. Cheong, D. L. Maslov, and G. Blumberg, Chiral Spin Mode on the Surface of a Topological Insulator, Phys. Rev. Lett. 119, 136802 (2017).
  * Cardona and Peter [2005] M. Cardona and Y. Y. Peter,  _Fundamentals of semiconductors_ , Vol. 619 (Springer, 2005) pp. 427–428.
  * Sarkar  _et al._ [2023] S. Sarkar, A. Lee, G. Blumberg, and S. Maiti, Spin-orbit interaction enabled electronic Raman scattering from charge collective modes, (2023), arXiv:2306.11240 [cond-mat.mes-hall] .
  * Pletyukhov and Gritsev [2006] M. Pletyukhov and V. Gritsev, Screening in the two-dimensional electron gas with spin-orbit coupling, Phys. Rev. B 74, 045307 (2006).
  * Maiti  _et al._ [2015] S. Maiti, V. Zyuzin, and D. L. Maslov, Collective modes in two- and three-dimensional electron systems with Rashba spin-orbit coupling, Phys. Rev. B 91, 035106 (2015).
  * Bawden  _et al._ [2015] L. Bawden, J. M. Riley, C. H. Kim, R. Sankar, E. J. Monkman, D. E. Shai, H. I. Wei, E. B. Lochocki, J. W. Wells, W. Meevasana, T. K. Kim, M. Hoesch, Y. Ohtsubo, P. Le Fèvre, C. J. Fennie, K. M. Shen, F. Chou, and P. D. C. King, Hierarchical spin-orbital polarization of a giant Rashba system, Science Advances 1, 10.1126/sciadv.1500495 (2015).
  * Maaß  _et al._ [2016] H. Maaß, H. Bentmann, C. Seibel, C. Tusche, S. V. Eremeev, T. R. Peixoto, O. E. Tereshchenko, K. A. Kokh, E. V. Chulkov, J. Kirschner, _et al._ , Spin-texture inversion in the giant Rashba semiconductor BiTeI, Nature communications 7, 11621 (2016).
  * Hurd [2012] C. Hurd,  _The Hall effect in metals and alloys_ (Springer Science & Business Media, 2012) pp. 8–10.
  * Goldstein  _et al._ [2002] H. Goldstein, C. Poole, and J. Safko, Classical mechanics (2002).
  * Loudon [1965] R. Loudon, Theory of the resonance Raman effect in crystals, Journal de Physique 26, 677 (1965).
  * Blumberg  _et al._ [1994] G. Blumberg, M. Klein, L. Börjesson, R. Liang, and W. Hardy, Investigation of the temperature dependence of electron and phonon Raman scattering in single crystal YBa2Cu3O6.952, Journal of superconductivity 7, 445 (1994).
  * Makhnev  _et al._ [2014] A. A. Makhnev, L. V. Nomerovannaya, T. V. Kuznetsova, O. E. Tereshchenko, and K. A. Kokh, Optical properties of BiTeI semiconductor with a strong Rashba spin-orbit interaction, Optics and Spectroscopy 117, 764 (2014).
  * Dresselhaus  _et al._ [2008] M. S. Dresselhaus, G. Dresselhaus, and A. Jorio, Applications of Group Theory to the Physics of Solids (2008).
  * Basov and Timusk [2005] D. N. Basov and T. Timusk, Electrodynamics of high-${T}_{c}$ superconductors, Rev. Mod. Phys. 77, 721 (2005).