Searches for Ultra-High-Energy Photons at the Pierre Auger Observatory

The Pierre Auger Observatory, being the largest air-shower experiment in the world, offers an unprecedented exposure to neutral particles at the highest energies. Since the start of data taking more than 18 years ago, various searches for ultra-high-energy (UHE, $E\gtrsim10^{17}\,\text{eV}$) photons have been performed: either for a diffuse flux of UHE photons, for point sources of UHE photons or for UHE photons associated with transient events like gravitational wave events. In the present paper, we summarize these searches and review the current results obtained using the wealth of data collected by the Pierre Auger Observatory.


Introduction
The search for neutral particles-in particular photons and neutrinos-of cosmic origin at the highest energies has been for many years one of the major scientific objectives of the Pierre Auger Observatory. From the theory side, such searches are well-motivated: many models for the origin of ultra-high-energy (UHE) cosmic rays predict at least some neutral particles as by-products, either directly at the sources or during the propagation through the Universe (see, e.g., [1][2][3][4][5]). In fact, even though no UHE photons have been unambiguously identified so far, the upper limits on their incoming flux have been already used to severely constrain so-called top-down models for the origin of UHE cosmic rays involving e.g. topological defects or super-heavy dark matter (see, e.g., [6][7][8][9]). In addition, the recent observations of photons with energies up to 10 15 eV [10] further motivate searches for photons at even higher energies. An observation of such photons would also be key in completing the multi-messenger approach aimed at understanding the most extreme processes in the Universe, taking advantage of the fact that neutral particles directly point back at their production site. However, one has to take into account that UHE photons, unlike neutrinos, interact with the background photon fields permeating the Universe, reducing their attenuation length to about 30 kpc around 10 15 eV, which increases to the order of 10 Mpc around 10 19 eV [11]. In the present paper, we review the current state of such searches at the Pierre Auger Observatory. After a short introduction addressing the specificities of air showers initiated by photons (Sec. 2), we briefly describe the Pierre Auger Observatory (Sec. 3). We then focus first on the searches for a diffuse flux of UHE photons using the different detector systems of the Observatory (Sec. 4), before we describe the searches for UHE photons from point sources and transient events (Sec. 5). We close with a short outlook (Sec. 6), discussing the ongoing detector upgrade of the Pierre Auger Observatory, dubbed AugerPrime.

Photon-induced air showers
When a UHE photon enters the Earth's atmosphere, it may interact with a particle from the atmosphere, for example a nitrogen nucleus, and induce an extensive air shower, much in the same way as a charged cosmic ray does. Hence, a cosmic-ray observatory detecting air showers is also, by construction, a photon observatory-and even a neutrino observatory, highlighting the importance of such observatories for multimessenger astrophysics. In fact, since the incoming flux of UHE cosmic particles is so low (on the order of one particle per square kilometer and year and less), measuring the extensive air showers they initiate in the atmosphere with large detector arrays on ground is the only way to efficiently detect them. The challenge lies in distinguishing air showers induced by photons from the vast background of air showers that are initiated by charged cosmic rays, i.e., protons and heavier nuclei (for a review, see, e.g., [11]). The two main differences between photon-and nucleus-induced air showers are shown schematically in Fig. 1. The longitudinal development of an air shower, as a function of the slant depth X, is delayed for a primary photon with respect to primary nuclei, due to the lower multiplicity of electromagnetic interactions (compared to hadronic interactions) that dominate in a photon-induced air shower. The maximum of the shower development within the atmosphere, Xmax, is reached later. For example, at a primary energy of 10 19 eV, the difference is about 200 g cm −2 . Since the mean free path for photo-nuclear interactions is much larger than the radiation length, only a small fraction of the electromagnetic component in a photon-induced shower is transferred to the hadronic component and subsequently to the muonic component. Showers induced by photons are thus characterized by a lower number of muons. On average, simulations show that photoninduced showers have nearly one order of magnitude less muons than those initiated by protons or nuclei of the same primary energy. In one way or another, all searches for UHE photons using air-shower data exploit these two key differences: Xmax, for example, can be directly measured using the air-fluorescence technique. The number of muons cannot yet be directly measured using the current detector systems of the Pierre Auger Observatory. However, one can measure the lateral distribution of secondary particles from the air shower at ground level, which depends on both the number of muons and the longitudinal development. In particular, the steepness of the lateral distribution is sensitive to the type of the primary particle.
The difference in Xmax between photon-and proton-induced air showers is amplified by the Landau-Pomeranchuk-Migdal (LPM) effect [12,13], i.e., the suppression of the bremsstrahlung and pair production cross sections at high energies. In addition, the preshower effect [14][15][16] also has to be taken into account: depending on its incident direction with respect to the local geomagnetic field, a UHE photon can initiate an electromagnetic cascade in the Earth's magnetic field even before entering the atmosphere-a preshower. The observed air shower is then a superposition of the individual cascades of the photons and electrons/positrons from the preshower, leading on average to a smaller measured Xmax than for non-preshowering photons of the same primary energy.

