A SEARCH FOR POINT SOURCES OF EeV PHOTONS

A direct way to identify the origins of cosmic rays is to find fluxes of photons (gamma rays) coming from discrete sources. This method has been used to identify several likely sources of cosmic rays up to about 100 TeV in the Galaxy (Fermi-LAT Collaboration 2013). At sufficiently high energies, such photons must be produced primarily by π⁰ decays, implying the existence of high-energy hadrons that cause the production of π⁰ mesons at or near the source. It is not known whether the Galaxy produces cosmic rays at EeV energies (1 EeV = 10¹⁸ eV). An argument in favor is that the "ankle" of the cosmic-ray energy spectrum near 5 EeV is the only concave upward feature, and the transition from a Galactic power-law behavior to an extragalactic contribution should be recognizable as just such a spectral hardening (Hillas 1984). The ankle can be explained alternatively as a "dip" due to e^± pair production in a cosmic-ray spectrum that is dominated by protons of extragalactic origin. In that case, detectable sources of EeV protons would not be expected in the Galaxy (Berezinsky et al. 2006).

Protons are known to constitute at least a significant fraction of the cosmic rays near the ankle of the energy spectrum (The Pierre Auger Collaboration 2012a, 2010a, 2013). Those protons are able to produce photons with energies near 1 EeV by pion photoproduction or inelastic nuclear collisions near their sources. A source within the Galaxy could then be identified by a flux of photons arriving from a single direction.

The search here is for fluxes of photons with energies from 10^17.3 eV up to 10^18.5 eV. The energy range is chosen to account for high event statistics and to avoid additional shower development processes that may introduce a bias at highest energies (Homola & Risse 2007).

The Pierre Auger Observatory (The Pierre Auger Collaboration 2004) has excellent sensitivity to EeV photon fluxes due to its vast collecting area and its ability to discriminate between photons and hadronic cosmic rays (The Pierre Auger Collaboration 2009). The surface detector array (SD) (The Pierre Auger Collaboration 2008a) consists of 1660 water-Cherenkov detectors spanning 3000 km² on a regular grid of triangular cells with 1500 m spacing between nearest neighbor stations. It is located at latitude −352 in Mendoza Province, Argentina. Besides the surface array, there are 27 telescopes of the air fluorescence detector (FD; The Pierre Auger Collaboration 2010b) located at five sites on the perimeter of the array. The FD is used to measure the longitudinal development of air showers above the surface array. The signals in the water-Cherenkov detectors are used to obtain the secondary particle density at ground measured as a function of distance to the shower core. The analysis presented in this work uses showers measured in hybrid mode (detected by at least one FD telescope and one SD station). The hybrid measurement technique provides a precise geometry and energy determination with a lower energy detection threshold compared to SD only measurements (The Pierre Auger Collaboration 2010b). Moreover, multiple characteristics of photon-induced air showers can be exploited by the two detector systems in combination, e.g., muon-poor ground signal and large depth of shower maximum compared to hadronic cosmic rays of the same energy. Several photon–hadron discriminating observables are defined and combined in a multivariate analysis (MVA) to search for photon point sources and to place directional upper limits on the photon flux over the celestial sphere up to declination +20°.

The sensitivity depends on the declination of a target direction. For the median exposure, a flux of 0.14 photons km⁻² yr⁻¹ or greater would yield an excess of at least 5σ. This corresponds to an energy flux of 0.25 eV cm⁻² s⁻¹ for a photon flux following a 1/E² spectrum, similar to energy fluxes of sources measured by TeV gamma-ray detectors. This is relevant because the energy flux per decade is the same in each energy decade for a source with a 1/E² spectrum, and Fermi acceleration leads naturally to such a type of spectrum (cf. Section 8). The Auger Observatory has the sensitivity to detect photon fluxes from such hypothetical EeV cosmic-ray sources in the Galaxy.

At EeV energies, fluxes of photons are attenuated over intergalactic distances by e^± pair production in collisions of those photons with cosmic-background photons. The e^± can again interact with background photons via inverse-Compton scattering, resulting in an electromagnetic cascade that ends at GeV–TeV energies. The detectable volume of EeV photon sources is small compared to the Greisen–Zatsepin–Kuz'min sphere (Greisen 1966; Zatsepin & Kuz'min 1966), but large enough to encompass the Milky Way, the Local Group of galaxies, and possibly Centaurus A, given an attenuation length of about 4.5 Mpc at EeV energies (Homola & Risse 2007; De Angelis et al. 2013; De Domenico & Settimo 2013).

