An Unexpected Dip in the Solar Gamma-Ray Spectrum
The solar disk is a bright source of multi-GeV gamma rays, due to the interactions of hadronic cosmic rays with the solar atmosphere. However, the underlying production mechanism is not understood, except that its efficiency must be greatly enhanced by magnetic fields that redirect some cosmic rays from ingoing to outgoing before they interact. To elucidate the nature of this emission, we perform a new analysis of solar atmospheric gamma rays with 9 years of Fermi-LAT data, which spans nearly the full 11-year solar cycle. We detect significant gamma-ray emission from the solar disk from 1 GeV up to GeV. The overall gamma-ray spectrum is much harder () than the cosmic-ray spectrum (). We find a clear anticorrelation between the solar cycle phase and the gamma-ray flux between 1–10 GeV. Surprisingly, we observe a spectral dip between 30–50 GeV in an otherwise power-law spectrum. This was not predicted, is not understood, and may provide crucial clues to the gamma-ray emission mechanism. The flux above 100 GeV, which is brightest during the solar minimum, poses exciting opportunities for HAWC, LHAASO, IceCube, and KM3NeT.
The solar disk is a bright, hard-spectrum gamma-ray source due to its constant bombardment by hadronic cosmic rays, which interact inelastically with gas in the solar atmosphere (including deep under the optical photosphere) and produce gamma rays. If magnetic fields are ignored, the expected gamma-ray flux is small, as only cosmic rays skimming the solar surface and directed towards Earth (the solar limb) produce unabsorbed gamma rays Zhou:2016ljf (). However, the pioneering work of Seckel, Stanev, and Gaisser Seckel:1991ffa () (hereafter SSG1991) showed that solar atmospheric magnetic fields can redirect some cosmic rays from ingoing to outgoing before they interact. This greatly enhances the gamma-ray flux, as the entire Earth-facing solar disk can emit gamma rays. These gamma rays from the solar disk, which we denote as the solar atmospheric gamma rays (SA), are a new probe of cosmic rays in the solar system Potgieter:2013pdj (), solar atmospheric magnetic fields 2014A&ARv..22…78W (), and potentially new physics Leane:2017vag (); Arina:2017sng ().
In addition to SA, two other processes produce gamma rays from the solar position. Inverse-Compton (IC) scattering between cosmic-ray electrons and sunlight Orlando:2006zs (); Moskalenko:2006ta (); Orlando:2013pza () produces a diffuse gamma-ray halo that is most luminous just outside of the solar disk. Also, particles accelerated during energetic solar events, such as solar flares and coronal mass ejections Kafexhiu:2018wmh (), can produce GeV gamma rays. These, however, are transients events with much lower maximum energy (4 GeV (Fermi-LAT:2013cla, ; Pesce-Rollins:2015hpa, ; Ackermann:2017uer, )) compared to the SA flux.
The first evidence of SA were reported in Ref. Orlando:2008uk () using EGRET data, following limits from earlier studies ROG:ROG56 (); JGRA:JGRA13592 (). Shortly after the launch of the Fermi Gamma-Ray Space Telescope, the Sun was detected with the LAT (Large Area Telescope) (Abdo:2011xn (), hereafter Fermi2011) from 0.1 GeV to 10 GeV, with much better background separation. Importantly, this analysis found that the SA flux is significantly higher than even the magnetically enhanced flux predicted by SSG1991. Later, a new analysis (Ng:2015gya (), hereafter NBPR2016), using six years of Fermi data, confirmed this result, extended the SA detection up to 30 GeV, and found hints of a signal up to 100 GeV, with a spectrum that is significantly harder than predicted by SSG1991. In addition, NBPR2016 also found that the SA flux anticorrelates with the solar activity with a large amplitude, a result that was hinted at by Refs. Orlando:2008uk (); Abdo:2011xn (), but was not theoretically predicted.
In Refs. Ng:2015gya (); Zhou:2016ljf (), we showed that powerful new tests of SA will be possible with TeV gamma-ray experiments. The Sun is a unique target for large ground-based gamma-ray experiments, such as ARGO-YBJ Aielli:2006cj (), Tibet AS- Hibino:1988er (), HAWC Abeysekara:2013tza (); Abeysekara:2017mjj (), and LHAASO (Zhen:2014zpa (); He:2016del (), under construction). The TeV regime is important for SA, as it likely contains the transition between the low-energy regime, where solar magnetic fields enhance the SA production, and the high-energy regime, where such effects terminate. Even a limit on the solar flux at TeV energies will provide critical information on the magnetic environment that is responsible for SA production Ng:2015gya (); Zhou:2016ljf ().
The production of SA necessitates the production of solar atmospheric neutrinos Moskalenko:1991hm (); Seckel:1991ffa (); Moskalenko:1993ke (); Ingelman:1996mj (); Hettlage:1999zr (); Fogli:2006jk (); Arguelles:2017eao (); Ng:2017aur (); Edsjo:2017kjk (); Masip:2017gvw (). The Sun can be detected by large neutrino telescopes Ng:2017aur (), such as IceCube, where a search is ongoing icecubesolar (), and KM3NeT Adrian-Martinez:2016fdl (), in the future. A detection of these neutrinos will provide important and complementary information concerning the physics of solar cosmic-ray interactions. Additionally, these neutrinos are backgrounds for solar dark matter searches, and could constitute an impenetrable sensitivity floor Arguelles:2017eao (); Ng:2017aur (); Edsjo:2017kjk (). It is therefore important to have a robust prediction for the solar atmospheric neutrino flux, which can only be achieved with an understanding of the gamma-ray flux.
The continuous operation of Fermi provides a steady stream of data probing the evolution of solar gamma rays over the entire solar cycle. In this work, we revisit solar observations with an improved data set (Pass 7 to Pass 8) and an improved analysis technique, aiming to address several unanswered questions from NBPR2016. In particular, we aim to extend the detection of SA above 30 GeV, in search of a spectral softening that may reveal the termination of magnetic field enhancements from the SSG1991 mechanism. Second, we utilize the additional 3 years of data provided by Fermi to definitively test the anticorrelation between the SA flux and solar activity. Addressing both issues will lay the foundation for building a better model of cosmic-ray interactions in the Sun.
The paper is structured as follows. We describe the analysis method in Sec. II and the physics results in Sec. III. We investigate an apparent spectral dip in Sec. IV, and study the SA flux during solar minimum in Sec. V. Lastly, we discuss the implications of this work in Sec. VI, and then conclude in Sec. VII. Here we focus on the SA spectrum above 1 GeV and its time variation. In a companion work Linden:2018 (), we pay special attention to the highest part of the energy spectrum (above 100 GeV) and, for the first time, study resolved images of the solar disk in gamma rays.
Ii Fermi data Analysis
ii.1 Data selection
Our analysis procedure is similar to that of NBPR2016, except that we use an improved photon-position calculation. Specifically, we utilize 468 weeks of Pass 8 data collected between 2008-8-7 and 2017-7-27. The data are processed with the Fermi Science Tools version 10r0p5, and we utilize events belonging to the Source event class in our default analysis. In several later analyses and cross-checks we restrict our analysis to UltraCleanVeto class events, to further reduce cosmic-ray contamination.
To track the solar position and perform the event selection, we produce 4.2-hr temporal bins (40/week), during which the Sun moves degrees. We utilize gtselect to extract photons within a 9-degree region of interest (ROI) around the Sun, and then calculate the exact instantaneous distance of each photon from the solar position using the sunpos.py program included in the Fermi tools. This final step is a significant improvement with respect to NBPR2016, which utilized average positions in each temporal bin. In particular, calculating the exact solar distance of each photon is necessary to take advantage of the 0.1 angular resolution provided by Fermi at high energies.
We process the data using standard cuts in gtmktime, setting keyword DATAQUAL==1 and LATCONFIG==1 to ensure that the data quality are suitable for scientific analysis. The maximum zenith angle is set to be 90 to avoid the bright Earth limb. To reduce the background due to bright point sources and the diffuse emission from the Galactic plane, we select data when the Sun is at latitude . Given the size of the ROI, this effectively removes photons coming from . We conservatively removed data from any week (selected using the weekly photon data) that contains a significant solar flare (detected by Fermi at ) flare_list (). Compared to NBPR2016, only one more flare, in week labeled 326, is removed.
We calculate the Fermi exposure with gtltcube and gtexpcube2 using the P8R2SOURCEV6 instrumental response function with pixel size 0.10.1. The exposure is calculated in identical 4.2-hr temporal bins, which are summed to determine the total solar exposure. The fine temporal binning ensures that the Sun moves only 0.2 within each exposure bin. This is negligible compared to the angular dependence of the Fermi-LAT effective area. We additionally apply the same solar flare and Galactic plane cuts described above to the exposure calculation.
ii.2 Likelihood analysis
To determine the SA signal from the data extracted above, we follow a similar procedure as NBPR2016. In our default analysis, we focus on gamma rays above GeV, dividing the data into 4 equal logarithmic bins per decade. The 1 GeV lower bound is motivated by the Fermi-LAT point spread function; at lower energies, it is difficult to separate the SA from background emission (see Fermi2011 for an analysis down to 0.1 GeV).
There are two background components in the ROI that need to be taken into account for the SA analysis. The first is the IC halo from cosmic-ray electrons scattering with solar photons Orlando:2006zs (); Moskalenko:2006ta (); Orlando:2013pza (). The second is the collective gamma-ray emission from point sources and diffuse emission near the ecliptic plane, which all merge into a single diffuse component due to the motion of the Sun. We thus denote this background simply as the “diffuse background.”
To separate the signal and the two background components, we exploit their different angular distributions in the ROI. The SA flux, by definition, comes from the atmosphere of the Sun, and has an apparent angular radius , where is the angular distance from the center of the Sun.
For the IC halo, we assume the intensity falls linearly with . This follows from the assumption that the cosmic-ray electron density in the solar system is homogeneous, which was found to be reasonable Abdo:2011xn () and is conservative for our purpose. The IC intensity thus peaks at small angles, except for , where we set the IC component to 0, due to the strong suppression of IC emission stemming from the anisotropy of the solar photons and the shadow of the Sun Orlando:2006zs (); Moskalenko:2006ta (); Orlando:2013pza ().
|Energies||Total cts.||SA cts.||SA Flux|
Due to the motion of the Sun, the diffusion emission is strongly smeared out, and assumed to be constant over the ROI. This is supported by estimates of the background using the “fake Sun” method (see also Appendix A). We perform the same analyses and cuts for 3 ROIs that are displaced from the true time by , , and days, which is roughly , , and away from the true Sun position. As a result, we can accurately estimate the diffuse background intensity and avoid contamination from the IC halo. The diffuse background flux obtained from these three fakes Suns is used directly in our SA analysis.
We utilize a binned profile likelihood analysis to extract the SA signal from both background contributions (see NBPR2016 for a technical description). We divide the ROI into annuli of equal angular width. The widths of the annuli are chosen to be and for photons between – GeV and GeV, respectively. These widths are chosen such that the first angular bin is larger than the Fermi angular resolution in the energy range, and thus the SA signal is fully contained in the first angular bin to better than 95%. In each energy bin, we fit the SA component and the background components using their distinct angular distributions. Because we are only interested in the SA signal, we treat the IC halo and diffuse background normalizations as nuisance parameters. For the diffuse background, we add a Gaussian term with a 20% systematic error to incorporate the information we obtained from the fake-Sun studies. The normalization of the IC halo is completely free, thus our SA result also includes the uncertainty in the IC halo flux. The resultant SA flux, its uncertainties and significance are then obtained through the standard profile likelihood procedure Rolke:2004mj (); Cowan:2010js (). We note that the best fitting diffuse flux is consistent with the expectation from studies of the diffuse Fermi-background (Ackermann:2014usa, ).
Iii Sa Analysis Results
iii.1 Time-Averaged Spectrum
Table 1 summarizes the results of our SA analysis. The 9-year averaged SA flux is detected above 5 significance in all bins between 1–100 GeV, except for the 32–56 GeV bin, which is 4.1. This substantially improves on NBPR2016, where the SA flux was only detected at in the 32–56 GeV and 56–100 GeV bins.
Above 100 GeV, only a handful of events are observed in our data, and the SA detection is only marginally significant. However, in a companion paper (Linden:2018, ), we have shown that the emission above 100 GeV is highly significant during solar minimum, while no emission is detected over the remaining solar cycle. Thus, averaging over these two different time-periods artificially produces a softer flux in the full dataset. We compare these analyses in detail in Section V. Above 316 GeV, we observe zero events in our default data set, and thus derive the 90% upper limit using the total exposure within from the Sun, conservatively assuming no background events.
In Figure 1, we compare the SA flux with previous theoretical and observational results. We confirm the key findings of Fermi2011 and NPBR2016, which indicated that the SA flux significantly exceeds the theoretical predictions of SSG1991. Given that their prediction is already a significant enhancement over the limb-only flux Zhou:2016ljf (), it is surprising that the observed data are even brighter, and have a harder spectrum than that of the parent cosmic rays. We extend these results to higher energies, finding that the SA spectrum has a significantly harder flux than SSG1991 models up to an energy of at least 100 GeV. We additionally compare the SA flux with the estimated contributions from IC and diffuse backgrounds. The theoretically predicted IC flux is calculated using a modified version of the StellarICs code (Orlando:2013pza, ; Orlando:2013nga, , see Ref. Zhou:2016ljf () for a detailed description of inputs and modifications), while the predicted diffuse background is based on an identical analysis centered on the positions of three “fake Suns” that trail the position of the real sun by 90, 180 and 270 days. Our analysis clearly indicates that the SA flux is significantly brighter and has a harder spectrum than either the IC or diffuse components. Finally, we find evidence for a spectral dip in the 32–56 GeV bin, which we discuss in detail in Section IV.
iii.2 Time Variability
In Figure 1, we note that the extracted SA flux is significantly dimmer than the analysis of Fermi2011 and slightly dimmer than NPBR2016. This is not due to any tension between the analyses, but rather due to the intrinsic time-variability of the SA flux. This time variation was first robustly detected in NBPR2016, which found that the SA flux anti-correlates with solar activity. With 9 years of data, covering most of the solar cycle, we can now study the time variation in more detail.
We repeat the likelihood analysis described in Sec. II, but divide the gamma-ray emission into individual years, which we define to be 52 weeks, starting from the beginning of Fermi-LAT operations. Because the Sun traverses the ecliptic once in each period, every analysis window includes the same average diffuse background. To focus on the time-variation in each analysis period, we combine the SA flux into two logarithmic bins per decade spanning from 1–100 GeV (see Appendix G for an analysis with finer energy bins).
In Figure 2, we show the results of this analysis. The flux in each time period is shown as a fraction of the time-averaged flux at each energy. We overlay the averaged sunspot number 111http://www.sws.bom.gov.au/Solar/1/6 as a tracer for the solar activity. (See Appendix G for time variations of the total flux, defined in the next section.)
Between 1–10 GeV, we find that the SA flux clearly anti-correlates with the sunspot number, peaking during solar minimum and reaching a minimum flux near solar maximum. To formally test the statistical significance of this time-variation in each energy bin, we calculate the fit obtained under the assumption that the flux is time-independent, finding values of 91, 35, 8 and 6 with 9 degrees of freedom in each respective energy bin. These results thus provide clear evidence of time-variability below 10 GeV. While there is no statistically significant detection of time-variability at higher energies, we note that this is primarily due to insufficient statistics, rather than any significant shift in the degree of time-variability. At low energies, we find that the amplitude of this time variability is approximately 50% of the time-averaged flux, and find no evidence that the time variability decreases as a function of the gamma-ray energy. This final observation is surprising, given that higher-energy cosmic rays are less likely to be affected by solar magnetic fields. We investigate this further in the discussion.
Last but not least, as we have recently shown in Ref. Linden:2018 (), the degree of time-dependence is very different for the highest energy gamma rays (above 100 GeV), where 100% of the gamma-ray flux is observed during the first 1.5 yr of solar minimum activity. This degree of temporal variation is significant at nearly the 5 level.
Iv Investigating the 30–50 GeV Dip
The most striking feature of Figure 1 is the peculiar dip in the SA flux in the 32–56 GeV energy bin, which falls more than a factor of two below adjacent energy bins (see Table 1). The feature appears to be statistically significant (4 compared to adjacent bins), and may correspond to a dip in the underlying gamma-ray spectrum between 30–50 GeV. In Section VI, we discuss possible interpretations of this feature. Here, we attempt to examine its robustness and quantify its statistical significance.
In order to more sensitively and straightforwardly ascertain the statistical significance of the spectral dip, we refine our standard analysis method in three ways. First, we more finely bin the solar gamma-ray flux into eight logarithmic energy bins per decade, in order to more accurately determine the spectrum of the dip. Second, we simplify our analysis routine by examining the total gamma-ray flux within 0.5 from the Sun, rather than utilizing the SA analysis techniques that fit the gamma-ray components based on their angular profiles. This eliminates potential systematics in the likelihood analysis stemming from the assumed background emission morphologies. At high energies the effect of this choice is marginal, because the IC and diffuse backgrounds are highly subdominant to the true SA flux above 10 GeV (see, e.g., Figure 1).
Third, we enhance our analysis by including a new independent dataset from photons recorded when the Sun lies in the latitude range 5 (applying equivalent solar flare cuts). This latitude cut avoids only the brightest regions of the Galactic plane. While regions near the Galactic plane would suffer from significant diffuse contaminations at low energies, the contamination at energies above 10 GeV is minimal due to the soft spectrum of diffuse emission and the smaller angular acceptance window above 10 GeV. We note that the 5 sample has roughly 60% of the exposure of the sample, and thus provides a significant boost in statistics to the combined sample.
In Figure 3, we show the gamma-ray spectrum contained within 0.5 of the solar position in eight logarithmic energy bins per decade. Spectra are shown for photons recorded for both our default latitude cut, and the low-latitude analysis, along with the combined spectrum from both event selections. We find evidence for the spectral dip in both datasets and note that the gamma-ray dip now spans two distinct energy bins, spanning the same energy range of 32–56 GeV as above. This energy range is not altered by utilizing more finely binned energy spectra. We note that the energy range of the spectral feature is significantly larger than the 7% energy resolution of the Fermi-LAT at 40 GeV, and is thus consistent with a real spectral feature. We quantitatively test the effect of the Fermi-LAT energy resolution in Appendix D.
We apply several different techniques to calculate the statistical significance of this dip, and describe each in Appendix C. Here, we note that by fitting a broken power law to the remaining dataset (i.e., removing the two “dip” datapoints), we can calculate the fit of the dip to this model. We find that the dip is statistically significant at 4.6 in the dataset, 2.1 in the data, and 5.7 in the combined dataset.
In Figure 4, we examine the time evolution of the dip feature by dividing the combined () data into three periods spanning three years (156 weeks) of data. This corresponds to the periods 8/08-7/11, 8/11-7/14, and 8/14-7/17. We note that the dip is observed in each temporal bin, and the location of the dip does not appear to vary in energy. We find that the dip is most significant during years 1–3 (), but is weak in years 4–6 () and years 7–9 (). This may provide an indication that the dip is stronger during periods of reduced solar activity, a phenomenon that can be tested by upcoming solar-minimum data.
Finally, we note that there are several systematic issues intrinsic to solar analyses. In particular, the position of the Sun in Fermi-LAT instrumental coordinates is unique, as the x-axis of the satellite is oriented towards the Sun in order to maximize the solar panel efficiency 2012ApJS..203….4A (). It is possible that this could imbue a systematic feature in the response of the instrument, and cause such a dip feature. In the Appendix B, we further check the robustness this feature against several potential instrumental effects, including the azimuthal dependence of the LAT effective area and energy reconstruction, additional checks regarding the fake-Sun background estimation, and finally the event selection. We have not found any systematic and instrumental effects that can produce the spectral dip. (We do find a spectral feature, previously identified by the Fermi-LAT collaboration fermifl8y (); fermicaveat () that affects solar events 10 GeV, but does not influence the observed spectral dip feature.) We thus move on and discuss possible interpretations in Sec. VI.
V The solar minimum SA flux
In this section, we focus on the SA spectrum during solar minimum. As first pointed out in Ref. Linden:2018 (), the solar minimum SA spectrum is much harder and extends to higher energies than the 9-year averaged spectrum. This has important implications for the results of this study. For example, the soft gamma-ray spectrum GeV in our default 9-year analysis is actually due to the combined contributions of hard-spectrum gamma-ray emission during solar minimum, and very soft-spectrum gamma-ray emission over the remaining solar cycle.
Motivated by this result, we apply our SA analysis to the solar-minimum period, defined here as the 76 weeks of data between 2008-8-7 and 2010-1-21. In this section, we again make several modifications to our default event selection in order to better represent the solar-minimum data. First, we again choose to utilize a more-lenient solar latitude cut of to maximize the statistics during this smaller time window. Second, we allow the diffuse background to have a free normalization, rather than utilizing fake-Sun constraints that were optimized for the analysis. Lastly, we apply the UltraCleanVeto event selection cuts to minimize cosmic-ray contamination, a cut which reduces the exposure by 20–30%. This cut is employed to remove a systematic spectral feature that we (and the Fermi-LAT collaboration fermifl8y (); fermicaveat ()) identify in finely binned solar data at an energy of 10 GeV. Unlike the 30–50 GeV feature discussed here, the 10 GeV feature appears to be due to cosmic-ray contamination, and disappears entirely in the UltraCleanVeto class events. We discuss this feature in detail in Appendix B.2.
In Figure 5, we show the resulting solar minimum SA spectrum, noting three important results. First, we find close agreement with the Fermi2011 result between 1–10 GeV, verifying that differences between the SA flux in these analyses is due to the intrinsic time variability of the solar flux. Second, we find close agreement with the results of Ref. Linden:2018 () above 10 GeV, which is independent of this analysis and utilizes a slightly different data selection. Third, we find that the solar minimum analysis provides a more significant detection of SA emission above 100 GeV (3.7) than is found in our default analysis. The higher statistical significance comes from the utilization of gamma-rays during periods when the Sun resided at 530, and not including the non-Solar minimum period in the analysis, where no photons are observed while the expected background increases with time. If we also narrow the angular bin size to be , the detection significance reaches 5.2, in agreement with the findings in Ref. Linden:2018 ().
Finally, we note that the dip between 30–50 GeV is extremely significant (8.4) in this analysis. At first, this result may be confusing, as it exceeds the 5.7 significance of the dip over the full 9-year dataset. This disparity is partly due to the harder gamma-ray flux above the spectral dip during the solar minimum period, which provides support for a harder gamma-ray spectrum (without a significant cutoff) on both sides of the spectral dip, and prevents the broken power law from absorbing a portion of the spectral dip. The high significance is also partly due to the dip being less significant in periods outside the solar minimum (see Figure 4), which lowers the dip amplitude in the full 9-year dataset.
v.1 Implications for the SA production efficiency
The hard SA spectrum observed at solar minimum has important implications for the underlying SA production mechanism. Even without a detailed model, it is informative to compare the observed SA flux with the maximum flux allowed from cosmic-ray interactions in the solar photosphere. This extremal prediction can be robustly calculated by maximizing several relevant parameters in the gamma-ray production mechanism. These include the fraction of the solar surface that includes magnetic fields capable of reversing the direction of incoming cosmic-rays, the fraction of cosmic-rays that effectively have their directions reversed before undergoing a hadronic interaction, and the fraction of cosmic-rays which subsequently undergo a hadronic interaction before escaping the solar surface. We describe this calculation in detail in Appendix E of Ref. Linden:2018 (). By setting these three terms to unity, we can calculate an upper bound on the gamma-ray flux, which can be approximately parameterized as above 1 GeV.
In Figure 5, we compare the solar minimum SA spectrum to this maximal gamma-ray flux, finding that the observed SA flux has a significantly harder spectrum () relative to the cosmic-ray spectrum (), allowing it to approach the upper bound at high energies. This suggests two seemingly contradictory results. On one hand, the SA production efficiency must be near unity at high energies to produce the observed flux. On the other hand the SA production efficiency must be increasing rapidly at high energies in order to maintain the hard SA spectrum. Moreover, while it is not surprising that the SA production efficiency is energy dependent, it is somewhat surprising that the efficiency increases with energy, given that it is more difficult to confine and reflect higher energy particles in the solar magnetic field. It is also important to note that the SSG1991 model, if continued to high energy, follows the cosmic-ray spectrum. Thus the gamma-ray data implies qualitatively different physical effects compared to SSG1991.
The proximity of the SA flux at solar minimum to the maximal cosmic-ray flux also presents significant issues. Above 100 GeV, the SA flux is close to 30% of the cosmic-ray upper bound. It is not clear how to achieve such an efficiency in a realistic setup. In Ref. Linden:2018 (), we discuss several possible methods of further enhancing the SA flux, including magnetic focusing from the large scale magnetic field, long-term cosmic-ray confinement in the coronal magnetic field, anisotropies in the SA emission mechanism, or some unknown intrinsic solar gamma-ray production mechanism.
Finally, because the highest-energy gamma-ray photon from a typical hadronic interaction has an energy of 0.1, the observation of a 100 GeV SA flux suggests that the enhanced gamma-ray production mechanism applies to even TeV cosmic rays, at least during solar minimum. This observation is partially supported by the cosmic-ray shadow measurements by Tibet AS measurement Amenomori:2013own (), which showed that solar coronal magnetic fields can affect cosmic rays 10 TeV. These will have important implications for searches of solar atmospheric neutrinos and dark matter signals from the Sun; we discuss further below.
v.2 Detection Prospects at TeV Energies
The high-energy SA flux during solar minimum has immediate implications for TeV solar gamma-ray searches. In this regime, only ground-based air shower experiments and water Cherenkov telescopes, such as ARGO-YBJ Aielli:2006cj (), Tibet AS-gamma Hibino:1988er (), HAWC Abeysekara:2013tza (); Abeysekara:2017mjj (), and LHAASO Zhen:2014zpa (); He:2016del () are capable of solar observations. Figure 5 also compares the SA flux with the 1-year differential sensitivity of HAWC and LHAASO. If the hard spectrum observed during solar minimum continues to the TeV range, then it will fall within the reach of HAWC’s 1-year sensitivity. Given that HAWC is already operating, the current solar minimum (Cycle 25, beginning presently) provides the perfect opportunity for HAWC to test the TeV SA flux. We also show the preliminary upper limits of SA from ARGO-YBJ argo:2016 () and HAWC Nisa:2017xpf (). The latter is obtained with only 1 month of data, and the sensitivity is expected to improve significantly.
In the near future, LHAASO is expected to provide even better sensitivity and a lower energy threshold. Excitingly, the projected sensitivity is already reaching the solar minimum spectrum at around 200 GeV. LHAASO may even be able to probe the considerably softer SA flux outside of the solar minimum.
Vi Discussion and Implications
vi.1 Interpretations of the 30–50 GeV dip
In this work we discover a spectral dip between approximately 30–50 GeV, with a local significance exceeding . This is surprising given that the cosmic-ray energy spectrum is incredibly smooth. Moreover, though SSG1991 predicts a spectral cutoff between 30–2000 GeV, a hard spectral component above the cutoff is not expected. While we cannot rule out any unknown systematic effects that may produce this feature, it has survived several sensitive tests detailed in the Appendix. In Appendix D, we find that the spectral dip is consistent with approximately a 50% reduction in the gamma-ray flux over an energy range of approximately 20 GeV, which is then smeared by the Fermi-LAT energy resolution. These numbers guide our search for a potential source of the spectral dip feature.
In this subsection, we consider and quickly investigate several potential interpretations of the dip feature, with additional details presented in the appendix.
vi.1.1 Hypothesis I: Two Cosmic-Ray Components
One model capable of producing a spectral dip includes two separate cosmic-ray proton components, which dominate SA emission below and above the spectral dip. Given that the interstellar cosmic-ray spectrum is smooth, and the efficiency of SA production is nearly maximized at high energies, this split would be most reasonably produced by magnetic-field structures that inhibit cosmic-ray interactions within a specific, and narrow, energy range. For example, it is possible that one class of magnetic field structures are responsible for gamma-ray production below 30 GeV (e.g. flux tubes, as suggested by SSG1991), which become ineffective at high energies and cause a cutoff. The rise in the gamma-ray flux above 50 GeV would then be caused by a separate class of magnetic structures that are only accessible to high-energy cosmic-rays.
However, we find that it is difficult to produce both the depth and the width of the observed spectral dip with a two-component proton spectrum. This is due to the fact that the gamma-ray kernel for a given proton interaction is broad, which we detail in Appendix. F. The kernel may be sharpened, if only a narrow kinematic region is considered (e.g., a tiny angular cone from the parent proton direction). However, it is not clear how such a scenario can be realized in the Sun, where the cosmic-ray direction is expected to be roughly isotropic, and the gamma-ray efficiency must be nearly maximal, in order to explain the high SA flux.
vi.1.2 Hypothesis II: Gamma-Ray Absorption
Another possible mechanism capable of producing the spectral dip is the preferential absorption of gamma-rays in the 30–50 GeV energy range. Gamma rays can be absorbed by gas through the interaction , with an absorption length of , which is comparable to that of p-p interactions Patrignani:2016xqp (). However, above the threshold ( 1 MeV), the cross section is smooth, and thus unlikely to cause a dip. While (e.g., ) could produce a dip feature, the cross sections are too small Patrignani:2016xqp (). Similar processes occur for electron targets with higher energy thresholds, but the cross-sections are even smaller. Absorption on photon targets, i.e. interaction, may produce a dip feature. The reaction threshold is ) GeV, taking the typical target photon energy at the photosphere ( 1 eV). The threshold can be lower and thus match the energy of the dip if the interaction occur below the photosphere, where the Sun is hotter. However, in this region the matter density would be orders of magnitude higher than the photon density, and thus matter effects should provide the dominant source of gamma-ray absorption. Therefore, we find that it is unlikely for gamma-ray absorption to produce the observed dip feature.
vi.1.3 Hypothesis III: Solar Neutrons
One intriguing possibility is that the dip is produced by solar neutrons, rather than solar gamma-rays. Unlike the case of more distant sources, solar neutrons can travel unimpeded to Earth before decaying. Assuming that the SA flux is produced by hadronic processes, the solar neutron flux recorded by the Fermi-LAT should be similar to the gamma-ray flux (see e.g., SSG1991). Solar neutrons have, in fact, been observed by Earth-bound observatories and correlated with gamma-ray emission during solar flares 2016SoPh..291.1241M (); Muraki:2017mhm ().
While the Fermi-LAT event reconstruction algorithms are finely tuned to eliminate the large contaminating flux from charged cosmic-rays, the efficiency in rejecting high-energy neutrons is not publicly available. Additionally, it may be difficult to definitively ascertain the efficiency of this cut, due to the very low flux of solar neutrons from non-solar or terrestrial sources. Furthermore, while analysis cuts of gamma-rays and hadrons have been carefully tuned to avoid inducing spectral features in the resulting gamma-ray spectrum, sharp features may potentially be induced from neutron sources. In Appendix E, we test several potential event-level quantities that may potentially distinguish neutron and -ray events, finding no obvious evidence for a neutron component. However, without utilizing non-public instrument-level data, we do not have the ability to rigorously test this hypothesis.
Assuming for a moment that solar neutrons are present in the Fermi-LAT data, we note that these events would be reconstructed utilizing analysis algorithms tuned to accurately determine the energy of gamma-ray showers. Thus, the energy reconstruction of neutron events may be highly biased. In this scenario, the observed “dip” in the solar gamma-ray spectrum may result from a population of 30–50 GeV neutrons, which are mis-identified by the Fermi-LAT analysis algorithm as 50–80 GeV gamma-rays. This would explain several features of the gamma-ray dip, including (1) its sharp spectral fall-off and rise, and (2) the large gamma-ray flux at energies immediately above the spectral dip. Additionally, if the neutron rejection efficiency were found to fall smoothly as a function of energy, it may explain the hard-spectrum of observed SA emission, although the neutron rejection efficiency would need to be fine-tuned in order to replicate power-law observations over more than two decades in energy.
We note, however, that this explanation does not explain the time-variation of SA emission, or the significant spectral differences observed from SA emission above 100 GeV during solar minimum. As pointed out by Ref. (Linden:2018, ), photons observed during solar minimum have unique spectral and morphological characteristics compared to photons observed during the remaining solar cycle. These morphological shifts, in particular, cannot be modeled by misidentified neutrons, and instead indicate a significant shift in the solar gamma-ray generation mechanism.
vi.1.4 Hypothesis IV:
Statistical fluke or Time-dependent instrumental effects
We note that the dip is most significant during the solar minimum and becomes much weaker in later time periods. While this may be a hint that the dip is related to solar minimum activity, we cannot fully rule out the possibly of an unlikely statistical fluke, given that only one solar minimum was observed by Fermi.
More importantly, we note that the first solar minimum corresponds to data taken just after the launch of the Fermi-LAT. Though contrived, we cannot rule out the possibility of a time-dependent instrumental effect that causes a dip only in the instrumental phase space of the Sun. The time-dependent effect could be caused by changes of Fermi operation or event reconstruction during its lifetime. For example, in September 2009, the rocking angle of Fermi was increased from to 2012ApJS..203….4A (). However, it is not clear show such changes would produce a spectral dip feature. Fortunately, these two possibilities can be tested by Fermi observations of the Sun during the on-going solar minimum.
vi.1.5 Hypothesis V: Intrinsic Solar Emission
Lastly, we must also note that we are unaware of any intrinsic solar mechanism that can produce such energetic gamma-ray emission. But in light of a two-component SA scenario, we cannot rule out any intrinsic solar component interpretation.
vi.2 Time variations of SA flux
In this work, we have observed moderate time variations () in the SA flux between 1–10 GeV. Interestingly, the amplitude of this variation does not seem to vary significantly with energy. Naively, one would expect that the amplitude of the time variation decreases with energy, as seen in the solar modulation of cosmic rays observed at Earth Potgieter:2013pdj (). However, simple models, where the amplitude of the time variation is proportional to either or , are ruled out by our SA analysis. The relatively energy-independent time variation amplitude could be a hint that additional solar modulation from within the Earth’s orbit plays a relatively minor role in affecting the cosmic-ray spectrum at the solar surface. Alternatively, the observed relationship could be produced via two competing effects. For example, if the efficiency of solar atmospheric magnetic fields in producing outgoing gamma rays from incoming cosmic rays increases with energy, this could offset the effect of solar modulation in decreasing the cosmic ray flux. Finally, the observed time variation could be dominated by other time-varying solar properties, such as a time-dependent scale height for cosmic-ray reversal that affects cosmic rays equally at all energies.
We note two potential observables that could illuminate these possibilities. First, a detailed analysis of the IC halo may be able to probe the effect of cosmic-ray electron propagation inside the Earth’s orbit, isolating any effects due to the solar atmosphere. Second, while the solar minimum SA flux is bright up to 200 GeV, no GeV events were observed outside the solar minimum. This extreme time variation may provide important clues for the underlying gamma-ray production mechanism.
vi.3 Implication for solar atmospheric neutrinos
Understanding the production mechanism of the SA flux is extremely important for solar atmospheric neutrino searches by IceCube and KM3NeT. Normally, the calculations of solar atmospheric neutrinos assume that the Sun has zero magnetic fields Ingelman:1996mj (); Arguelles:2017eao (); Edsjo:2017kjk (), implying that Earth-bound neutrinos are produced from the far side of the Sun. Importantly, neutrino absorption in the Sun starts to be important above 100 GeV. Even in this case, the harder solar atmospheric neutrino spectrum (compared to the Earth’s atmospheric neutrino background) makes the Sun a potentially detectable source for large neutrino telescopes above 1 TeV Ng:2017aur (). It was, however, not clear how magnetic fields could affect the solar atmospheric neutrino flux.
The observed GeV SA flux during the solar minimum implies that solar magnetic fields must be able to significantly affect cosmic rays TeV to boost the gamma-ray production. If there is a large component of reflected TeV cosmic rays that interacts in the atmosphere, this may enhance the solar atmospheric neutrino flux by producing neutrinos in lower density regions that have a longer charged pion interaction length, and where the neutrinos suffer negligible absorption.
Importantly, the observed solar minimum SA flux provides a model-independent estimate of the minimal solar atmospheric neutrino flux at this energy range. This estimate assumes no gamma-ray absorption. If gamma rays are efficiently absorbed by the solar atmosphere, it would increase the predicted neutrino flux.
We take the solar minimum SA flux above GeV to be , with . The corresponding muon neutrino flux is about a factor of smaller, thus . This takes into account that the number of photons is roughly the same as each flavor of neutrino (including antineutrinos) from hadronic interaction, but photons have twice of the averaged energy Kistler:2006hp (). Approximating the neutrino effective area of IceCube as Aartsen:2016zhm (), and considering from 100 GeV to 1 TeV, the number of events expected by IceCube is about 1 event per 1.5 year. It is thus extremely interesting to see if IceCube can detect these neutrino in the upcoming solar minimum.
If the high-energy SA flux is transient in nature (e.g., correlating with the formation and destruction of specific magnetic field structures), then a coordinated search using IceCube, Fermi and HAWC will be extremely beneficial. In particular, employing variability information can improve the sensitivity of IceCube by allowing it to relax some data selection criteria to potentially lower the energy threshold, while maintaining a high signal-to-background ratio within a smaller time window. It is fortunate that IceCube and HAWC are located in different hemispheres, as IceCube is significantly more sensitive to solar neutrinos when the Sun is below horizon.
vi.4 Implication for solar dark matter searches
The detection of GeV SA flux could significantly affect solar dark matter searches with neutrino telescopes. The fact that magnetic fields strongly affect cosmic rays above 1 TeV adds considerable uncertainty to current solar atmospheric neutrino flux calculations. If the neutrino flux were similarly enhanced, compared to the calculations without magnetic fields, this could raise the dark matter sensitivity floor Arguelles:2017eao (); Ng:2017aur (); Edsjo:2017kjk (), and reduce the dark matter parameter space that can be probed by neutrino telescopes Bagnaschi:2017tru (); Costa:2017gup ().
Most dark matter models do not produce a significant neutrino flux above 1 TeV, due to neutrino absorption within the solar core. However, if dark matter accumulated in the solar core annihilates first into some metastable mediator, that mediator could subsequently decay outside the solar core and produce observable neutrinos, gamma rays, and other messengers at TeV energies (for some examples, see Refs. Bell:2011sn (); Feng:2016ijc (); Leane:2017vag (); Arina:2017sng ()). Correlated high-energy neutrino and gamma ray observations were proposed to search for these metastable mediator dark matter model, as the cosmic-ray induced gamma-ray flux was expected to be negligible at TeV energies Leane:2017vag (). If solar magnetic fields can enhance both gamma-ray and neutrino production above TeV, it would be difficult to differentiate this emission from metastable mediator dark matter signals, at least in the TeV range.
Lastly, while it is tempting to associate the SA component above the 30–50 GeV dip with dark matter or another new physics origin, it is not clear how exotic physics can accommodate the extreme time variation and morphological shift Linden:2018 () observed in SAevents GeV. On the other hand, the extreme time variation in the high-energy SA flux potentially provides a method for removing the astrophysical contribution to reveal any underlying dark matter component.
Vii Conclusions and Outlook
In this work, we analyzed 9 years of Fermi data to detect the SA flux and study its spectrum and time variation. We detect significant SA emission from 1 GeV up to 200 GeV, extending our previous results in NPBR2016 that detected SA emission only up to 30 GeV.
Using data covering most of the solar cycle, we find that the SA flux anticorrelates with solar activity between 1–10 GeV. The SA flux peaks during the 2008 solar minimum, and decreases as the Sun becomes more active. Additionally, we observe, for the first time, the increase of SA flux corresponding to the decreasing solar activity since 2013. We find that the amplitude of this time variation does not change significantly with energy between 1–10 GeV, which may shed light on the underlying gamma-ray production mechanism.
Interestingly, we discover a significant spectral dip between 30–50 GeV in multiple datasets. We find no evidence that this dip is a systematic or instrumental artifact, and find no obvious theoretical interpretation. The dip appears to be more significant when the Sun is less active, thus Fermi observations during the upcoming solar minimum may shed light on the nature of this feature.
Motivated by our companion paper Linden:2018 (), we compute the SA flux during the Cycle 24 solar minimum. The solar minimum SA flux exhibits a hard and bright spectrum up to at least 200 GeV. This hard SA component could be tested by Fermi and HAWC during the upcoming Cycle 25 solar minimum, which is beginning presently. In the future, LHAASO can probe the SA flux at TeV energies with even better sensitivity.
The hard spectrum of the solar minimum SA flux suggests that solar magnetic fields affect cosmic rays above 1 TeV in the solar atmosphere, which suggests that magnetic fields may significantly affect the expected solar atmospheric neutrino flux. In addition, the observed SA flux during solar minimum implies a minimal detectable neutrino flux, which may be constrained by IceCube. A coordinated search of IceCube, HAWC, and Fermi data may improve the sensitivity of these observations and test whether this component is transient in nature.
The results of this work lay the foundation for building a more accurate model for cosmic-ray induced high-energy emission from the Sun. This is crucial for a robust prediction of high-energy solar neutrinos, which is, in turn, important for solar dark matter searches. Ultimately, high-energy emission from the Sun may provide a novel probe of the solar atmosphere, cosmic rays in the solar system, and new physics.
We thank members of the Fermi-LAT collaboration, and in particular Regina Caputo, for discussions concerning potential instrumental systematics near the solar position. We also thank Markus Ackermann, Keith Bechtol, Ofer Cohen, and Carsten Rott for helpful discussions. TQW is supported by the National Natural Science Foundation of China under grants No. 11547029. KCYN is supported by Croucher Fellowship and Benoziyo Fellowship. TL, BZ, and AHGP are supported in part by NASA Grant No. 80NSSC17K0754. BZ is also supported by a University Fellowship from The Ohio State University. JFB is supported by (and BZ is partially supported by) NSF grant PHY-1714479.
Here we provide more details concerning the analysis techniques utilized throughout the main text, and perform systematic checks to evaluate the possibility that instrumental systematics affect our determination of the SA flux and spectrum. In Appendix A, we investigate potential systematics stemming from the non-standard data analysis techniques required to evaluate the flux and spectrum of the Sun, which is rapidly moving compared to diffuse backgrounds. In Appendix B, we investigate potential systematic issues in the Fermi-LAT effective area or energy reconstruction stemming from the unique position of the Sun in instrumental phase space. In Appendix C, we describe the methodology used to quantify the significance of the 30–50 GeV dip in the SA spectrum. In Appendix D, we check the effect of instrumental energy dispersion regarding the gamma-ray absorption hypothesis of the dip. In Appendix E, we check the quality of the photons regarding the solar neutron hypothesis. And in Appendix F, we investigate the potential two-component cosmic-ray interpretation of the dip. Lastly, in Appendix G, we show the time variation of the SA flux in finer energy bins.
Appendix A Systematic Effects in the Data Analysis Methodology
a.1 Background Estimation Using “Fake Suns”
In this section, we examine the gamma-ray flux and spectrum from three “fake Sun” models, which trail the position of the Sun by 3 months, 6 months, and 9 months, respectively. This test serves two purposes. First, it provides an accurate estimate of the diffuse background at the solar position, because each fake Sun moves (over the course of a year) through the same array of coordinates as the real Sun. Because the vast majority of gamma-ray emission occurs in steady state, this provides an exposure-weighted calculation of the background flux. Second, this test allows us to evaluate our data extraction methodology, to ensure that our analysis of a quickly moving solar region does not produce any spectral features in the position of power-law backgrounds.
In Figure A 1, we compare the SA flux calculated in the main text with the calculated gamma-ray flux from each fake-Sun position. We immediately note two conclusions: (1) the diffuse gamma-ray emission accounts for approximately half of the observed gamma-ray flux near the solar position at low-energies, but contributes only 10% of the total emission at high energies. These results are consistent with our SA analysis (shown in Table 1) (2) there is no spectral feature in the diffuse gamma-ray flux of the fake-Sun positions in the 30–50 GeV energy range, which argues against any contribution from either diffuse emission or our analysis routines to the observed solar spectral dip.
a.2 Background Estimation Using Stacked “Fake Suns”
At high energies, the statistical information that can be obtained from the fake-Sun positions analyzed in Section A.1 is minimal, due to the very small diffuse gamma-ray flux at high energies. In particular, the small statistical sample limits our ability to conclusively determine whether or not a dip in the gamma-ray spectrum exists between 30–50 GeV in our fake-Sun sample.
In this section, we increase the statistical sample by evaluating the total emission from 177 fake sky positions, including “fake-Suns” that trail the position of the real sun by between 6–360 days in two-day increments. We choose the two degree increment because the Sun moves approximately 1 per day, and remove positions within 5 of the Sun due to the non-negligible solar ICS component close to the real Sun. Because this analysis is focused on the origin of the 30–50 GeV dip, we use a consistent radial cut of 1.0 at all energies in this analysis.
Unfortunately, due to the large number of solar positions we analyze, it is not computationally feasible to calculate the Fermi-LAT exposure at each fake-sun position. Thus, we plot the spectrum in units E dN/dE, utilizing the result from Section B.1 which indicates that their is no peculiar energy dependence in the Fermi-LAT effective area near the position of the gamma-ray dip.
In Figure A 2 we show the resulting spectrum from our analysis, finding no evidence for a unique spectral feature between 30-50 GeV. The combination of this test, alongside the fluxes of individual fake suns in Appendix A.1, indicate that our analysis technique does not introduce any spectral dips in the 30-50 GeV energy range.
Appendix B Systematic Effects in the Instrumental Response Near the Solar Position
b.1 The -dependence of the Effective Area
We utilize the Fermi-LAT tools gtltcube and gtexpcube2 to calculate the energy-dependent effective area and exposure of solar observations. In our default analysis, we follow the standard procedure, which finely bins the effective area in the angular coordinate denoting the distance of the observed photon from the Fermi-LAT boresight (), but does not finely bin the effective area in the second angular coordinate, which describes the recorded angle of photons with respect to the Fermi solar panels (). For standard target analyses, these choices are correct - because the effective area varies strongly with and only marginally with .
However, because the Fermi-LAT solar panels are typically aligned with the Sun, it occupies a unique position in -space. Thus, it is likely that any small -dependent systematic will appear first in solar observations. In Figure A 3, we plot the distribution of solar photons that pass our angular and energy reconstruction cuts in instrumental phase-space, compared to a randomly selected sample of photons located more than 5 from the solar position. We find that the Sun has a slightly-skewed distribution in -space, but note that we carefully take this into account by weighting the effective area in bins of cos() = 0.025. However, the -distribution of solar photons differs drastically from the average, and this is not taken account in our default analysis. In fact, the distribution is somewhat more peaked than portrayed in Figure A 3, as nearly all solar events are located within 3 of = 0 (which gives cos() 0.998).
While the Fermi-LAT collaboration notes that the overall change in the effective area at 10 GeV as a function of is only 1% fermi_performance (), they do not discuss the energy-dependence of this effect. Here, we utilize the hidden parameter phibins in gtltcube to directly calculate the -weighted energy-dependent exposure of solar photons. This accounts for any known systematic issue affecting the Fermi-LAT effective area or energy reconstruction near the solar position that was not taken into account in our standard analysis. In Figure A 4, we plot the Fermi-LAT effective area as a function of energy for models where we set phibins=1 (our default assumption) as well as phibins=10 and phibins=30. While even the choice of 30 -bins produces 12 bins that are somewhat larger than the distribution of solar photons in -space, we found that this was the largest computationally feasible measurement. Moreover, because events near = 0 pass directly through the flat side of the instrument, we expect the -dependence to be mild at values very close to = 0 itself.
We note two effects. First, the calculated effective area increases by about 4% when the analysis is increased to 30 bins. This result is similar to the Fermi-LAT collaboration analysis, and implies that the calculated solar gamma-ray flux is decreased by about 4% relative to the values calculated in the main portion of the text. Second, the modification to the Fermi-LAT effective area is relatively energy-independent, inducing spectral changes only at the level of 1%. The result of this analysis indicates that no known systematic issue in the Fermi-LAT event reconstruction at values near = 0 contributes to the dip in the solar gamma-ray spectrum.
b.2 Spectral Features in Events Near =0
In the previous section, we ruled out the possibility that “known” features in the -dependence of the Fermi-LAT effective area could account for the observed spectral-dip in the Solar gamma-ray flux. However, it is also possible that “unknown” features — those not represented in models of the Fermi-LAT effective area or energy dispersion — could account for the spectral dip observed in the solar gamma-ray flux. Because most sources have relatively uniform exposures in , it is possible that a significant effect has gone unnoticed in previous analyses of the Fermi-LAT data – only to appear first in high-energy analyses at the solar position.
We develop the following test of this hypothesis. We first extract every gamma-ray photon that passes our default event class and temporal cuts with the following characteristics: (1) a recorded energy exceeding 1 GeV, (2) a distance of at least 5 from the contemporaneous solar position, (3) an observed -angle within 3 of = 0. These cuts produce a dataset of 223,707 events with similar instrumental reconstructions as solar photons, but which do not have a solar origin. For each event, we select a single random “partner” event with the following characteristics: (1) an energy exceeding 1 GeV, (2) an observed -angle exceeding 3 from = 0.0, and (3) a small angular separation from the original event in r.a. and dec.
The final selection criteria ensures that the average spectrum of all source and partner events is similar, as they are produced in equivalent regions of the Fermi-LAT sky. We pick angular separations following a tiered approach, we first randomly select any partner event within 0.1 of a given 0.0 event, noting that this works for 90% of our population and the separation is much smaller than the Fermi-LAT angular resolution, indicating that the origin of these two events is nearly identical. If no such photon is found, we widen our angular acceptance to 0.2 and then 0.5 to ensure that every photon has a partner. In combination with our observation (in the previous section) that the -dependence of the instrumental effective area is energy-independent to within 1%, this implies that the counts spectrum of gamma rays from both our source events and our random “partner” events should be identical.
First, we verify that this hierarchical event-selection procedure does not induce any spectral features into the 0.0 and partner events. We test this by repeating the analysis routine 10 times, with randomly chosen seed events. We find no spectral differences between the parent and partner events at the level of 1%, indicating that our hierarchical event selection technique does not introduce any systematic biases into the energy-distribution of events.
In Figure A 5 (left) we show the results of this analysis. Overall, we find that the spectra of the 3 and 3 event selections are similar. However, there is a highly statistically significant feature extending from approximately 3 GeV to nearly 25 GeV in energy. In this range, we see an initially smaller population of 3 gamma rays shift to a larger population of 3 gamma rays at an energy of 8 GeV. We note that, because there are an equal number of events in each selection, an excess of photons at one energy must correspond to a deficit elsewhere. The magnitude of this effect is nearly 20% at 10 GeV. In the center panel, we show that this effect disappears when we select only the subset of gamma rays that pass the UltraCleanVeto event selection, while in the right panel, we show that this effect is more extreme for gamma rays that pass the Source, but not the UltraCleanVeto event classes. This indicates that the spectral feature is likely due to differences in the cosmic-ray rejection (or cosmic-ray energy reconstruction) between events observed near = 0 and events observed elsewhere in instrumental coordinates.
We note that this spectral feature could contribute to the observed solar gamma-ray spectrum at 5-15 GeV, shown in Figure 1. In particular, there are hints of a rise in the observed solar gamma-ray spectrum between 5-8 GeV. However, we do not see any evidence in Figure A 5 for a statistically significant spectral change between 30–50 GeV for photons observed near = 0. Because the photon count is equivalent in both datasets, we note that spectral changes should appear as changes in the gamma-ray ratio between selection cuts, rather than in the overall normalization of this selection ratio.
One potential explanation for the peculiar spectrum of gamma rays observed near = 0 is the influence of the Sun itself. This would be a particularly concerning possibility, as it would have an even more acute effect on our solar analysis compared to = 0 events, many of which are not located near the Sun. While the mechanism for such an effect is not known, the influence of solar radiation may potentially affect the performance of the Fermi-LAT. In Figure A 6 (left), we test this possibility by comparing the spectrum of photons observed between 5–15 from the contemporaneous solar position with photons observed more than 90 from the solar position, finding no statistically significant evidence for such an effect. We note that this observation tightly constraints the systematic shown in Figure A 5, because many photons located between 5–15 from the Sun have rather small values. Thus the systematic effect is strongly peaked only toward events with 3.
Additionally, in the center panel of Figure A 6, we examine the distribution of events recorded near = 180, finding a similar spectral feature at 10 GeV. Because = 180 events tend to point away from the contemporaneous solar position, but also are observed in an orientation normal to the solar panel plane, this implies that the systematic issue may appear in certain locations of instrumental -space. Finally (right), we test photons recorded near = 135, which occupies a corner of the Fermi-LAT. We do not find the same systematic issue, but do observe a small dip in the gamma-ray spectrum between 10-25 GeV with an amplitude of 5%.
The conclusion of this analysis is uncertain. There does appear to be a statistically significant effect, with an amplitude of nearly 20%, affecting the reconstruction of Source class (but not UltraCleanVeto class) events near = 0. The fact that this event disappears at angles between 5–15 from the Sun, suggests that this effect disappears very rapidly for larger values of . While this effect may influence our Solar spectrum near an energy of 10 GeV (where a similar feature is observed), there does not appear to be any lingering effect between 30–50 GeV. Moreover, the amplitude of this effect is insufficient to explain the much larger spectral dip observed in our data. While more analysis is necessary to understand this effect – at present we find no indication that it explains the observed data. In the next section, we verify this result by re-running our analysis technique on the UltraCleanVeto event class, which shows no evidence for any peculiar -dependent effects.
b.3 Analysis of UltraCleanVeto Events
In the last subsection, we showed that a systematic spectral feature appears to be induced for events near 10 GeV in the Source event class, but does not appear to be present in events which pass the more stringent UltraCleanVeto event reconstruction. Here, we re-run our analysis on events that pass the more stringent UltraCleanVeto cuts. This decreases the size of our statistical sample, but likely removes any remaining effects due to cosmic-ray misidentification.
In Figure A 7 we show the resulting gamma-ray spectrum, utilizing the full angular ROI (5) with solar flare cuts implemented. In order to minimize cosmic-ray background, we repeat the same analysis routine, restricting our analysis to photons which pass the UltraCleanVeto analysis cuts. This event class places a significantly more stringent rejection cut on possible charged cosmic-ray contamination, compared to the Source Event class utilized in the main text – but does so at the cost of a 20–30% rejection of gamma-ray events. This latter effect can be accounted for in the re-calculation of the instrumental effective area. Our analysis shows that the spectral dip remains a resilient feature in this event class, providing strong evidence that charged cosmic-rays events are unlikely to be responsible for the gamma-ray spectral dip.
b.4 Analysis of Events Immediately Outside of the Solar Disk
Thus far, we have found a significant systematic artifact in the reconstruction of Fermi-LAT events near = 0, which affects photons near 10 GeV, but does not affect events near the 30–50 GeV dip. Here, we produce another test, by calculating the spectrum of gamma-rays observed immediately outside of the solar disk. In particular, we analyze the gamma-ray spectrum in an annulus between 1–3 from the center of the Sun. This lies outside the Fermi-LAT PSF for all but the lowest energy gamma-rays, and is thus relatively free of contamination from the solar disk. Instead, these gamma-rays are expected to be produced by a combination of solar IC emission and diffuse background sources. Thus, this emission would not be expected to show any matching spectral features at either 10 GeV or between 30–50 GeV.
In Figure A 8 we show the resulting gamma-ray spectrum between 1–3 from the solar center, using Source class photons over the full 9 yr exposure. We find clear evidence for the 10 GeV spectral feature found in Figure A 5, providing additional evidence for its systematic origin. On the other hand, we find no evidence of a 30–50 GeV spectral dip in this spectrum. The statistical sample of these events is large enough that the observed emission is statistically inconsistent with the spectral dip observed across the solar disk. This implies that if the 30–50 GeV spectral dip is a systematic feature, the effect is confined to within 1 of the solar position, which we believe to be unlikely. Thus, this argues against a systematic origin for the 30–50 GeV feature.
b.5 Analysis of the Earth Limb
The brightest -ray source observable by the Fermi-LAT is the Earth limb, which produces gamma-ray emission through the interaction of cosmic-rays with the solar atmosphere. Because this emission source is known to be hadronic and feature-less, it serves as an optimal test location for instrumental artifacts. In this section, we compare the spectrum of Earth-limb events (those recorded with zenith angles exceeding 110) from with events at larger azimuthal angles .
In Figure A 9, we show the results of this analysis, which provides no evidence for a spectral feature between 30–50 GeV in either dataset. We do note a slight (but highly statistically significant) spectral softening in the data below 10 GeV, compared to background events. This spectral feature may result from the different exposure of = 0 events towards the Earths surface, stemming from the preferred orientation of the solar panels towards the Sun. While we are unable to conclusively test this hypothesis, we stress that this feature is extremely smooth, most prominent at relatively low energies, and does not have an amplitude capable of significantly affecting the measurements presented in the main text.
In this section, we have discussed both systematic issues relating to the instrumental effective area of Fermi-LAT events recorded near the solar position, as well as systematic issues relating to the data-analysis procedure utilized in the text. While we have found a peculiar systematic issue involving events recorded near = 0 and near an energy of 10 GeV, we have found no plausible systematic issue which might produce a dip in the solar gamma-ray spectrum in the 30-50 GeV energy range under investigation. In addition, we do not detect any 30–50 GeV dip in photons recorded very close to the contemporaneous solar position. We stress that the results here do not conclusively rule out the possibility that a systematic issue produces the peculiar 30–50 GeV spectral dip, as we are unable to adequately account for the extremely unique role of the Sun in both the coordinate system and the operation of the Fermi-LAT. Further studies are necessary on an event-reconstruction level in order to verify the presence and characteristics of the 30-50 GeV spectral feature.
Appendix C Statistical significance of the dip
In this section, we describe our calculation for the significance of the spectral dip between 30–50 GeV.
We first consider the case of the total flux (Figure 3 and Figure 4), using the 9-year averaged flux with the and selection as an example. For this data set, we only consider events GeV, as point-spread-function effect causes the spectrum to depart from the power-law behavior at low energies. This choice is conservative. To estimate the effect of adding the low energy information, we also consider the selection with photons GeV, where the power-law behavior is restored, with is also a reasonably approximation to the SA flux.
In each case, we first fit the data with a broken power-law plus a negative Gaussian line to approximate the dip feature (BPL w/ line). Because the energy, amplitude and width of the line is allowed to vary, the negative Gaussian completely absorbs the two points belonging to the spectral dip, yielding values of nearly zero for the two dip points. Next, we compute the using the obtained power-law parameters, but without the negative line (BPL w/o line). The significance of the line can be obtained by the difference in the value.
In this case, and for the and selection, respectively. The local significance of the dip can then be obtained by using the distribution of two degrees of freedom (width and normalization of the line), which are 5.7 and 6.0 for these two cases. The global significance can be obtained by considering three degrees of freedom (including the energy of the dip), which are 5.4 and 5.7, respectively. For the solar minimum SA flux (Figure 5), we apply the same procedure to the full spectrum above 1 GeV, and obtain , which corresponds to local (global) significance of 8.4 (8.2). These fits are shown in Figure A 10. In the main text, we quote the local significance of the dip with this procedure, with the total flux using the and GeV selection.
However, this definition of the significance is not unique. We can alternatively compare the BPL w/ line with the best-fitting broken power-law (BPL), which is fit to data that includes the dip points. In this case, the best-fit broken power-law adjusts to fit the dip points, reducing the local (global) significance of the dip to 3.1 (2.7) for the 9-year total flux, 2.6 (2.3) for the 9-year total flux, and 5.9 (5.6) for the solar minimum SA flux. Intriguingly, the drop in the significance of the 9-year flux is significantly larger than that for the solar minimum flux. This is due to the fact that the hard solar minimum spectrum prevents the broken-power law from falling quickly near the position of the dip in order to account for the dip points. On the other hand, the 9-year averaged spectrum is softer, allowing the best-fit broken power law to more easily adjust to accommodate the dip.
Appendix D The Width of the Spectral Dip Compared to the Instrumental Energy Dispersion
In the previous section, we have shown that the dip is a statistically significant feature. We now verify that the dip spectrum is not unphysically steep, given the finite energy-resolution of the Fermi-LAT. This is an important diagnostic test, as an extremely steep energy-spectrum would provide a clear indication of instrumental systematic effects. In fact, such a diagnostic was employed in searches of the 130 GeV gamma-ray line to argue for a systematic origin of the line, as the line-feature was more narrow than could be explained by the Fermi-LAT energy dispersion (Ackermann:2013uma, ).
In Figure A 11, we produce three models of a two component solar gamma-ray emission spectrum. The first component is a simple dN/dE E spectrum spanning from 0.1–300 GeV. We then modify this spectrum by placing a “hole” in the gamma-ray emission centered at an energy of 40 GeV. We test three different combinations of hole widths and depths, including a 10 GeV range with 100% of the total gamma-ray emission removed, a 20 GeV range with 50% of the total gamma-ray emission removed, and a 40 GeV range with 25% of the gamma-ray emission removed. In each case, we smear the injected spectrum by the Fermi-LAT energy dispersion on a photon by photon basis, utilizing distributions in both and instrumental energy-dispersion class that are identical to the photons recorded in our main analysis. This allows us to recover the parameters of the observed spectral dip and compare with our observations. We repeat this process several hundred times to build up a large statistical sample of possible dip configurations, and show the average results.
We note two important results. First, our steepest dip (10 GeV range with 100% of the injected emission removed), produces an observed spectral feature which is both sharper and deeper than Fermi-LAT observations. This indicates that the dip is not physically impossible, given the finite accuracy of Fermi-LAT energy reconstructions. Second, we note that the dip is best fit not by an emission component which blocks 100% of the emission in a given energy range, but instead by a component that includes approximately a 50% inhibition of gamma-ray emission in a certain energy range, and which effects solar photons over an energy range of 20 GeV.
Appendix E Examining Event-Level Photon Quality Near the Spectral Dip
In Section VI.1.3, we discussed the possibility that the observed spectral-dip may be due to a contamination of the gamma-ray spectrum by solar neutrons produced in the same hadronic interactions that produce the observed gamma-ray emission. Because neutrons undergo hadronic interactions within the Fermi-LAT detector, the energy-reconstruction of these events may be inaccurate, potentially producing unexpected spectral features.
We note that this possibility could be directly tested by comparing the shower reconstruction of events near the position of the 30-50 GeV spectral dip. Unfortunately, the full event-level shower reconstruction information is not publicly available to those outside the Fermi-LAT collaboration. We note that the re-emergence of the solar spectrum above 50 GeV stands as the most puzzling portion of the solar spectrum, as the sharp falloff near 30 GeV could otherwise be interpreted as as a simple exponential suppression, as predicted in the SSG1991 mechanism. In this scenario, the emission above 50 GeV would be produced predominantly by neutrons with inaccurate energy reconstructions.
In this section, we instead examine several event-level qualities that may be correlated with hadronic interactions within the Fermi-LAT detector. However, we stress here that these qualities are “black boxes”, and we have no pure neutron source to examine in order to determine whether any of these qualities would, in fact, provide reasonable neutron/gamma-ray separation.
In Figure A 12, we examine events with energies between 50–100 GeV, that were recorded after the end of solar minimum (after 2010-01-21). The second cut ensures that no high-energy gamma-rays (observed primarily during solar minimum) were likely to be present in the event analysis. We cut our photon selection between Solar Disk events, recorded within 0.75 of the contemporaneous solar position, and background events, recorded between 0.75 and 5 from the Sun. Because these background events were recorded relatively close to the solar position (but were not produced by the solar disk itself), they will have similar distributions in Fermi-LAT instrumental coordinates, implying that their event reconstructions should be similar.
For each set of events, we examine three possible parameter spaces that might differentiate neutron and gamma-ray interactions within the Fermi-LAT detector. The first (left) is the parameter Tkr1FirstLayer, which specifies the first silicon layer to record the observed event. Because neutron interaction length is much longer than the gamma-ray radiation length in Tungsten, we would expect the majority of neutron events to be observed predominantly in the final layers of the Fermi-LAT instrument. The second (center) is a comparison of the parameters WP8CTAllProb and WP8CTCalTkProb, which compare the probability of cosmic-ray contamination using, or ignoring, information from the Anti-Coincidence detector, respectively. Because the anti-Coincidence detector is the primary rejection technique for charged cosmic-rays, but is unlikely to be triggered by neutrons, we might expect neutron events to have abnormally large cosmic-ray probabilities when information from the Anti-Coincidence detector is ignored. The third (right) is a comparison of the reconstructed energy of the event, compared to the parameter CalEnergyRaw, which denotes the energy deposited in the calorimeter. We may expect any event (neutron or gamma-ray) that has an overestimated energy, to have abnormally large variations between the recorded and reconstructed gamma-ray energies.
However, for each test, we find no statistically significant difference between source and background events, indicating that there are no obvious differences between solar disk events and nearby background regions. However, we again caution the reader that these tests do not definitively rule out the neutron hypothesis, as we cannot calculate the sensitivity of any of these tests to a true neutron signal.
Appendix F Testing dip from modified proton spectrum
In this section, we consider the possibility that the observed gamma-ray dip between 30–50 GeV is produced by two different proton spectra via hadronic interactions. We consider various configurations of proton spectrum that may produce the gamma-ray spectrum, with the proton-proton Kelner:2006tc () interaction optimistically set in the thin regime (interaction length 1).
Ignoring the dip, the gamma-ray spectrum below 100 GeV is roughly , which can be obtained by a input proton spectrum . This is shown by the black line in Figure A 13 (no gap). To test the two component hypothesis, we artificially make a infinite-deep gap in the proton spectrum. The proton gap thus mimics a rapid decay and a rapid rise for the low- and high-energy proton component, respectively. The resultant spectra, corresponding to proton gaps in – GeV, – GeV, and – GeV, are shown in Figure A 13. We find that in all cases, the dip in the computed spectra are too shallow and wide to explain the observed dip feature. This is because the -to- kernel (i.e. gamma-ray spectrum per proton), which is peaked near , is far too wide (see, e.g., Kelner:2006tc ()), making the high-energy component rises slowly.
We also test the case that the high-energy proton spectrum has a harder spectral index. However, we find that the changing the proton shape cannot make gamma-ray dip sharper. This is because the gamma-ray spectral shape just below the cutoff only depends on the shape of -to- kernel, not the shape of proton spectrum.
Lastly, we note that -to- kernel can be sharper if only contributions from a small forward cone is considered. This is because the low energy gamma rays mostly emitted at larger emission angle. However, it is not clear how such a configuration are achieved in a realistic solar atmospheric environment.
From these studies, we find that it is difficult to produce the depth and the width of the observed gamma-ray dip using two different populations of proton, unless taking extreme kinematic regime.
Appendix G Time variations
In the main text, we calculated the time variation of the SA flux, using only data with 30 and two logarithmic energy bins per decade, in order to present the main result with minimal statistical fluctuations. In Figure A 14, we show the same result, instead using four logarithmic bins per decade, and additionally showing the raw (not background subtracted) results for data both in the latitude range 30 and 5. These latter datasets have fewer systematic assumptions (because no modeling is done to remove either diffuse or IC gamma-rays), but show variations with smaller amplitudes (because these background sources are in steady state).
- (1) B. Zhou, K. C. Y. Ng, J. F. Beacom, and A. H. G. Peter, Phys. Rev. D96, 023015 (2017), 1612.02420.
- (2) D. Seckel, T. Stanev, and T. K. Gaisser, Astrophys. J. 382, 652 (1991).
- (3) M. Potgieter, Living Rev. Solar Phys. 10, 3 (2013), 1306.4421.
- (4) T. Wiegelmann, J. K. Thalmann, and S. K. Solanki, A&A Rev. 22, 78 (2014), 1410.4214.
- (5) R. K. Leane, K. C. Y. Ng, and J. F. Beacom, Phys. Rev. D95, 123016 (2017), 1703.04629.
- (6) C. Arina, M. Backović, J. Heisig, and M. Lucente, Phys. Rev. D96, 063010 (2017), 1703.08087.
- (7) E. Orlando and A. Strong, Astrophys. Space Sci. 309, 359 (2007), astro-ph/0607563.
- (8) I. V. Moskalenko, T. A. Porter, and S. W. Digel, Astrophys. J. 652, L65 (2006), astro-ph/0607521, [Erratum: Astrophys. J. 664, L143 (2007)].
- (9) E. Orlando and A. Strong, (2013), 1307.6798.
- (10) E. Kafexhiu, C. Romoli, A. M. Taylor, and F. Aharonian, (2018), 1803.02635.
- (11) Fermi-LAT Collaboration, M. Ajello et al., Astrophys. J. 789, 20 (2014), 1304.5559.
- (12) M. Pesce-Rollins et al., Astrophys. J. 805, L15 (2015), 1505.03480.
- (13) Fermi-LAT, M. Ackermann et al., Astrophys. J. 835, 219 (2017), 1702.00577.
- (14) E. Orlando and A. W. Strong, Astron. Astrophys. 480, 847 (2008), 0801.2178.
- (15) J. F. Dolan and G. G. Fazio, Reviews of Geophysics 3, 319 (1965).
- (16) D. J. Thompson, D. L. Bertsch, D. J. Morris, and R. Mukherjee, Journal of Geophysical Research: Space Physics 102, 14735 (1997).
- (17) Fermi-LAT Collaboration, A. A. Abdo et al., Astrophys. J. 734, 116 (2011), 1104.2093.
- (18) K. C. Y. Ng, J. F. Beacom, A. H. G. Peter, and C. Rott, Phys. Rev. D94, 023004 (2016), 1508.06276.
- (19) Argo-YBJ, G. Aielli et al., Nucl. Instrum. Meth. A562, 92 (2006).
- (20) TIBET ASgamma Collaboration, K. Hibino et al., Nucl. Phys. Proc. Suppl. 10B, 219 (1989).
- (21) A. U. Abeysekara et al., Astropart. Phys. 50-52, 26 (2013), 1306.5800.
- (22) A. U. Abeysekara et al., Astrophys. J. 843, 39 (2017), 1701.01778.
- (23) Z. Cao, Frascati Phys. Ser. 58, 331 (2014).
- (24) LHAASO Collaboration, H. He, PoS ICRC2015, 1010 (2016).
- (25) I. V. Moskalenko, S. Karakula, and W. Tkaczyk, Astron. Astrophys. 248, L5 (1991).
- (26) I. V. Moskalenko and S. Karakula, J. Phys. G19, 1399 (1993).
- (27) G. Ingelman and M. Thunman, Phys. Rev. D54, 4385 (1996), hep-ph/9604288.
- (28) C. Hettlage, K. Mannheim, and J. G. Learned, Astropart. Phys. 13, 45 (2000), astro-ph/9910208.
- (29) G. L. Fogli, E. Lisi, A. Mirizzi, D. Montanino, and P. D. Serpico, Phys. Rev. D74, 093004 (2006), hep-ph/0608321.
- (30) C. A. ArgÃ¼elles, G. de Wasseige, A. Fedynitch, and B. J. P. Jones, JCAP 1707, 024 (2017), 1703.07798.
- (31) K. C. Y. Ng, J. F. Beacom, A. H. G. Peter, and C. Rott, Phys. Rev. D96, 103006 (2017), 1703.10280.
- (32) J. Edsjo, J. Elevant, R. Enberg, and C. Niblaeus, JCAP 1706, 033 (2017), 1704.02892.
- (33) M. Masip, Astropart. Phys. 97, 63 (2018), 1706.01290.
- (34) IceCube, S. In and C. Rott, Solar atmospheric neutrino search with IceCube, in Proceedings, 35th International Cosmic Ray Conference (ICRC 2017): Bexco, Busan, Korea, July 12-20, 2017, 2017, 1710.01194.
- (35) KM3Net, S. Adrian-Martinez et al., J. Phys. G43, 084001 (2016), 1601.07459.
- (36) T. Linden et al., (2018), 1803.05436.
- (38) W. A. Rolke, A. M. Lopez, and J. Conrad, Nucl. Instrum. Meth. A 551, 493 (2005), arXiv:physics/0403059.
- (39) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011), arXiv:1007.1727.
- (40) Fermi-LAT, M. Ackermann et al., Astrophys. J. 799, 86 (2015), 1410.3696.
- (41) E. Orlando and A. W. Strong, Nucl. Phys. Proc. Suppl. 239-240, 266 (2013), 1303.5491.
- (42) http://www.sws.bom.gov.au/Solar/1/6.
- (43) M. Ackermann et al., ApJS 203, 4 (2012), 1206.1896.
- (46) HAWC, M. U. Nisa, Probing Cosmic-ray Propagation with TeV Gamma Rays from the Sun Using the HAWC Observatory, in Proceedings, 35th International Cosmic Ray Conference (ICRC 2017): Bexco, Busan, Korea, July 12-20, 2017, 2017, 1708.03732.
- (47) Argo-YBJ, Z. Li, S. Z. Chen, and H. H. He, 7th Workshop on Air Shower Detection at High Altitude .
- (48) Tibet ASgamma, M. Amenomori et al., Phys. Rev. Lett. 111, 011101 (2013), 1306.3009.
- (49) Particle Data Group, C. Patrignani et al., Chin. Phys. C40, 100001 (2016).
- (50) Y. Muraki et al., Solar Physics 291, 1241 (2016), 1508.04923.
- (51) Y. Muraki et al., (2017), 1706.09082.
- (52) M. D. Kistler and J. F. Beacom, Phys. Rev. D74, 063007 (2006), astro-ph/0607082.
- (53) IceCube, M. G. Aartsen et al., Eur. Phys. J. C77, 146 (2017), 1612.05949.
- (54) E. Bagnaschi et al., (2017), 1710.11091.
- (55) J. C. Costa et al., Eur. Phys. J. C78, 158 (2018), 1711.00458.
- (56) N. F. Bell and K. Petraki, JCAP 1104, 003 (2011), 1102.2958.
- (57) J. L. Feng, J. Smolinsky, and P. Tanedo, Phys. Rev. D93, 115036 (2016), 1602.01465.
- (59) Fermi-LAT, M. Ackermann et al., Phys. Rev. D88, 082002 (2013), 1305.5597.
- (60) S. R. Kelner, F. A. Aharonian, and V. V. Bugayov, Phys. Rev. D74, 034018 (2006), astro-ph/0606058, [Erratum: Phys. Rev. D79, 039901 (2009)].