The Pierre Auger Observatory
The Pierre Auger Observatory [17], located near the town of Malargüe in the Argentinian Pampa Amarilla, is the largest cosmic-ray observatory to date, offering an unprecedented exposure to UHE photons. A key feature of the Pierre Auger Observatory is the hybrid concept, combining a Surface Detector array (SD) with a Fluorescence Detector (FD), see Fig. 2. The SD consists of 1600 water-Cherenkov detectors arranged on a triangular grid with a spacing of 1500 m, covering a total area of about 3000 km 2 . The SD is overlooked by 24 fluorescence telescopes, located at four sites at the border of the array. The SD samples the lateral shower profile at ground level, i.e., the distribution of particles as a function of the distance from the shower axis, with a duty cycle of ∼100 %, while the FD records the longitudinal shower development in the atmosphere above the SD. The FD can only be operated in clear, moonless nights, reducing the duty cycle to ∼15 %. Through combining measurements from both detector systems in hybrid events, a superior accuracy of the air-shower reconstruction can be achieved than with just one system [18]. In the western part of the SD array, 50 additional SD stations have been placed between the existing SD stations, forming a sub-array with a spacing of 750 m and covering a total area of about 27.5 km 2 . With this sub-array, air showers of lower primary energy (below 10 18 eV) with a smaller footprint on ground can be measured. To allow also for hybrid measurements in this energy range, where air showers develop above the field of view of the standard FD telescopes, three additional High-Elevation Auger Telescopes (HEAT) have been installed at the FD site Coihueco, overlooking the 750 m SD array. The HEAT telescopes operate in the range of elevation angles from 30 • to 60 • , complementing the Coihueco telescopes operating in the 0 • to 30 • range. The combination Figure 1: Schematic depiction of the main differences between photon-induced air showers and those initiated by primary nuclei (protons or heavier nuclei).
of the data from both HEAT and Coihueco ("HeCo" data) enables fluorescence measurements of air showers over a large range of elevation angles.