The present study targets all exposed celestial directions without prejudice. It is a "blind" search to see if there might be an indication of a photon flux from any direction. One or more directions of significance might be identified for special follow-up study with future data. Because there is a multitude of "trials," some excesses are likely to occur by chance. A genuine modest flux would not be detectable in this kind of blind search. The possible production of ultra-high energy photons and neutrons has been studied extensively in relation to some directions in the Galaxy (Medina Tanco & Watson 2001; Bossa et al. 2003; Aharonian & Neronov 2005; Crocker et al. 2005; Gupta 2012).

The Auger Collaboration has published stringent upper limits on the diffuse intensity of photons at ultra-high energies (The Pierre Auger Collaboration 2007, 2008b, 2009; Settimo 2011). Those limits impose severe constraints on "top-down" models for the production of ultra-high energy cosmic rays. At the energies in this study, however, the limits do not preclude photon fluxes of a strength that would be detectable from discrete directions.

The paper is organized as follows. In Section 2, mass composition-sensitive observables are introduced, exploiting information from the SD as well as from the fluorescence telescopes. These observables are combined in a MVA explained in Section 3, before introducing the data set and applied quality cuts in Section 4. A calculation of the expected isotropic background contribution is described in Section 5. A blind search technique and an upper limit calculation are explained in Sections 6 and 7, respectively. Finally, results are shown and discussed in Section 8.

The strategy in searching for directional photon point sources is based on the selection of a subset of photon-like events, to reduce the isotropic hadronic background. Such a selection relies on the combination of several mass composition-sensitive parameters, using an MVA.

Once the MVA training is defined, the photon-like event selection is optimized direction-wise, accounting for the expected background contribution from a given target direction to take into account the contribution of different trigger efficiencies.

Profiting from the hybrid nature of the Auger Observatory, we make use of FD- and SD-based observables, which provide complementary information on the longitudinal and lateral distributions of particles in the showers, respectively. By means of Monte Carlo (MC) simulations, five observables are selected to optimize the signal (photon) selection efficiency against the background (hadron) rejection power. The selected observables are described below in detail.

2.1. FD Observables

2.2. SD Observables

When observed at ground, photon-induced showers have a generally steeper lateral distribution than nuclear primaries because of the almost absent muon component. It is worth noting that, as a consequence of the trigger definition in the local stations and of the station spacing in the array (The Pierre Auger Collaboration 2010c), the SD alone is not fully efficient in the energy range used in this work. Thus, as opposed to previous work based on SD observables (The Pierre Auger Collaboration 2008b), we adopt here observables that are defined at the station level and which do not necessarily require an independent reconstruction in SD mode. Such observables are related to an estimator (S_b ) of the lateral distribution of the signal or to the shape of the flash analog digital converter (FADC) trace in individual stations.

The S_b parameter is sensitive to different lateral distribution functions, due to the presence/absence of the flatter muon component (Ros et al. 2011), and has already been used in previous studies (Settimo 2011). It is defined as

$$\begin{equation} S_b = \sum _{i=1}^N \left[ S_i \cdot \left(\frac{r_i}{1000 \,{\rm m}} \right)^b \right], \end{equation} \tag{ 1 }$$

where the sum extends over all N triggered stations, S_i expresses the signal strength of the ith SD station, r_i is the distance of this station to the shower axis, and b is a variable exponent. It has been found that, in the energy region of interest, the optimized b for photon–hadron separation is b = 3 (Ros et al. 2013).

As a result of both the smaller signal in the stations, on average, and the steeper lateral distribution function, smaller values of S_b are expected for photon primaries. To prevent a possible underestimate of S_b (which would mimic the behavior of a photon-like event), due to missing stations during the deployment of the array or temporarily inefficient stations, events are selected requiring at least four active stations (fully operational, but not necessary triggered) within 2 km from the core.

Other observables, containing information on the fraction of electromagnetic and muonic components at the ground, are related to measurements of the time structure derived from the FADC traces in the SD. The spread of the arrival times of shower particles at a fixed distance from the shower axis increases for smaller production heights, i.e., closer to the detector station. Consequently, a larger spread is expected in case of deep developing primaries (i.e., photons). Here we introduce the shape parameter, defined as the ratio of the early arriving to the late arriving integrated signal as a function of time measured in the water-Cherenkov detector with the strongest signal:

$$\begin{equation} {\rm Shape\,P}(r,\theta) = \frac{S_{\rm early}(r,\theta)}{S_{\rm late}(r,\theta)}. \end{equation} \tag{ 2 }$$

The early signal S_early is defined as the integrated signal over time bins less than a scaled time $t_i^{\rm scaled} \le 0.6$ μs, beginning from the signal start moment. The scaled time varies for different inclination angles θ and distances r to the shower axis and can be expressed as

$$\begin{equation} t_i^{\rm scaled}(r,\theta)=t_i \cdot \frac{r_0}{r} \cdot \frac{1}{c_1 + c_2 \cdot \cos (\theta)}, \end{equation} \tag{ 3 }$$

where t_i is the real time of bin i and r₀ = 1000 m is a reference distance. c₁ = −0.6 and c₂ = 1.9 are scaling parameters to average traces over different inclination angles. Correspondingly, the late signal S_late is the integrated signal over time bins later than $t_i^{\rm scaled} > 0.6$ μs, until signal end.

Table 1. Overall Separation of the Observables Using BDTs (All) and the Remaining Separation if Excluding a Single Observable from the MVA

Observables	Separation
All	0.668
No Sb	0.438
No Xmax	0.599
No $\chi ^2_{\rm Gr} / {\rm dof}$	0.662
No EGr/EGH	0.664
No ShapeP	0.667

Download table as:  ASCIITypeset image

The selected discriminating observables are combined by a MVA technique to enhance and maximize their photon–hadron separation power. In particular, the analysis was developed by using a Boosted Decision Tree (BDT) as classifier (Breiman et al. 1984; Schapire 1990). Several other classifiers were also tested, but the BDT stands out due to the simplicity of the method, where each training step involves a one-dimensional cut optimization, in conjunction with high-performance photon–hadron discrimination (Kuempel 2011). Another advantage of BDTs is that they are insensitive to the inclusion of poorly discriminating variables. The observables were selected from a larger sample of SD and hybrid observables by looking at their individual discrimination power in different energy and zenith bins, their strength in the BDT, and the stability of the results. One benchmark quantity that assess the performance of BDTs is the separation. It is defined to be zero for identical signal and background shapes of the output response, and it is one for shapes with no overlap. The overall separation as well as the separation if excluding a specific observable from the analysis are listed in Table 1. The most significant variables contributing to the BDT are X_max and S_b . With these variables alone we achieve a separation of 0.654. During the classification process the BDT handles the correlation of the observables to energy and zenith angle of the primary particle by including them as additional parameters.