Searches for a diffuse flux of UHE photons
First, we focus on the searches for a diffuse-i.e., direction-independent, unresolved-flux of photons. At the Pierre Auger Observatory, such searches have been performed in the past using hybrid data (with the analysis based on Xmax only [19,20], or based on a combination of Xmax and additional SD-related quantities [21]) as well as SD-only data [22]. In the following, we briefly summarize the three most up-to-date publications [23][24][25] in different energy ranges and the corresponding results.
4.1 A search for photons with energies above 2×10 17 eV using hybrid data from the low-energy extensions of the Pierre Auger Observatory We first discuss the photon search targeting the lowest energy range, which starts from 2×10 17 eV [23]. In this energy range, data from the low-energy extensions of the Pierre Auger Observatory, i.e., from the 750 m SD array combined with HeCo data, can be used to efficiently search for photons. Three observables are used in the analysis: Xmax, measured directly with the fluorescence telescopes, is used together with the SD quantities S b and Nstations. S b [26] is defined through where Si denotes the measured signal in the i-th SD station at a perpendicular distance Ri to the shower axis and the parameter b has been chosen as b = 4 to optimize the photon-hadron separation. By construction, S b is sensitive to the lateral distribution, which in turn depends on the depth of the air-shower development in the atmosphere and the number of muons, as stated in Sec. 2. The third observable, Nstations is the number of triggered SD stations, as it has been shown previously that it can significantly improve the overall performance of the analysis [21].
These three quantities are combined in a multivariate analysis (MVA) using the Boosted Decision Tree (BDT) method. To take into account energy and zenith angle dependencies, the photon energy Eγ (defined as the calorimetric energy obtained through the integration of the longitudinal profile plus a missing-energy correction of 1 % appropriate for primary photons [21]) and the reconstructed zenith angle θ are also included  [17]; each dot represents one SD station; also shown are the four FD sites at the border of the SD array. Top right: the fluorescence telescopes at the FD site Los Leones; even though the picture was taken during the daytime, the shutters were had been opened for maintenance. Bottom right: a single SD station in the Pampa Amarilla.
in the MVA. A large sample of simulated events has been used to study the photon/hadron separation of the observables mentioned before, to train the multivariate analysis, and to evaluate its performance. For these samples, both primary photons (as a signal sample) and primary protons (as a conservative, "worst-case" assumption for the background sample) have been simulated using CORSIKA and the Auger Offline Software Framework. The analysis is eventually applied to hybrid data collected by the Coihueco and HEAT telescopes and the 750 m SD array between 1 June 2010 and 31 December 2015, in total more than 500,000 events. A number of selection criteria is applied to both the data sample and the simulated samples to select only well-reconstructed, reliable events. These selection criteria are described in detail in [23]. After all criteria have been applied, 2,204 events remain in the data sample with a photon energy Eγ above 2×10 17 eV.
In Fig. 3, the normalized distributions of the discriminating observables Xmax, S b and Nstations are shown for the simulated samples as well as the data sample. In addition, the corresponding distributions of the output from the BDT β, which is used as the final discriminator for separating photon-induced air showers from the hadronic background, are displayed. A more detailed discussion of these distributions can be found in [23]. Here, we only note that for β, the photon and proton distributions are well-separated. The background rejection at a signal efficiency of 50 %, i.e., the fraction of proton-induced events with β larger than the median of the photon distribution-which is used as the photon candidate cut, marked with the dashed line-is (99.91 ± 0.03) % for energies Eγ ≥ 2×10 17 eV.
Zero events from the data sample have a β value above the candidate cut value, hence no photon candidate events are identified in this analysis. The final results of this study are therefore given in terms of upper limits on the integral flux of photons Φ C.L. γ, U.L. (Eγ>E0). The integrated, efficiency-weighted exposure for photons needed to calculate these upper limits is obtained from simulations. In the energy range of interest between 2×10 17 eV and 10 18 eV, the weighted exposure varies between 2.4 and 2.7 km 2 sr yr under the assumption of a power-law spectrum ∝ E −2 . Upper limits on the integral photon flux are placed at threshold energies of 2, 3, and 5×10 17 eV, as well as 10 18 eV, at a confidence level of 95 %. At these threshold energies, the upper limits are 2.72, 2.50, 2.74, and 3.55 km −2 sr −1 yr −1 , respectively. Using the energy spectrum of cosmic rays measured by the Pierre Auger Observatory [27], the upper limits on the integral photon flux can be translated into upper limits on the integral photon fraction. At a confidence level of 95 %, these are 0.28 %, 0.63 %, 2.20 % and 13.8 % for the same threshold energies as above.