For the classification process, BDTs are trained and tested using MC simulations. Air showers are simulated using the CORSIKA v. 6.900 (Heck et al. 1998) code. A total number of ∼30,000 photon and ∼60,000 proton primaries are generated according to a power-law spectrum of index −2.7 between 10^17.2 eV and 10^18.5 eV, using QGSJET-01c (Kalmykov & Ostapchenko 1989a, 1989b) and GHEISHA (Fesefeldt 1985) as high- and low-energy interaction models, respectively. The impact of different hadronic interaction models is discussed in Section 8. The detector response was simulated using the simulation chain developed within the offline framework (Argiro et al. 2007), as discussed in (Settimo 2012). The same reconstruction chain and the same selection criteria as for data (discussed in Section 4) were then applied. During the classification phase, photon and proton showers are reweighted according to a spectral index of −2.0 and −3.0, respectively. The impact of a changing photon spectral index on the results is discussed in Section 8. The distribution of observables for photon and proton simulations for a specific energy and zenith range is shown in Figure 1. The MVA output response value is named β and shown in Figure 2 for the training and testing samples. Note that the β distribution is by construction limited to the range [ − 1, 1].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The search for photon point sources is performed on the sample of hybrid events collected between 2005 January and 2011 September, under stable data-taking conditions. Events are selected requiring a reconstructed energy between 10^17.3 eV and 10^18.5 eV, where the energy is determined as the calorimetric one plus a 1% missing energy correction associated with photon primaries (Barbosa et al. 2004). To ensure good energy and directional reconstruction, air showers with zenith angle smaller than 60° and with a good reconstruction of the shower geometry are selected. For a reliable profile reconstruction we require: a reduced χ² of the longitudinal profile fit to the Gaisser–Hillas function smaller than 2.5, a Cherenkov light contamination smaller than 50%, and an uncertainty of the reconstructed energy less than 40%. Overcast cloud conditions can distort the light profiles of EAS and influence also the hybrid exposure calculation (Chirinos 2013). To reject misreconstructed profiles, we select only periods with a detected cloud coverage ⩽80% with a cut efficiency of 91%. In addition, only events with a reliable measurement of the vertical optical depth of aerosols are selected (BenZvi et al. 2007). As already mentioned, at least four active stations are required within 2 km of the hybrid-reconstructed axis to prevent an underestimation of S_b . To enrich our sample with deep showers, we do not require that X_max has been observed within the field of view. This cut is usually applied to assure a good X_max resolution (e.g., The Pierre Auger Collaboration 2009, 2010a, 2010d), but for this analysis we focus on a maximization of the acceptance for photon showers. Profiles, for which only the rising edge is observed are certain to have a deep X_max below ground level (∼840/cos (θ) g cm⁻²). Shallow showers, which are also enriched by the release of this cut, can be easily removed in later stages of the analysis by the MVA β cut (cf. Section 6). The photon X_max resolution with this type of selection increases from 39 to 53 g cm⁻², but the photon acceptance is increased by 42%. The energy resolution is about 20%, independently of the primary mass. These resolutions do not affect significantly the analysis since the trace of photons from a point source is an accumulation of events from a specific direction, and the event direction is well reconstructed also with the relaxed cuts: as shown in Figure 3, the angular resolution is about 07. We also verified that the separation power of the MVA is not significantly modified by the weaker selection requirements.

Zoom In Zoom Out Reset image size

After selection, the final data set consists of N_data = 241,466 events with an average energy of 10^17.7 eV. In fact, the energy distribution of these events expresses a compensation effect of the energy spectral index and trigger inefficiencies at low energies. A discussion of the hybrid trigger efficiency for hadrons in the energy range below 10¹⁸ eV is given in (Settimo 2012). In Section 7 this discussion is extended to the case of photons. The average number of triggered stations in the current data set is 2 at 10^17.5 eV, where the bulk of events is detected, generally increasing with zenith angle and with energy (up to 4 between 10¹⁸ eV and 10^18.5 eV).

The contribution of an isotropic background is estimated using the scrambling technique (Cassiday et al. 1990). This method has the advantage of using only measured data and takes naturally into account detector efficiencies and aperture features. Therefore, it is not sensitive to the (unknown) cosmic ray mass composition in the covered energy range.

As a first step, the arrival directions (in local coordinates) of the events are smeared randomly according to their individual reconstruction uncertainty. In a second step, N_data events are formed by choosing randomly a local coordinate and, independently, a Coordinated Universal Time from the pool of measured directions and times. This procedure is repeated 5000 times. The mean number of arrival directions within a target is then used as the expected number for that particular sky location. As each telescope has a different azimuthal trigger probability, events are binned by telescope before scrambling. The number of events observed in each telescope varies between 4358 and 14,100. Since the scrambling technique is less effective in the southern celestial pole region, ¹⁰³ declinations <−85° are omitted from the analysis.

Sky maps are pixelized using the HEALPix software (Gorski et al. 2005). Target centers are taken as the central points of a HEALPix grid using N_side = 256 (target separation ∼03), resulting in 526,200 target centers south of a declination of +20°. The treatment of arrival directions is based on an unbinned analysis, i.e., angular distances are calculated analytically. For each target direction, we use a top-hat counting region of 1° (selecting events within a hard cut on angle from the target center), motivated to account for low event statistics (cf. Alexandreas et al. 1993) and a possible non-Gaussian tail of the error distribution. ¹⁰⁴

The expected directional background contribution for the covered search period is shown in Figure 4. There is an azimuthal asymmetry in the expected background as a result of a seasonally dependent duty cycle, i.e., during austral summer, data taking using the fluorescence telescopes is reduced compared to austral winter (Settimo 2012; The Pierre Auger Collaboration 2011a).

Zoom In Zoom Out Reset image size

When performing the blind search analysis, we use only a subset of the recorded data, selected as "photon-like" according to the β distribution. The definition of "photon-like" (i.e., the β_cut position when selecting events with β ⩾ β_cut) is related to the MC photon and the data efficiencies, $\varepsilon _\gamma ^\beta$ and $\varepsilon _{\rm data}^\beta$, respectively, and to the expected number of background events n_b (α, δ) which is a function of the celestial coordinates α and δ. The efficiencies $\varepsilon _\gamma ^\beta$ and $\varepsilon _{\rm data}^\beta$ are shown in Figure 5 as a function of the multivariate cut β_cut. To estimate $\varepsilon _{\rm data}^\beta$ more accurately, a declination dependence is taken into account, $\varepsilon _{\rm data}^\beta =\varepsilon _{\rm data}^\beta ({\rm \delta })$, indicated by the red shaded area in Figure 5. The expectation of a purely hadronic composition is shown as a gray band. To improve the detection potential of photons from point sources, the cut on the β distribution is optimized, dependent on the direction of a target center or, more specifically, dependent on the expected number of background events n_b (α, δ). In this way the background contamination is reduced while keeping most of the signal events in the data set. This optimization procedure can be described as follows: the upper limit of photons n_s from a point source at a given direction is calculated under the assumption that $n_{\rm data} = n_b^{\beta }$, i.e., when the observed number of events (n_data) is equal to the expected number (cf. Section 5). The expected number of events after cutting on the β distribution can be estimated as $n_b^\beta (\alpha, \delta) = n_b(\alpha, \delta) \cdot \varepsilon _{\rm data}^\beta (\delta)$ and is typically less than four events for the β_cut values finally chosen. There are several ways to define an upper limit to the number of photons n_s , at a given confidence level (CL) in the presence of a Poisson-distributed background. Here the procedure of Zech (Zech 1989) is utilized, where n_s is given by

$$\begin{equation} P\left(\le n_b^\beta | n_b^\beta + n_s\right) = \alpha _{\rm CL} \cdot P\left(\le n_b^\beta | n_b^\beta\right), \end{equation} \tag{ 4 }$$

with α_CL ≡ 1 − CL = 0.05, and where the expected background contribution is $n_b^\beta$ (cf. The Pierre Auger Collaboration 2012b). The frequentist interpretation of the above equation is as follows: for an infinitely large number of repeated experiments looking for a signal with expectation n_s and background with mean $n_b^\beta$, where the background contribution is restricted to a value less than or equal to $n_b^\beta$, the frequency of observing $n_b^\beta$ or fewer events is α_CL. Since $n_b^\beta$ is not an integer in general, a linear interpolation is applied to calculate the Poisson expectation. To determine the optimized β_cut, the sensitivity is maximized by minimizing the expected upper limit by scanning over the entire range of possible β_cut, also taking into account the photon efficiency $\varepsilon _\gamma ^\beta$:

$$\begin{equation} \min \left(\frac{n_s(\beta _{\rm cut})}{\varepsilon _\gamma ^\beta (\beta _{\rm cut})} \right)\ {\rm with}\ \beta _{\rm cut} \in [-1,1]. \end{equation} \tag{ 5 }$$

The optimized mean β_cut is shown in Figure 6 as a function of the expected background contribution. The gray area indicates the declination-dependent variation of the optimization. The mean β_cut value used in this analysis is 0.22 resulting in an average background contribution after β_cut of 1.48 events. Applying the optimized β_cut to measured data reduces the data set to 13,304 events. The sky distribution of these events is shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When performing a blind search for photon point sources, the probability p of obtaining a test statistic at least as extreme as the one that was actually observed is calculated, assuming an isotropic distribution. The test statistic is obtained from the ensemble of scrambled data sets (cf. Section 5), assuming a Poisson-distributed background. This p value is calculated for a specific target direction as

$$\begin{equation} p = {\rm Poiss} \left({\ge} n_{\rm data}^\beta | n_b^\beta\right), \end{equation} \tag{ 6 }$$

where ${\rm Poiss} (\ge n_{\rm data}^\beta | n_b^\beta)$ is the Poisson probability to observe $n_{\rm data}^\beta$ or more events given a background expectation after β_cut of $n_b^\beta$. Note that the superscript "β" indicates the number of events after applying the optimized β_cut. The fraction of simulated data sets p_chance, in which the observed minimum p value p_min is larger than or equal to the simulated p value $p_{\rm min}^{\rm scr}$, is given by