A search for ultra-high-energy photons at the Pierre Auger Observatory exploiting air-shower universality
The photon search in the energy range between 10 18 and 10 19 eV [24] is performed by exploiting the hybrid configuration of the Pierre Auger Observatory, much like in the previous analysis. As before, Xmax can  Figure 3: Normalized distributions of the three discriminating observables X max , S b and N stations used in the photon search based on HeCo data. The (simulated) photon sample is shown in blue, the (simulated) proton sample in red, and the data sample in black. In addition, normalized distributions of the final discriminator β-which is based on a MVA combining X max , S b and N stations as well as the photon energy and the zenith angle-are displayed. The dashed line denotes the median of the photon test sample, which is used as the photon candidate cut. In all plots, only events with E γ >2×10 17 eV are shown. For more details, see [23].
be measured directly with the FD. The muon content of a measured air shower is accessed through the parameter Fµ, which is derived from the SD signals by using the air-shower universality concept [28]. A universality-based model [29] is used to predict the signals induced by the secondary particles in the individual SD stations. This model describes the total signal as the sum of four components: two electromagnetic components, one related to high-energy pions (Seγ) and one related to low-energy hadrons (S eγ(had) ), and two muon-related components, a pure muon component (Sµ) and one related to electrons and photons resulting from muon decays (S eγ(µ) ). The predicted signal, S pred , can then be expressed as where i runs over the four components. Each of the four signal components, S i comp , has a universal behavior, which can be parameterized as a function of the primary energy E, Xmax, and the geometry of the air shower. The relative contributions β i of each of the four components depend only on the mass of the primary particle through a parameter, Fµ, representing the number of muons in the air shower.
Following the approach developed in [30], the universality-based signal model is applied to the case of hybrid events measured by the FD and the SD. As the hybrid reconstruction provides E, Xmax and the shower geometry, the four components S i comp can be directly calculated for each SD station involved in a hybrid event. Thus, given the reconstructed signal, Srec, in a station of the SD, Fµ can be calculated for each station in an event by matching Srec to S pred from Eq. (2). To obtain an event-wise estimate of the muon content, the average Fµ of all SD stations is assigned to an event if more than one station is available.
Overall, Fµ provides a very good photon-hadron separation (see Fig. 4, top left). To fully exploit the hybrid approach, it is combined in this analysis with Xmax in a linear Fisher discriminant analysis [31]. of X max and F µ , i.e., the observables used in the hybrid search for photons using air-shower universality, for simulated primary photons (blue) and protons (red); The contour lines enclose 90 %, 50 % and 10 %, respectively, of the events. Top right: distributions of the Fisher discriminant f for simulated primary photons (signal, blue) and protons (background, red), and for the burnt sample (black); the dashed red line marks the tail of the proton distribution, the dashed blue line indicates the median of the photon distribution. Bottom: the tail of the distribution of f for the hybrid data sample (black dots); the dashed line represents the photon-candidate cut; the shaded blue regions show the 1σ, 2σ and 3σ uncertainty bands for background expectation. For more details, see [24].
The resulting distributions of the Fisher discriminant f are shown in the top right panel of Fig. 4, for simulations of primary protons (red) and photons (blue), as well as for 5 % of the data sample (black), which is used as a burnt sample. The distributions of f obtained for the proton and photon samples are well separated, resulting in a background rejection of about 99.90 % at a signal efficiency of 50 %, corresponding to the dashed blue line in the top right panel of Fig. 4. The burnt sample and the photon distributions are even more separated, therefore the events contained in the burnt sample can be considered as background events only, and can then be used to derive a data-driven estimate of the expected background. Due to the limited number of events in the burnt sample, we used in a first step the rightmost tail of the proton distribution, specifically only the events with f > −1.3 (red dashed line in Fig. 4, top right), to derive the functional form of the background distribution. Then, the expected background distribution was derived through a fit to the burnt sample and scaled up linearly with time to match the full data sample.
Finally, the analysis has been applied to hybrid events above 10 18 eV recorded by the Pierre Auger Observatory between 1 January 2005 and 31 December 2017. Of the total data set, which consists of approximately 32,000 events, the rightmost tail of the distribution of the Fisher discriminant f (black dots) is shown in the bottom panel of Fig. 4. The data distribution resembles the background expectation, shown as the shaded blue bands, representing the uncertainties in its estimation for different σ levels. After applying the photon selection cut (dashed vertical line), corresponding to the median of the photon distribution (dashed blue line in Fig. 4, top right), 22 photon-candidate events are selected. This number is consistent with the background expectation of 30 ± 15 false-positive candidate events. Since no significant excess with respect to the background has been found, the final results of this study are given in terms of upper limits on the integral flux of photons Φ C.L.
γ, U.L. (Eγ>E0). Five different energy thresholds (1, 2, 3, and 5×10 18 eV, as well as 10 19 eV) were considered. The number of photon-candidate events found for each energy threshold were 22, 2, 0, 0 and 0, respectively. The upper limits are determined taking into account the expected number of background events derived from the burnt sample for each threshold (30 ± 15, 6 ± 6, 0.7 ± 1.9, 0.06 ± 0.25 and 0.02 ± 0.06, respectively) as well as the integrated, efficiency-weighted exposure for photons, which was again determined from simulations. In the energy range between 10 18 and 10 19 eV, the weighted exposure increases from 420.7 to 1245.9 km 2 sr yr, under the assumption of a power-law spectrum In the energy range above 10 19 eV, UHE photons are searched for among the data collected with the 1500 m SD array of the Pierre Auger Observatory [25]. While the photon search using SD-only data can profit from the large exposure due to the high duty cycle of the SD, the lack of a corresponding fluorescence measurement for the bulk of the data poses some challenges. For example, there is no direct measurement of Xmax available. Also the primary energy can only be accessed indirectly, using S(1000)-the interpolated signal in the SD stations at a perpendicular distance of 1000 m from the shower axis-as a proxy.
Two observables are used in this analysis, one related to the thickness of the shower front at ground and one based on the steepness of the lateral distribution. These two properties of an air shower depend on the type of the primary particle initiating the shower, hence they can be used for photon-hadron separation. The first observable, ∆, is based on the risetime t 1/2 in the individual SD stations, which is defined as the time at which the integrated signal in the measured time trace rises from 10 % to 50 % of its total value. For showers of the same primary energy and zenith angle, t 1/2 is expected to be larger for primary photons with respect to primary nuclei due to the reduced muonic content, which implies larger scattering and attenuation of secondary particles, and to Xmax being closer to the ground. ∆ is defined as: which can be taken as the average deviation of the measured rise-times from a data benchmark t bench 1/2 , describing the average risetime of all of the SD data (assumed to be overwhelmingly constituted by primary nuclei) [32], in units of sampling fluctuations σt 1/2 . Details on the selection criteria for the SD stations can be found in [25]. By construction, ∆ is expected to average to zero for data and to be significantly positive for photon-induced air showers. As photon-induced air showers are also expected to have a steeper lateral distribution of the signals in the SD stations than the average of all SD data, a second observable LLDF is introduced to quantify the departure of the observed lateral distribution function (LDF) from the average of all SD data (see also [33]): where Si is the total signal measured in the i-th selected station and fLDF(ri) gives the average signal, obtained from all SD data, for a station at the same distance ri from the shower axis. The photon energy Eγ is determined for each measured air-shower event taking into account S(1000) and the reconstructed zenith angle. For this purpose, a look-up table has been constructed using a large simulation sample. Only non-preshowering photon events are used, which are weighted according to a reference spectrum ∝ E −2 . In Fig. 5, the distributions of the two observables are shown as a function of the photon energy. In particular ∆ shows a good separation between photons and data. Finally, the two variables are combined using a Fisher discriminant analysis with the burnt sample representing the background and photon simulations the signal.
The analysis is applied to SD data collected between 1 January 2004 and 30 June 2020. Only air-shower events with a reconstructed zenith angle between 30 • and 60 • are taken into account to ensure that the majority of possible selected photon-induced showers reach their maximum development before reaching ground level. A number of selection criteria are applied to ensure a reliable reconstruction of the two 19    and L LDF (right), i.e., the observables used in the search for photons based on SD-only data, as a function of the photon energy E γ for simulated primary photons (in blue) and a fraction of the data sample (in red) that is used as a burnt sample; the bands represent one standard deviation of the photon distributions. Bottom: distributions of the Fisher discriminant for the burnt sample (grey), the search sample (red) and simulated primary photons (non-preshowering in blue and preshowering in light blue), weighted with an E −2 spectrum; the search sample and the photon distributions are scaled as to have the same integral as the burn sample one; the vertical line indicates the value of the photon-candidate cut; the dashed line shows the result of the fit of an exponential to the 5% of events in the burnt sample with the largest values of the Fisher discriminant. For more details, see [25].
observables. These criteria are described in detail in [25]. Overall, the data sample (search sample) consists of 48,061 selected events with a photon energy Eγ ≥ 10 19 eV, excluding the burnt sample which consists of 886 events (1.8 % of the total number of selected events). The results of the analysis are shown in Fig. 5, bottom. 16, 2 and 0 events above energy thresholds of 1, 2 and 4×10 19 eV, respectively, had a value of the Fisher discriminant above the photon-candidate cut, which was fixed to the median of the Fisher distribution for non-preshowering primary photons (shown as the solid black line in Fig. 5, bottom). The number of observed candidate events is in statistical agreement with what is expected from the fit of an exponential to the tail of the distribution of the Fisher discriminant for the burnt sample. In addition, no peak-like features, which would indicate the presence of a photon population, are observed above the fall-off of the distribution. Overall, the results are consistent with the expectation for a background of UHE protons and nuclei, hence upper limits on the integral flux of UHE photons are determined. To calculate these upper limits, the signal efficiency of the analysis is required. The efficiency has been determined from simulations, and it increases from 0.26 for a threshold energy of 10 19 eV to 0.39 for 4×10 19 eV, under the assumption of a power-law spectrum ∝ E −2 . Upper limits on the integral photon flux are placed at threshold energies of 1,