$$\begin{equation} p_{\rm chance}\left(p_{\rm min}^{\rm scr} \le p_{\rm min}\right). \end{equation} \tag{ 7 }$$

This corresponds to the chance probability of observing p_min anywhere in the sky. The results when applying this blind search to the hybrid data of the Pierre Auger Observatory will be discussed in Section 8.

Here we specify the method used to derive a skymap of upper limits to the photon flux of point sources. The directional upper limit on the photon flux from a point source is the limit on the number of photons from a given direction, divided by the directional acceptance (cf. Section 5) from the same target at a confidence level of CL = 95%, and by a correction term:

$$\begin{equation} f^{\rm UL} = \frac{n_s^{\rm Zech}}{n_{\rm inc} \cdot \mathcal {E}_{\rm {\beta }}}. \end{equation} \tag{ 8 }$$

Here $n_s^{\rm Zech}$ is the upper limit on the number of photons obtained by using the β_cut definition in Figure 6, and applying the procedure of Zech (cf. Equation (4)) for the observed number of events in data $n_{\rm data}^\beta$:

$$\begin{equation} P\left({\le} n_{\rm data}^\beta | n_b^\beta + n_s^{\rm Zech}\right) = \alpha _{\rm CL} \cdot P\left({\le} n_{\rm data}^\beta | n_b^\beta\right). \end{equation} \tag{ 9 }$$

The expected signal fraction in the top-hat search region is n_inc = 0.9, and $\mathcal {E_\beta }$ is the total photon exposure. This latter exposure is derived as

$$\begin{equation} \mathcal {E}_{\beta }(\alpha, \delta) = \mathcal {E}(\alpha, \delta)\cdot \varepsilon _\gamma ^\beta, \end{equation} \tag{ 10 }$$

where $\mathcal {E}$ indicates the exposure before applying the multivariate cut β_cut (cf. Equation (11)), and $\varepsilon _\gamma ^\beta$ is the photon efficiency when applying a β_cut.

The exposure $\mathcal {E}(E)$ is typically defined as a function of energy E, cf. (The Pierre Auger Collaboration 2011a; Settimo 2012; The Pierre Auger Collaboration 2010d). In a similar way, the photon exposure $\mathcal {E}$ as a function of celestial coordinates α and δ is defined as

$$\begin{equation} \mathcal {E}(\alpha, \delta) = \frac{1}{c_E} \int _{E}\int _T\int _{S} E^\zeta \varepsilon (E,t,\theta,\phi,x,y) {{d}}S {{d}}t {{d}}E, \end{equation} \tag{ 11 }$$

where the coordinates α and δ are functions of the zenith (θ) and azimuth (ϕ) angles and of the time t; ε is the overall efficiency including detection, reconstruction and selection of the events and the evolution of the detector in the time period T. The integration over energy is performed assuming a power-law spectrum with index ζ = −2 and normalization factor c_E = ∫E^ζ dE. The area S encloses the full detector array and is chosen sufficiently large to ensure a negligible (less than 1%) trigger efficiency outside of it. The exposure for the hybrid detector is not constant with energy and is not uniform in right ascension. Thus, detailed simulations were performed to take into account the status of the detector and the dependence of its performance with energy and direction (both zenith and azimuth). For the exposure calculation applied here, time-dependent simulations were performed, following the approach described in (The Pierre Auger Collaboration 2011a; Settimo 2012). This takes into account the photon trigger efficiency, possible periods of overcast cloud conditions, and offers the possibility of also investigating systematic uncertainties below the EeV range. At the low energy edge of 10^17.3 eV, the hybrid trigger efficiency for photon-induced air showers is larger than 80%, rapidly increasing before reaching full efficiency at 10^17.8 eV. Compared to hadron induced air showers, the photon trigger efficiency is always larger as a consequence of a later shower development, on average. The derived directional photon exposure has a mean of 180 km² yr and varies between 50 km² yr and 294 km² yr. The impact of different photon spectral indices is discussed in Section 8. Directional upper limits on photons from point sources derived in this blind search analysis will be also discussed in Section 8.

In the following, results on p values and upper limits are given, based on the analysis method described in Sections 6 and 7.

The p values, as defined in Equation (6), refer to a local probability that the data is in agreement with a uniform distribution. The integral distribution of −log (p) values is shown in Figure 8. The corresponding sky map of −log (p) values is illustrated in Figure 9. The minimum p value observed is p_min = 4.5 × 10⁻⁶ corresponding to a chance probability that p_min is observed anywhere in the sky of p_chance = 36%. This blind search for a flux of photons, using hybrid data of the Pierre Auger Observatory, therefore finds no candidate point on the pixelized sky that stands out among the large number of trials. It is possible that some genuine photon fluxes are responsible for some of the low p values. If so, additional exposure should increase the significance of those excesses. They might also be identified in a future search targeting a limited number of astrophysical candidates. The present search, however, finds no statistical evidence for any photon flux.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Directional photon flux upper limits (95% confidence level) are derived using Equation (8) and shown as a celestial map in Figure 10. The mean value is 0.035 photons km⁻² yr⁻¹, with a maximum of 0.14 photons km⁻² yr⁻¹. Those values correspond to an energy flux of 0.06 eV cm⁻² s⁻¹ and 0.25 eV cm⁻² s⁻¹, respectively, assuming an E⁻² energy spectrum.

Zoom In Zoom Out Reset image size

Various sources of systematic uncertainties were investigated and their impact on the mean flux upper limit is estimated. The systematics on the photon exposure ranges between ±30% at 10^17.3 eV and ±10% above 10¹⁸ eV and are dominated by the uncertainty on the Auger energy scale. A systematic uncertainty of the Auger energy scale of +14% and −14% (Verzi 2013) changes the mean upper limit by about +8% and −9%, respectively. Variations in determining the fraction of photon ($\epsilon ^\beta _\gamma$) and measured events ($\epsilon ^\beta _{\rm data}$) passing a β_cut, introduced by, e.g., an additional directional dependency for photons, contribute less than 6%. A collection of ∼50,000 proton CORSIKA air shower simulations, using EPOS Large Hadron Collider (Pierog et al. 2013), were additionally generated to estimate the impact of using a different high-energy hadronic interaction model. The resulting change of the mean limit of the photon flux is −9%. Furthermore, the assumed photon flux spectral index of −2 could be incorrect. To estimate the impact of this, the analysis is repeated assuming a spectral index of −1.5 or −2.5. The mean upper limit changes by about −34% and +51%, respectively, whereas the dominant contribution arises from a changing directional photon exposure, i.e., assuming a flatter primary photon spectrum increases the photon exposure, while reducing the average upper limits, and vice versa.

The limits are of considerable astrophysical interest in all parts of the exposed sky. The energy flux in TeV gamma rays exceeds 1 eV cm⁻² s⁻¹ for some Galactic sources with a differential spectral index of E⁻² (Hinton & Hofmann 2009; H. E. S. S. Collaboration 2011). A source with a differential spectral index of E⁻² puts out equal energy in each decade, resulting in an expected energy flux of 1 eV cm⁻² s⁻¹ in the EeV decade. No energy flux that strong in EeV photons is observed from any target direction, including directions of TeV sources such as Centaurus A or the Galactic center region. This flux would have been detected with >5σ significance, even after penalizing for the large number of trials (using Equations (6) and (7)). Furthermore, an energy flux of 0.25 eV cm⁻² s⁻¹ would yield an excess of at least 5σ for median exposure targets. If we make the conservative assumption that all detected photons are at the upper energy bound, a flux of 1.44 eV cm⁻² s⁻¹ would be detectable. This result for median exposure targets is independent of the assumed photon spectral index, and implies that we can exclude a photon flux greater than 1.44 eV cm⁻² s⁻¹ with 5σ significance.

Results from the present study complement the blind search for fluxes of neutrons above 1 EeV previously published by the Auger Collaboration (The Pierre Auger Collaboration 2012b). No detectable flux was found in that search, and upper limits were derived for all directions south of declination +20°. A future study will look for evidence of photon fluxes from particular candidate sources and "stacks" of candidates having astrophysical characteristics in common. A modest excess may be statistically significant if it is not penalized for a large number of trials. Neutrons and photons arise from the same types of pion-producing interactions. The photon path length exceeds the path length for EeV neutron decay, so this study is sensitive to sources in a larger volume than just the Galaxy.