Summary of the searches for a diffuse flux of UHE photons
The upper limits on the integral photon flux derived through the three analyses discussed in the previous sections are compiled in Tab. 1 and shown in Fig. 6, together with upper limits published by other experiments. The Pierre Auger Observatory currently provides the most stringent limits over a wide energy range, spanning from 2×10 17 eV to the highest energies. In addition, the set of upper limits derived from HeCo data closed the gap between the upper limits at ultra-high energies (derived from hybrid and SD data) and those determined by smaller air-shower experiments such as KASCADE-Grande, leading to a full coverage of the aforementioned energy range. It is worth mentioning here that extensive systematic studies have been performed to test the robustness of the analyses and their results against various sources of uncertainties, for example in the hadronic interaction models used in the air-shower simulations or in the reconstruction of the different observables. Overall, the results proved to be very robust, more details on these studies can be found in [23][24][25].
For comparison, also the expected fluxes of UHE photons under different theoretical assumptions are shown in Fig. 6. First, we discuss briefly the expected fluxes resulting from interactions of UHECRs with the background photon fields permeating the Universe, most notably the cosmic microwave background [40,41]. In Fig. 6, the expectations for two different pure-proton scenarios [2,5] are shown, as well as a scenario involving a mixed composition at the sources [3]. While the experimental sensitivities reached above ≈ 3×10 18 eV start to approach or already constrain the most optimistic expectations of the cosmogenic photon flux from protons, they are about 1 to 1.5 orders of magnitude above those from the mixed-composition model. Another cosmogenic flux is that from the interactions of UHECRs with the matter traversed in the Galactic plane [4]. While this flux becomes comparable to the aforementioned ones below 10 18 eV, they are still two to three orders of magnitude below the current upper limits in this energy region. Finally, UHE photons could also result from the decay of super-heavy dark matter (SHDM) particles. It should be noted that previous upper limits on the incoming photon flux already severely constrained non-acceleration models in general, and SHDM models in particular, trying to explain the origin of cosmic rays at the highest energies (see, e.g., [21,22]). With the upper limits on the incoming photon flux, it is possible to constrain the phase space of mass and lifetime of the SHDM particles [7]. As an example, we show the expectations for three different assumptions: For a hadronic decay (X → qq), we show the expected fluxes according to [38] for a mass MX of the SHDM particles of 10 10 GeV and a lifetime τX of 3×10 21 yr, as well as for MX = 10 12 GeV and a lifetime τX = 10 23 yr. Both combinations are currently allowed. Since also a decay into leptons (X → νν) is possible, we show the expected flux according to [39] for MX = 10 10 GeV and a lifetime τX = 3×10 21 yr. As the sensitivity of current photon searches increases, it will be possible to further  [34], EAS-MSU (magenta triangles) [35]) and Telescope Array (green squares from [36] and turquoise squares from [37]). The ranges of expected GZK photon fluxes under the assumption of two different pure-proton scenario are shown as the red and gray bands (following [2] and [5], respectively). The green band shows the expected GZK photon flux assuming a mixed composition that would fit the Auger data [3], while the blue band denotes the range of photon fluxes that would be expected from cosmic-ray interactions with matter in the Milky Way [4]. In addition, the expected photon fluxes from the decay of super-heavy dark matter particles are included (decay into hadrons, X → qq, based on [38] constrain these values [7].

Searches for UHE photons from point sources and transient events
The analyses discussed in the previous section are searches for a diffuse, i.e., direction-independent flux of UHE photons. Naturally, the arrival direction of a cosmic particle carries important information. Photons, like neutrinos, are neutral particles. They are therefore not deflected by galactic and extragalactic magnetic fields and point right back at their production site. This can be used to search for point sources of UHE photons by simply looking for an excess of events from a certain direction [42], taking into account the angular resolution of the experiment. Complementary to a blind search over the full visible sky, the search for photons can also be restricted to the directions of putative sources to reduce the statistical penalty [43]. In the following sections, we briefly summarize these two analyses. In addition, we discuss follow-up searches for UHE photons from transient events, for example in association with gravitational-wave events [44].

A search for point sources of EeV photons
The blind search for point sources of UHE photons is performed in the energy range between 10 17.3 eV and 10 18.5 eV using hybrid data collected between January 2005 and September 2011. The data set covers a Figure 7: Celestial map, in Galactic coordinates, of upper limits on the incoming photon flux [42]. The white regions indicate regions of the sky that are either not in the field of view of the Pierre Auger Observatory (northern hemisphere) or omitted in this analysis (southern celestial pole). For more details, see [42].
declination range between −85 • and +20 • . The average angular resolution of this data set is 0.7 • . To reduce the contamination of (isotropically distributed) hadronic background events, photon-like air showers are selected using a BDT. The main input variables of the BDT are Xmax and S b , complemented by additional observables based on the fit of a Greisen function to the recorded longitudinal profile, and the ratio of the early-arriving to the late-arriving signal in the SD station with the highest signal [42]. The selection of photon-like events is optimized for each direction by taking into account the expected number of background events, which has been derived using the scrambling technique [45]. No evidence for an excess of photon-like events has been found for any direction within the declination band specified before. The resulting upper limits on the flux of UHE photons for each direction are shown in Fig. 7. The average upper limit on the particle flux is 0.035 km −2 yr −1 with a maximum of 0.14 km −2 yr −1 , corresponding to upper limits on the energy flux of 0.06 eV cm −2 s −1 (average) and 0.25 eV cm −2 s −1 (maximum), under the assumption of an energy spectrum following a power law with spectral index −2.

A targeted search for point sources of EeV photons with the Pierre Auger Observatory
The targeted search for point sources of UHE photons follows the same analysis logic as the blind search summarized before, albeit with a larger data set (January 2005 to December 2013). To reduce the statistical penalty of looking at all directions in the visible sky, the targeted search is restricted to 12 pre-defined target classes, containing in total 364 targets. Since the attenuation length of photons in the energy range considered here (10 17.3 to 10 18.5 eV, same as before) varies between 90 and 900 kpc [43], these target classes contain mostly galactic sources such as, e.g., millisecond pulsars, γ-ray pulsars, and low-mass and high-mass X-ray binaries as well as the Galactic center. In addition, two nearby extragalactic target sets are included: three powerful γ-ray emitters in the Large Magellanic Cloud and the core region of Centaurus A. The different target classes are listed in Tab. 2.
A p-value pi is assigned to each candidate source i of a target set, taking into account the observed number of events from this target direction as well as the expected number of background events. The pvalues of all targets in a set are combined with and without statistical weights, which take into account both the measured electromagnetic flux from the source (taken from astrophysical catalogs) and the directional exposure for photons, derived from simulations. The combined weighted probability Pw is the fraction of isotropic simulations yielding a weighted product that is not greater than the measured weighted product. The combined unweighted probability P is calculated similarly, but with equal weights for all targets. The results of the analysis for each of the 12 target sets are shown in Tab. 2, along with information about the target with the smallest p-value in each set. In addition, the penalized p-values p * are given, i.e., the chance probability that one or more of the targets in the set have a p-value less than p under the assumption of a uniform probability distribution. No combined p-value (weighted and unweighted) nor any individual  Table 2: Combined unweighted probabilities P and weighted probabilities P w for the 12 target sets analyzed in [43]. In addition, selected information on the most significant target from each target set is given: the unpenalized (p) and penalized (p * ) p-values and the derived upper limit on the photon flux at 95 % C.L.. More details on the most significant targets, e.g. the galactic coordinates and upper limits on the energy flux, can be found in [43].
p-value for a target has a statistical significance as great as 3σ. No target class therefore reveals compelling evidence for photon-emitting sources in the 10 18 eV range. There is also no evidence for one outstanding target in any target set.

Follow-up search for UHE photons from gravitational wave sources with the Pierre Auger Observatory
Since the first direct detection of gravitational waves in 2015 [46], the field of multimessenger astronomy has made tremendous progress. In the past years, the transient sources of gravitational wave events-compact binary mergers of black holes and/or neutron stars-have been analyzed by various astronomical instruments. With its unique exposure to UHE particles, the Pierre Auger Observatory has joined the global multimessenger campaign by searching for UHE neutral particles, in particular with follow-up searches for neutrinos (see, e.g., [47]) and photons [44] in association with gravitational wave events. In this context, the search for photons poses several challenges. Not only is the possible flux of UHE photons from any distant source expected to be heavily attenuated due to interactions with the cosmic background radiation fields, but also the non-negligible background of air-shower events with hadronic origin makes the unambiguous identification of primary photons challenging. The identification of photon candidate events is based on the standard search for photons at the Pierre Auger Observatory using the SD (see Sec. 4.3). In order to still maintain a high sensitivity towards a possible photon signal from a transient source despite the considerable background, a dedicated gravitational wave selection strategy has been developed which accepts only close or well-localized sources. Three classes of accepted gravitational wave events are defined in Fig. 8, left, in the space of the 50 % sky localization region and the luminosity distance. While close sources are the most promising candidates to yield a detectable flux of UHE photons, the detection of a photon in coincidence with a distant but well localized source would give a strong hint towards new physics. Out of all gravitational wave events published in the GWTC-1 and -2 catalogs [48,49], four events-including GW170817, which originated from the merger of two neutron stars [50] and after which a short gamma-ray burst was observed from the same direction in the sky [51]-were selected and analyzed for coincident UHE photons within the time period of one sidereal day after the gravitational wave event. No photon candidate events could be identified. Preliminary upper limits on the spectral fluence within the respective time windows are shown in Fig. 8, right.
A similar analysis ansatz can be used to search for UHE photons in association with other transient events, like flares of the blazar TXS 0506+056 [52]. Neither during the first period of enhanced neutrino activity observed by IceCube from October 2014 to February 2015 nor during the second from March 2017 to September 2017 [53] could coincident photon events be identified in the dataset collected with the SD. However, at Gpc scales, no flux of UHE photons could possibly be detected at Earth without new physics processes altering the attenuation of photons in the extragalactic medium. The three classes of selected gravitational wave sources in the follow-up search for photons in association with gravitational wave events, as defined by their 50 % localization region (Ω 50% ) and luminosity distance (D L ); the circled markers in the acceptance region mark the events which had at least some overlap with the field of view of the SD at any time. Right: Preliminary upper limits on the spectral fluence of UHE photons at Earth for each of the selected gravitational wave sources; the uncertainty bars include both the directional uncertainty of the gravitational wave event (blue) and the uncertainty due to the choice of the spectral index used to calculate the spectral fluence (red); for the second event, the uncertainty bars extend beyond the plotted range, since this source is located right at the edge of the field of view of the Pierre Auger Observatory. For more details, see [44].

Outlook
Future improvements to the searches for ultra-high-energy photons summarized in this review can naturally be expected from using larger datasets, profiting from the constant increase in exposure over time with which the upper limits scale inversely. Assuming an increase by a factor of two in exposure, this would translate directly into upper limits that are lower by the same factor of two, in the absence of photon candidate events. In addition, the ongoing detector upgrade of the Pierre Auger Observatory, dubbed AugerPrime [54,55] will play a major role in the future. A key part of this upgrade is the installation of scintillation detectors on top of the water-Cherenkov detector stations of the SD. The current photon searches discussed before already exploit the well-known differences in these components between photon-and hadron-induced air showers (see Sec. 2), albeit in a rather indirect way. AugerPrime will allow for a more direct access, which will lead to an overall better separation between photon-induced air showers and the vast hadronic background. In addition, the detector stations will be equipped with radio antennas to measure of the radio signals emitted by an air shower, which act as a proxy to the electromagnetic component and can therefore also be exploited in searches for UHE photons. All of these efforts combined will significantly improve the upper limits on the incoming photon flux or, in the best case, lead to the first unambiguous detection of photons at ultra-high energies.