The absence of detectable point sources of EeV neutral particles does not mean that the sources of EeV rays are extragalactic. It might be that EeV cosmic rays are produced by transient sources such as gamma ray bursts or supernovae. The Auger Observatory has been collecting data only since 2004. It is quite possible that it has not been exposed to neutral particles emanating from any burst of cosmic-ray production. Alternatively, it is conceivable that there are continuous sources in the Galaxy which emit in jets and are relatively few in number, and if so none of those jets are directed toward Earth. The protons would be almost isotropized by magnetic fields, but neutrons and gamma rays would retain the jet directions and would not arrive here. Another possibility is that the EeV protons originate in sources with much lower optical depth for escaping than is typical of the known TeV sources. The production of neutrons and photons at the source could be too meager to make a detectable flux at Earth.

The null results from this search for point sources of photons is nevertheless interesting in light of the stringent upper limits on cosmic-ray anisotropy at EeV energies (The Pierre Auger Collaboration 2012c, 2012d, 2011b). EeV protons originating near the Galactic plane, whether from transient sources or steady sources, are expected to cause an anisotropy that exceeds the observational upper limit. Those expectations are not free of assumptions about magnetic field properties away from the Galactic disk, so the case against Galactic EeV proton sources is by no means closed. However, evidence for such sources remains absent, despite a sensitive search for any flux of EeV photons.

The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe.

We are very grateful to the following agencies and organizations for financial support: Comisión Nacional de Energía Atómica, Fundación Antorchas, Gobierno De La Provincia de Mendoza, Municipalidad de Malargüe, NDM Holdings and Valle Las Leñas, in gratitude for their continuing cooperation over land access, Argentina; the Australian Research Council; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Financiadora de Estudos e Projetos (FINEP), Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ), São Paulo Research Foundation (FAPESP) grant Nos. 2010/07359-6, 1999/05404-3, Ministério de Ciência e Tecnologia (MCT), Brazil; MSMT-CR LG13007, 7AMB14AR005, CZ.1.05/2.1.00/03.0058 and the Czech Science Foundation grant 14-17501S, Czech Republic; Centre de Calcul IN2P3/CNRS, Centre National de la Recherche Scientifique (CNRS), Conseil Régional Ile-de-France, Département Physique Nucléaire et Corpusculaire (PNC-IN2P3/CNRS), Département Sciences de l'Univers (SDU-INSU/CNRS), France; Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Finanzministerium Baden-Württemberg, Helmholtz-Gemeinschaft Deutscher Forschungszentren (HGF), Ministerium für Wissenschaft und Forschung, Nordrhein Westfalen, Ministerium für Wissenschaft, Forschung und Kunst, Baden-Württemberg, Germany; Istituto Nazionale di Fisica Nucleare (INFN), Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR), Gran Sasso Center for Astroparticle Physics (CFA), CETEMPS Center of Excellence, Italy; Consejo Nacional de Ciencia y Tecnología (CONACYT), Mexico; Ministerie van Onderwijs, Cultuur en Wetenschap, Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), Stichting voor Fundamenteel Onderzoek der Materie (FOM), Netherlands; National Centre for Research and Development, grant Nos. ERA-NET-ASPERA/01/11 and ERA-NET-ASPERA/02/11, National Science Centre, grant Nos. 2013/08/M/ST9/00322 and 2013/08/M/ST9/00728, Poland; Portuguese national funds and FEDER funds within COMPETE - Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia, Portugal; Romanian Authority for Scientific Research ANCS, CNDI-UEFISCDI partnership projects nr.20/2012 and nr.194/2012, project nr.1/ASPERA2/2012 ERA-NET, PN-II-RU-PD-2011-3-0145-17, and PN-II-RU-PD-2011-3-0062, the Minister of National Education, Programme for research—Space Technology and Advanced Research—STAR, project number 83/2013, Romania; Slovenian Research Agency, Slovenia; Comunidad de Madrid, FEDER funds, Ministerio de Educación y Ciencia, Xunta de Galicia, Spain; The Leverhulme Foundation, Science and Technology Facilities Council, UK; Department of Energy, Contract No. DE-AC02-07CH11359, DE-FR02-04ER41300, and DE-FG02-99ER41107, National Science Foundation, grant No. 0450696, The Grainger Foundation, USA; NAFOSTED, Vietnam; Marie Curie-IRSES/EPLANET, European Particle Physics Latin American Network, European Union 7th Framework Program, grant No. PIRSES-2009-GA-246806; and UNESCO.