Multi-Messenger Physics with the Pierre Auger Observatory

An overview of the multi-messenger capabilities of the Pierre Auger Observatory is presented. The techniques and performance of searching for Ultra-High Energy neutrinos, photons and neutrons are described. Some of the most relevant results are reviewed, such as stringent upper bounds that were placed to a flux of diffuse cosmogenic neutrinos and photons, bounds placed on neutrinos emitted from compact binary mergers that were detected by LIGO and Virgo during their first and second observing runs, as well as searches for high energy photons and neutrons from the Galactic center that constrain the properties of the putative Galactic PeVatron, observed by the H.E.S.S.\ collaboration. The observation of directional correlations between ultra-high energy cosmic rays and either high energy astrophysical neutrinos or specific source populations, weighted by their electromagnetic radiation, are also discussed. They constitute additional multi-messenger approaches aimed at identifying the sources of high energy cosmic rays.


INTRODUCTION
With the discovery of neutrinos from SN1987A [1,2] arriving four hours before the light detected by conventional telescopes it became clear that there was a lot to learn from examining any type of particles and radiation coming from astrophysical objects, and that neutrino detectors could give early alerts that would facilitate the observation of the evolution of such transients from the earliest stages. Besides the different wavelengths of traditional astronomy, neutrinos, cosmic rays, very high energy gamma rays, and gravitational waves provide complementary information to study the most energetic objects of the Universe. While the SN1987A event might be said to mark the onset of multi-messenger astronomy, the term was only introduced in the late 1990s (see, e.g. [3,4]), and the final boost to the field took place quite recently with the emergence of both neutrino [5] and gravitational-wave astronomy [6]. Indeed the discovery of gravitational waves from the merging of a neutron star binary system by LIGO and Virgo [7] triggered a spectacular series of observations in the full electromagnetic spectrum from radio to the very high energy gamma rays [8] and searches for neutrino fluxes with ANTARES, IceCube, and the Pierre Auger Observatory [9]. The combined effort marks an unprecedented leap forward in astrophysics revealing many aspects of the Gamma-Ray Burst (GRB) induced by the merger and its subsequent kilonova. More recently an energetic neutrino candidate was detected in IceCube in coincidence with the powerful blazar TXS-0506+056 during a flare in the very-high-energy gamma-ray band, incompatible with a chance correlation at the 3σ-level [10]. After scanning archival data, possible evidence for enhanced neutrino emission (quantified at 3.5σ-level) was also found from this direction, during an independent 5-month period between September 2014 and March 2015 [11]. It is now quite clear that there is much potential in the combined analysis of data from different experiments and multi-messenger astronomy has excellent prospects of making many more significant contributions in the near future.
The Pierre Auger Observatory is the largest and most precise detector of ultra-high energy air showers, such as those regularly induced by cosmic rays. While cosmic rays are expected to have curved trajectories in their path to Earth and thus lose time correlation with light emission, at the highest energies they could still keep directional information on the as-yet unidentified sources where they are produced. The study of the arrival directions of these particles has revealed deviations from isotropy [12,13] that, combined with other particle searches, could be of relevance for multi-messenger astronomy. Indeed spatial correlations have been searched for between the arrival directions of the highest energy cosmic rays and classes of objects that have been proposed as sources in the very high-energy regime such as Active Galactic Nuclei (AGN) and Starburst Galaxies (SBG) [14]. Notably, it has also been shown that it is possible to detect neutral particles of sufficiently high energy with the Pierre Auger Observatory. By looking at the depth development of the showers, it is relatively easy to identify neutrinos which interact in the lower layers of the atmosphere. Photons can also be discerned from the background of cosmic rays because the produced showers have a reduced number of muons and they develop deeper than cosmic rays in the atmosphere. Finally, there is no known possibility to separate neutron-induced showers from the charged cosmic rays on the basis of the shower development but, since neutrons are directional, it is in principle possible to identify sources of nearby neutrons by looking at an excess from given directions or to exploit potential time and directional correlations. This procedure could also be applied to any type of neutral particles that induce a shower in the atmosphere such as photons.
In this article, we review the capability of the Observatory to search for signals of such particles and discuss the contributions that have been made.

THE PIERRE AUGER OBSERVATORY
The Pierre Auger Observatory is designed to detect extensive air showers produced by primary cosmic rays above 0.1 EeV. It is located near the city of Malarge, Argentina, at a latitude of 35.2 • S, a longitude of 69.2 • W and at an approximate altitude of 1400 m above sea level, or equivalently an atmospheric depth X ground = 880 g cm −2 . It comprises a surface detector (SD) array of 1660 water-Cherenkov stations deployed over a triangular grid of 1.5 km spacing and a system of 27 telescopes grouped in four sites forming the fluorescence detector (FD). The telescopes are erected at the periphery of the Observatory to observe the atmosphere over the full area of 3000 km 2 covered by the SD array. The SD stations sample the density of the secondary particles of the air shower at the ground and are sensitive to the electromagnetic, muonic and hadronic components. The FD observes the longitudinal development of the air shower by detecting the fluorescence and Cherenkov light emitted during the passage of the secondary particles of the shower in the atmosphere. Unlike the SD, the fluorescence telescopes are operated only during clear and moonless nights, for an average duty cycle of about 14 %. The hybrid feature of the Auger Observatory combining these two well-established techniques has proved to be extremely rewarding in making it the most precise instrument to reconstruct the energy, mass, and direction of Ultra-High Energy Cosmic Rays (UHECRs). Details about the different detector components and their performances can be found in [15,16].
Besides its general hybrid capabilities, a feature that makes the Pierre Auger Observatory a very versatile and powerful multi-messenger observatory is the specific design of the SD stations. They are constructed as cylinders of 3.6 m diameter and 1.2 m height filled with 12 tonnes of purified water, each. Charged particles entering a station induce emission of Cherenkov light which is reflected at the walls by a diffusive Tyvek liner, and collected by three 9-inch photomultiplier tubes (PMT) at the top surface and in optical contact with the water. The PMT signals are sampled by flash analog-to-digital converters (FADC) with a time resolution of 25 ns [16] 1 . This provides good discrimination of electrons and muons entering the detector station from the top, a feature which is important to identify muon-poor air showers induced by photons, not only with the FD but also with the SD. Moreover, the detector stations also provide a large cross section for inclined and up-going air showers, a feature that is of key-importance for the detection of neutrino-induced air showers. Both of these aspects will be discussed in more detail below.
The modular structure of the surface detector array and fluorescence telescopes allowed the data taking to start in 2004 in a partial configuration. In 2008 the detector was completed and since then data have been taken, practically without interruption. Once installed, the local stations work practically without external intervention. There are two types of trigger conditions: a local trigger at the level of an individual station (second order or T2 trigger), and a global trigger (third order or T3 trigger) to register events. The T2 trigger condition is satisfied when either the signal exceeds the equivalent of 3.2 Vertical Equivalent Muons (VEM) in at least one time bin of the signal trace in the FADC -the so-called "threshold trigger"or when it exceeds a much lower threshold (0.2 VEM) in at least 13 bins within a 3 µs time window (i.e., 120 bins) -the so-called Time-over-threshold (ToT) trigger. The ToT condition was designed to trigger on signals broad in time, characteristic of the electromagnetic component dominant over the first 1000 g cm −2 of the extensive air shower, and it is crucial for neutrino identification as explained below. The data acquisition system receives the local T2 triggers and builds a global T3 trigger requiring a relatively compact configuration of at least three local stations compatible in time, each satisfying the ToT trigger, or four triggered stations within a time window slightly larger than light travel time with any type of T2 trigger [15]. With the completed array, the global T3 trigger rate is about two events per minute, one third being showers with energies above 3 × 10 17 eV.

NEUTRINOS
Neutrinos can travel in straight paths from the confines of the Universe, are capable of going through large matter depths without absorption, and are excellent messengers from extragalactic sources, where protons or nuclei are thought to be accelerated to high energies. Neutrinos are also produced from cosmic-ray interactions with the cosmic microwave background. Besides pointing to the position of their sources, they can also provide valuable information concerning hadronic acceleration processes, composition, the local environment at the sources, and their evolution with redshift [17].
The idea of using inclined showers to search for neutrino interactions is old [18], and the large potential of the Pierre Auger Observatory for neutrino detection was already recognized in its design stages [19]. Inclined showers induced by cosmic rays were observed in Haverah Park [20], in the early days of extensive air shower arrays. It was at the beginning of the 2000s when the muon patterns of the showers at ground level were sufficiently understood to allow shower reconstruction [21] (see also [22]). These showers traverse large atmospheric depths, and their electromagnetic component gets almost completely absorbed before the shower particles reach ground level. As a result, for zenith angles exceeding ∼ 60 • , the shower front reaching the SD detector stations consists mainly of muons that typically have energies between 20 and 200 GeV and travel tens of km without decaying. These muons are little affected by interactions, except for continuous energy loss and deflections in the Earth's magnetic field [21]. As their radius of curvature is a few thousands km, they do not deviate much from their initial trajectories, keeping the timing of the shower front sharp.
As opposed to cosmic rays, inclined neutrinos can interact deep in the atmosphere because their interaction length exceeds the matter depth of the atmosphere for any zenith angle, θ. The resulting "young" showers are rich in electrons and photons at ground level. These showers induce signals typically spread over much wider time intervals than cosmic rays with the same θ because of multiple scattering of the electrons and photons in the shower front. The broader digitized signal traces recorded at the particle detectors of the SD array resemble those produced by vertical air showers and it is thus relatively easy to single them out from the sharply arriving fronts produced by inclined protons or nuclei. Tau neutrinos interact in the Earth just below the surface and produce tau leptons that escape into the atmosphere, decaying in flight and producing up-going air showers [23,24]. For neutrino energies exceeding 100 PeV, the combined conversion and exit probability is maximal for nearly horizontal ("Earth-skimming") directions that exit with an elevation angle between 1 • and 2 • above the horizontal direction [25]. These showers develop close and almost parallel to the ground. As the shower particles, particularly electrons and photons, spread laterally they can reach the SD stations, inducing characteristic signals which are also easily identified because of their broader time spread.

Selection and identification
The event selection, the neutrino identification, and the exposure calculations are performed separately for Earth-skimming (ES) and down-going (DG) showers [26]. To select ES showers, all triggered events involving at least three stations are considered. The ellipsoidal signal pattern at the ground is required to have a large eccentricity (characterized by a ratio of major and minor axes greater than 5), and the direction of the major axis indicates the azimuthal angle of the event. The apparent speed of the signal at the ground, measured for each pair of stations and averaged over all of the pairs, must be in the interval [0.967, 1.034]c 2 , and the r.m.s. spread must be less than 0.267c 3 . For the DG selection [27] a minimum of four triggered stations is required to reduce noise from random triggers, and two zenith angle groups are made:"Low" (DGL) and "High" (DGH) respectively for 60 • < θ < 75 • and 75 • < θ < 90 • . The selection of DGL showers is just made requiring the reconstructed zenith angle to be in the 58.5 • -76.5 • range allowing for uncertainties. It is obtained fitting a plane to the positions of the triggered stations and the start time of the signals. For the selection of DGH showers, the reconstructed θ must exceed 75 • , the apparent average speed of the signal at the ground must be smaller than 1.044c 4 , its spread must be smaller than 0.08c, and the eccentricity must exceed a value of 3.
The ES (DG) neutrino identification has been optimized using extensive shower simulations for different energies, arrival directions, decay (interaction) depths, and impact parameters, with standard software such as AIRES [28] and CORSIKA [29] and including the SD response of the individual detector stations using specific software [15,30] based on GEANT4 [31]. A single variable has been chosen for the identification A. Aab et al.  procedure in both searches. Its distribution for neutrino simulations is compared to that of regular showers obtained from a small sample of the data assumed to be exclusively cosmic rays. To select neutrino candidates, a cut is made on this variable at a value that would just yield one cosmic-ray in 50 years according to the extrapolated background distribution. For instance, in the ES case, this variable is the average area over peak (AoP ) value of the signal traces of each event. The area is the integrated charge of the signal trace and the peak is its maximum value, both normalized to one for vertical throughgoing muons used for on-line calibration. The cut obtained is AoP > 1.83, ensuring sufficient time spread in the signal to reject cosmic-ray showers [26]. Events with just three stations are further required to have AoP exceeding 1.4. For DG showers, a Fisher discriminant approach combines several variables based on the AoP of a few selected stations. For the DGH case, the multivariate Fisher method [32] uses nine variables obtained with the AoP values of the four earliest stations and a tenth variable accounting for the asymmetry in the time distribution of the earliest and latest stations of the event. For the DGL case, the Fisher variable combines five or six variables obtained from the AoP of the four or five stations closest to the shower core and their product. For DGL showers, 75% of the triggers are also required to be of special "ToT" kind [26]. For optimization purposes of DG neutrinos, further subdivisions have been made using the number of triggered stations for the DGH set and the reconstructed θ for the DGL set. The condition for a neutrino candidate is slightly different in each subset.

Efficiency and Aperture
The effective aperture of the detector (aperture from now on) is defined so that when it is multiplied by the neutrino spectral flux, it gives the energy distribution of the detection rate. For point sources in the DG case, this aperture is simply calculated integrating the interaction and detection probabilities over the SD area. Three separate apertures are calculated according to interaction type and neutrino flavor: 1) neutral current for all flavors and charged current muon neutrinos, where the shower is produced by the nuclear debris (typically carrying about 20% of the neutrino energy); 2) charged current for electron neutrinos where all the energy is converted to the air shower (80% induced by the electron and the rest by nuclear fragments); 3) charged current for tau neutrinos where the shower can be produced either by the nuclear debris (20% of the energy) or, deeper in the atmosphere, by the tau decay (80% of the energy). The three apertures can be added under the assumption of equal fluxes for each neutrino flavor, as expected because of neutrino flavor oscillations in the characteristic long baselines provided by extragalactic distances to plausible sources [34,35].
On the other hand, only tau neutrinos can be efficiently detected in the ES mode. When muon neutrinos convert to high energy muons that exit the Earth, practically all of them also exit the atmosphere without decaying because of their longer lifetime and when electron neutrinos convert to electrons the shower develops rapidly in the Earth before exiting. In the ES case, the aperture is calculated as an integral over the Earth's area (transverse to the flux direction), weighted by the probability of a tau neutrino interacting and producing a tau which subsequently exits the Earth, decays in the atmosphere and induces a shower that triggers the SD and is effectively selected and identified. This is done in two steps. The differential probability of a tau lepton of given energy exiting the Earth has been calculated as a function of θ using simulations of tau neutrino interactions in rock that include regeneration [36]. The tau-exit probability must be integrated over decay distance weighted by the survival probability, the decay probability per unit distance and the probability of detection with the SD. The detection probability for both DG and ES showers includes trigger, selection and identification and it is also obtained from the simulations.
The exit and detection probabilities in the aperture calculation depend strongly on θ. A given source located at right ascension α, and declination δ in equatorial coordinates, is observed with a zenith angle that oscillates with time according to: Here λ = −35.2 • is the latitude of the detector, t is the sidereal time, and T the duration of a sidereal day. Naturally, the aperture is strongly dependent on source position and time. At any given instant the field of view of the Observatory for each of the ES, DGH and DGL channels is limited to the bands corresponding to the zenith angle range of the channel, as displayed in Figs. 2 (left) and 3 for specific examples discussed below.
Typically, the search for neutrinos from point sources is performed on a pre-selected time interval chosen to match plausible emission times according to theoretical models describing the mechanisms acting within the sources. The expected event rate at the detector is calculated by integrating the aperture over the neutrino flux and, if no candidate events are found in a given time interval, a 90% C.L. upper limit to the flux normalization is obtained matching the expected number of detected events over the time interval to 2.39 [37]. An assumption about the energy spectrum of the flux has to be made. It is customary to again assume a canonical E −2 ν . In this case, 90% of the detected events at the Observatory are between 100 PeV and 25 EeV [26]. The most effective channel for neutrino detection is ES because the target mass provided by the Earth is large compared to the atmosphere and also because the threshold energy is low. For the ES channel algorithms have been found to identify events with just three triggered stations. It has been shown that Earth-skimming events can only be observed for zenith angles between 90 • and 95 • , corresponding to declination angles −55 • < δ < 60 • . The exposure is maximal for δ −53 • and δ 55 • . The apertures provided by the DGH and DGL are progressively smaller because, as the zenith angle decreases, the atmosphere offers less target mass in which the neutrinos can interact and be identified. As the search is effectively performed only for θ between 60 • and 95 • , the field of view of the Observatory for neutrinos is limited to equatorial declinations between ∼ −85 • and ∼ 60 • .
When the pre-selected time interval is short, much less than a day, the aperture can be very different depending on the source position and the observation time. The instantaneous aperture reaches its maximum  The selection and identification procedures have been applied to the data registered by the Observatory. As no neutrino candidates have been found, upper limits to the ultra-high energy (UHE) neutrino fluxes have been obtained. Upper limits to the diffuse flux with important implications for some models of UHECR production were first obtained in the ES channel [16,38] and subsequently extended to include the DGH [39] and the DGL [26] channels. Upper limits to point source fluxes were obtained for the ES and DGH channels [40]. Upper limits on neutrino fluxes from point sources are shown in Fig.1 (right) as a function of declination for the search period up to 31 March 2017. Articles with further updates and details on both diffuse and point source searches are to appear soon.

Correlated searches of neutrinos
The detection of gravitational waves from the merging of binary systems has marked the birth of gravitational wave astronomy. The first signals from the merging of two Black Hole (BH) systems GW150914 [6] and GW151226 [42] were reported in 2016. The luminosity distances of the two events were deduced to be 410 +160 −180 and 440 +180 −190 Mpc and the radiated energies 3.0 +0.5 −0.5 and 1.0 +0.1 −0.2 solar masses respectively. They triggered the search for the emission of electromagnetic radiation and neutrinos [43,44], even though the black hole systems are not expected to emit any other type of radiation unless matter debris  Figure 3. Field of view of the Observatory in Earth-skimming and down-going channels at the instant of the detection of neutron star merger GW170817. The small red contour marks the event localization obtained by the Ligo-Virgo collaborations [41] and the black cross is the position of NGC 4993, later correlated to the event by optical telescopes [52].
and magnetic fields can be found in their neighbourhood, possibly remaining from the BH progenitors [45,46].
Neutrinos were searched for with the Pierre Auger Observatory in two periods of time: ±500 s around the UTC times at which the mergers were observed and one day following their occurrence [44]. These intervals were motivated by the association of the merging of compact objects to Gamma-Ray Bursts (GRBs) [47][48][49]. Neutrinos have been postulated to be produced by accelerated cosmic rays interacting with the GRB photons in the prompt phase and with the lower energy photons of the afterglow. The 1000 s time window is an upper bound of the duration of the prompt GRB phase [50,51] while the 1-day window is a conservative bound of the duration of the afterglow [49].
No neutrino candidates were observed in coincidence with these events in either of the time windows. Assuming a flux k GW E −2 ν , where k GW is a normalization constant, 90% C.L. limits on the neutrino emission in the EeV from these events were reported as fluence limits. The spectral fluence is defined as Here φ is the spectral flux and ∆t is the chosen time interval over which the emission is assumed to take place and during which it is also assumed to be constant. In these conditions the spectral fluence is related to the total energy emitted in neutrinos. For the assumed spectrum it is trivial to convert the spectral fluence to an upper limit for the total energy emitted in neutrinos in a given energy band, which can be compared to the energy radiated in gravitational waves. The reported limits were obtained for a one-day period. Due to the poor localization of these events the upper limits, shown in Fig.2 (left) for GW150914, are reported as a function of declination. A similar limit was obtained for GW151226 since the distances to both mergers are quite similar. The most restrictive limit assuming emission during a single day for GW150914 (GW151226) is obtained for δ = −53 • (δ = 55 • ) and would correspond to a total energy radiated in UHE neutrinos smaller than 7.7 × 10 53 (9.7 × 10 53 ) erg when integrated from 100 PeV to 25 EeV. This corresponds to 14.3% (44.1%) of the radiated energy in gravitational waves. The limits found by IceCube and ANTARES apply to lower neutrino energies and give more restrictive limits on the total energy radiated in these neutrinos [43].  [51] in the case of exact alignment of the line of sight to the rotation axis and for selected off-axis viewing angles. The bounds obtained for a 14-day period on the bottom panel are compared to models of longer lived emission [53]. All models have been scaled to 40 Mpc, the distance to the host galaxy NGC 4993. (Modified from [9]).
Models have been proposed where black hole mergers accelerate UHECR [46], and the energy budget of these collapsing events could account for all of the UHECR with a modest fraction of the gravitational wave energy (of order 3%) going to UHECR acceleration. Assuming that a similar fraction goes into UHE neutrinos, about 0.5 events could have been expected at the Observatory in coincidence with GW150914 in the most optimistic scenario. The upper limits were obtained averaging the instantaneous aperture over a day. If the emission time was shorter than one day, more stringent limits would be obtained provided the source was in the neutrino field of view during the neutrino emission, especially if it was in the ES field of view. Due to the poor localization of both events, it is not possible to know if the events were in the neutrino field of view at the detection time. Given the 90% C.L. contour of GW151226 (GW150917) the overlap with the field of view of the neutrino search is 68.9% (13%) for a time window of 1000 s. The overlap over a whole sidereal day is 100% in both cases. The detection of black hole mergers closer to us and with more directional precision could provide more stringent constraints for such models.
The observation of a neutron star merger GW170817 by the LIGO-Virgo collaboration [54] less than two seconds before the short GRB 170817A detected by Fermi-LAT [55,56] from a coincident direction, has had outstanding implications. The succesful alert systems triggered observations and detection in practically all bands of the electromagnetic spectrum in an unprecedented way, marking the onset of a new era of Multimessenger Astronomy [8]. Early optical observations allowed a precise localization slightly off-centered in the galaxy NGC 4993 at equatorial coordinates α(J2000.0) = 13 h 09 m 48 s , δ(J2000.0) = −23 • 22 53. 343 [52]. Searches for neutrino events were also performed at neutrino telescopes and at the Pierre Auger Observatory [9]. At the time of the GW detection, the source was located at a zenith angle of 91.9 • at the Observatory site, just below the horizon and extremely close to the sweet-spot for Earth-skimming neutrinos, see Fig. 3. When considered in a time interval of ±500 s about the detection (93.3 • < θ < 90.4 • ), the EeV exposure is larger than that of dedicated neutrino telescopes and provides the most stringent upper limit to the neutrino fluence at 90% C.L. in the 100 PeV to 25 EeV interval [9], complementary to IceCube and ANTARES, as illustrated in Fig. 4. The analysis was also made for a longer time window of 14 days matching predictions from longer-lived emission processes [57,58] and also displayed in Fig. 4. Neutrino bounds obtained with the Pierre Auger Observatory and other neutrino telescopes are within the range predicted in some models of GRB from neutron star mergers. IceCube bounds allowed constraints to be placed on the orientation angle relative to the jet axis (assumed to be coincident to the rotation axis) for one of the most optimistic models.
Acceleration of hadrons in astrophysical objects inevitably leads to photon and neutrino production with similar energy fluences. On September 17th of 2017, a 290 TeV neutrino was detected at the IceCube telescope [59] in Antarctica and subsequent correlated gamma-ray observations pointed for the first time to a powerful blazar in a flaring state, TXS 0506+056, as a candidate source of the neutrino [10]. This result was further supported by a neutrino burst from the same direction of the sky obtained in a time window of about 158 days when searching for time correlations from the same direction [11]. A search for possible coincident signals at the Auger Observatory was also performed when the results were made public. This source is 21% of the sidereal day in the field of view of the neutrino search but it was not in the field of view at the exact time of the neutrino detection. No neutrinos were found [60]. The details of this search and the bounds obtained are to be reported in a separate article.

Photon identification
Ultra-high energy (UHE) photons are among the possible particles contributing to the flux of cosmic rays and may arise from a number of processes. Firstly, UHE photons are expected from the GZK-process just like UHE neutrinos and can be used to probe the parameter space of UHECR sources. Secondly, they may be produced by hadronic interactions of cosmic rays within the sources or within their local environment. In these cases, the photons are produced on average with around 10 % of the energy of the primary incident proton. Thirdly, a large flux of UHE photons has been predicted in top-down models with UHECR originating from the decay of supermassive particles such as topological defects or Super Heavy Dark Matter (SHDM) particles (see e.g. [61,62]).
As opposed to neutrinos, photons undergo interactions with the extragalactic background light (EBL) inducing electromagnetic cascades. This process makes photons sensitive to the extragalactic environment (e.g. EBL and magnetic fields) but it also limits the volume from which EeV photons may be detected. It is small compared to the GZK-sphere, 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 [63]. The secondary electrons and positrons created by pair production in the photon fields can again interact A. Aab et al. Nine year SD data sample [73]. Previous data from Auger as well as data from TA, AGASA, Yakutsk, and Haverah Park are included for comparison. The lines and shaded regions give the predictions for top-down models and GZK photon fluxes, respectively, assuming different parameters (references can be found in [72]).

Multi-Messenger Physics with the Pierre Auger Observatory
with background photons via inverse-Compton scattering, resulting in an electromagnetic cascade that ends at GeVTeV energies. Comparing these expected diffuse GeV-TeV fluxes measured with instruments sensitive in this energy range provides another important cross-link to γ-astronomy (see e.g. [64]).
The search for UHE photon primaries is based on the different development and particle content of electromagnetic and hadronic air showers. The induced electromagnetic cascades develop slower than hadronic ones so that the shower maximum X max is reached closer to the ground. Proton and photon simulated showers have average X max values that differ by about 200 g cm −2 in the EeV energy range. This difference is enhanced at energies above 10 EeV because of the Landau-Pomeranchuk-Migdal (LPM) effect [65,66]. Above 50 EeV, photons also have a non-negligible probability to convert in the geomagnetic field [67][68][69] producing a bunch of low-energy electromagnetic particles, called pre-shower, entering the atmosphere. The X max of the pre-showered cascades is smaller than for non-converted ones and the separation between the average X max for photons and proton correspondingly reduced. The fluorescence detector (FD) with its high resolution of about 15 g cm −2 and systematic uncertainties of less than about 8 g cm −2 for X max [70,71] is an ideal instrument to discriminate photon-from hadron-induced air showers with high sensitivity even in single events. Thus the FD is able to provide strong constraints on the photon fraction of UHECR.
It is also possible to search for photon showers using the SD because the shower development and nature of the primary cosmic ray affect the content and the shape of the distribution of particles at the ground as a function of the distance from the shower axis (Lateral Distribution Function, LDF). Photon-induced showers generally are expected to have a steeper LDF compared to hadron primaries because of the subdominant role played by the flatter muonic component. The high-energy effects (LPM and pre-showering) affect the average muon content of air showers by typically less than 20-50 % dependent on zenith angle, which is small compared to the factor 5 to 8 difference in muons with regards to hadronic showers [69,74]. However, the different stage of shower development at the ground (relative to X max ) leads to a modification of the observed LDF. Given the steeper LDF and the muon-driven SD triggers, the footprint at the ground, and consequently the number of triggered stations, N stat , is typically smaller for electromagnetic showers [39]. These features can be combined into a single observable S b , that discriminates the photon-and hadron-induced air showers [75]. Discrimination on a shower by shower basis is less efficient compared to the FD, but given the much larger statistics of the SD, strong constraints on the photon fluxes can be provided. For hybrid observations, the discriminating observables X max , S b , and N stat are injected into a Boosted Decision Tree (BDT) algorithm that takes into account also the energy and angular dependencies of the discriminating observables. To identify photons, a cut is defined at the median of the BDT response distribution for simulated photons, as depicted in Fig. 5 (left). We note that the discrepancy between the data and the proton simulations is in agreement with the current experimental indications of a composition varying from light to heavier in the EeV range [76] and with the muon deficit observed in simulations with respect to the Auger data [75,77]. Applying the cut in this conservative way, the signal efficiency remains constant independently of the composition and hadronic model assumptions. Events having a BDT response larger than the median cut are selected as photon candidates. This yields a background contamination of ∼ 0.14 % for proton showers using QGSJET-II-04. This background level overestimates the one expected in data because the composition is generally heavier than that of pure protons and because the interaction models underestimate the muon number in EAS, making the showers look more photon-like. If the multivariate analysis is performed with a mixture of 50 % proton and 50 % iron as input to the training phase, the background contamination reduces to ∼ 0.04 % with the main contribution coming from the smaller values of the simulated X max . At present, such probabilities are to be considered systematic uncertainty to the background of potential photon candidates, i.e. such background events would dilute the photon limits.  [72]. Eight candidates were found in the first two energy intervals. These events are close to the applied photon cut and detailed simulations of hadronic showers will need to be performed to judge the probability that they are caused by UHECR rather than by photons. For now, they are conservatively considered background rather than signal so that the upper limits of the first two energy bins become less stringent. The differential flux limit between 10-30 EeV is found to be 6.80 · 10 −11 GeV cm −2 s −1 sr −1 at 90 % C.L. Some top-down scenarios proposed to explain the origin of trans-GZK cosmic rays (dashed lines) are illustrated and rejected by previous bounds on the photon flux. A recent super-heavy dark matter proposal (SHDM II) developed in the context of an inflationary theory is shown as a long-dashed line [62]. Constraints on the lifetime-and-mass parameter space of SHDM particle will be imposed by current and future limits on the photon flux, as obtained for example in [78].

Photon searches
A prime interest in the search for photons is to identify the first UHE photon point sources or -in case of non-observations -to provide relevant upper limits thereby constraining the source characteristics. As already mentioned, the horizon is limited to only a few Mpc which reaches out to Cen A but not much farther [79]. If the energy spectra of TeV γ sources measured by atmospheric Cherenkov telescopes [80] extend to EeV energies, it is plausible that photons and neutrons from these sources are detectable also by the Auger Observatory. Sources that produce particle fluxes according to an E −2 energy spectrum inject equal energy into each decade. Thus, a measured energy flux of 1 eV cm −2 s −1 in the TeV decade, as is A. Aab et al.  Figure 6. Gamma-ray spectrum from the Galactic center region as measured by the H.E.S.S. collaboration (red points) [81]. The measured photon flux is extrapolated into the EeV range, given the quoted spectral index and its uncertainties (blue shaded region). The Auger limit [79] is indicated by a green line (the green band reflects again the spectral uncertainties of the gamma-ray spectrum). A spectral index with cutoff energy E cut = 2.0 · 10 6 TeV is indicated by the dashed line.

Multi-Messenger Physics with the Pierre Auger Observatory
found for a number of Galactic sources, would result in the same energy flux in the EeV decade if the spectrum continues to such high energies and energy losses during propagation are negligible. A source of specific interest is the galactic center for which the H.E.S.S.-collaboration measured a gamma-ray flux up to about 50 TeV without any observation of a cutoff or a spectral break, suggesting that the Galactic center hosts a peta-electron volt accelerator, called PeVatron [81].
It is still debated whether these photons are produced in hadronic processes. An interesting test for this is provided by a search for neutrons from this direction, because neutrons would necessarily be produced in charge exchange interactions of protons at the source. This will be discussed in section 5. The ratio between photon and neutron emissivities from pp collisions at the same primary proton energy depends primarily on the spectral index of the protons at the source. Assuming a continuation of the parent proton spectrum with a spectral index Γ p 2.4 well beyond the energies observed in the photon spectrum, the photon emissivities are expected to dominate over those of neutrons [82].
To search for photons from a list of predefined target sources [79], more relaxed cuts are applied to the observed events [83]. For each candidate source direction an optimized cut in the multivariate output distribution is determined which depends on the expected number of isotropic background events in that direction. This number is calculated by applying a scrambling technique that takes into account detector efficiencies and aperture features. For each target direction, a top-hat counting region of 1 • is assumed [83]. Averaging over all considered target directions, the multivariate cut is expected to retain 81.4 % of photons while rejecting 95.2 % of background hadrons. In none of the sources and source classes could EeV photons be detected. As an example, the result for the Galactic center is illustrated in Fig. 6. Conservatively, the extrapolation of the H.E.S.S. data to EeV energies does not take into account the increase of the pp cross-section by about 40 % relative to TeV energies. The current upper limit of 0.044 eV cm −2 s −1 (95 % C.L., Γ = 2) at energies above 0.2 EeV can already constrain the allowed parameter space for a flux continuation to EeV energies.
The observation of photon fluxes from individual sources or from stacked sets of targets would have proved that EeV protons are being accelerated at those considered sources within the galaxy or its neighborhood. However, the null results leave open the possibility that protons observed at EeV energies are of extragalactic origin. This is in fact suggested by the large scale anisotropies observed with the Auger Observatory [12]. Moreover, the absence of detectable photon fluxes, as reported here, does not exclude the production of EeV protons within the Galaxy because the derived flux limits are time-averaged values. EeV photons might be produced in transient sources, such as gamma-ray bursts or supernovae, or be aligned in jets not pointing to us. Extending the searches to bursting sources of EHE photons is a goal of ongoing multi-messenger analyses within Auger, briefly addressed in Section 7.

NEUTRONS
Like photons and neutrinos, neutrons travel in straight lines, undeflected by magnetic fields. They produce air showers that are indistinguishable from those produced by protons. Thus, a flux of neutrons from a discrete source would cause an excess of cosmic-ray events around the direction to the source, clustered within the angular resolution of the Observatory. Since free neutrons undergo β-decay with a mean lifetime of about 880 s at rest, the mean travel distance for relativistic neutrons is 9.2 kpc×E n /EeV. The distance from Earth to the Galactic center is about 8.3 kpc, and the radius of the Galaxy is approximately 15 kpc. Thus, sources in part of the Galactic disk, including the Galactic center, should be detectable via neutrons above 1 EeV. Above 2 EeV, the volume for detectable neutron emitters includes most of the Galaxy.
Neutron production is governed by charge exchange interactions of high energy cosmic-ray protons with ambient photons, protons, or nuclei leading to the creation of a π + -meson. The π + takes the positive charge of the proton and a leading neutron emerges with most of the energy that the proton had. The production of neutrons via creation of π + -mesons is necessarily accompanied by photons originating from decay of similarly produced neutral pions. However, photons resulting from the decay of neutral pions acquire only a small fraction of the proton energy, so that the production of neutrons exceeds the hadronic production of photons of the same energy, provided the accelerated proton spectrum falls approximately like 1/E 2 or more steeply with energy. This makes searches of high energy neutrons a highly relevant and sensitive probe for Galactic hadronic cosmic ray accelerators.
Similarly to the targeted search of EeV photon sources, discussed in Section 4, a flux of neutrons from astrophysical sources in the Galaxy can be detected in the Pierre Auger Observatory as an excess of cosmic-ray air showers arriving from the direction of the sources. To avoid the statistical penalty for making many trials, classes of objects were tested in combinations as target sets. Those target sets include msec pulsars (PSR), pulsar wind nebulae (PWN), microquasars, and magnetars (the full list of targets can be found in [84]). In addition, a search for a neutron flux from the Galactic center and from the Galactic plane was performed. Within a target set, each candidate source is weighted in proportion to its electromagnetic flux, its exposure to the Auger Observatory, and its flux attenuation factor due to neutron decay.
None of the searches provided evidence for a neutron flux from any class of candidate sources. Based on the first nine years of data (Jan. 1, 20041, to Oct. 31, 2013, the upper limits on the energy flux from these candidate sources, including the Galactic center, are mostly at a value of 0.10-0.15 eV cm −2 s −1 [84], which is about a factor of 10 below the energy fluxes detected from TeV gamma-ray sources in the Galaxy [80] and at about the level of sensitivity reached with the targeted photon searches, discussed above. If those sources were accelerating protons in the same environment to EeV energies with the 1/E 2 dependence expected for Fermi acceleration, then the energy flux of neutrons in the EeV energy decade would even exceed the energy flux in TeV gamma rays. Similarly, the energy flux of neutrons from the Galactic plane could be derived to less than 0.56 eV cm −2 s −1 (95 % C.L.). This provides a stringent constraint on models for continuous production of EeV protons in the Galaxy because the injection rate of protons into the Galactic disk must be sufficiently strong to balance their escape from the Galaxy. The concomitant neutron emission rate is model dependent. It could exceed the proton emission rate if protons are magnetically bound to the sources and only the produced neutrons escape, yielding EeV protons by their later decays. More likely, however, the neutron luminosity at any fixed energy is less than the proton luminosity. Based on an estimate of the proton emission rate, the limits on the neutron flux from the Galactic plane could be used to derive an upper limit on the ratio η = neutron luminosity/proton luminosity. It was found to be η UL 0.006 [84], which is a significant constraint on models for continuous production of EeV protons in the Galaxy.

DIRECTIONAL CORRELATIONS OF UHECRS WITH NEUTRINOS AND SOURCE CANDIDATES
Even the most energetic cosmic rays can be subject to significant deflection and corresponding energydependent time delays. The strength and properties of the extragalactic magnetic field (EGMF) causing these effects remain largely unknown, and UHECR observed from bursting and continuous sources have already helped to provide constraints to the EGMF (see e.g. [85][86][87]). In section 3.3 we have discussed the time-correlated observation of neutrinos from merging binary systems. In fact, black hole mergers are also expected to produce UHECR and may provide sufficient luminosity to power the UHECR acceleration up to the highest energies [46]. On the other hand, 100 EeV protons from an extragalactic bursting source located at a distance of 30 Mpc are expected to arrive with time delays relative to photons and neutrinos by O(1 − 1000) yrs [85] and may even reach out to 10 6 yrs in extreme cases, so that UHECR could not be included in those combined directional and temporal searches in a straightforward way.

Search for UHECR-Neutrino Correlations
Purely directional correlations between high-energy neutrinos and UHECR still provide an interesting target of opportunity because UHECRs accelerated in astrophysical sources are naturally expected to produce high-energy photons and neutrinos in interactions with the ambient matter and radiation. If neutrinos result from the decays of pions produced in pγ or pp processes, they would carry about 35 % of the proton energy. Hence, neutrinos observed with energies of 30 TeV to 2 PeV, such as observed by IceCube, would have been produced by protons with energies in the 1 -100 PeV range. These energies are much smaller than those registered by the Auger Observatory, and it is possible that the sources that produce the PeV neutrinos do not accelerate CRs up to ultra-high energies. In that case no correlations would be expected. On the other hand, the sources that produce the UHECR can also be expected to produce lower energy CRs. If the observed neutrinos come from the same sources as the UHECR, some degree of correlation in the arrival directions of the highest energy cosmic rays and the observed neutrinos could be expected depending on composition and on the strength of magnetic fields. It should also be noted that neutrinos can reach us from cosmological distances while the UHECR horizon is limited to only a few 100 Mpc. Thus, only a small fraction of the detected neutrinos could be expected to have originated from visible UHECR sources so that correlations, if they were found, cannot be very prominent.
After some initial cross-correlation efforts by the ANTARES collaboration [88,89] different analysis strategies have been devised in an intercollaborative effort of IceCube and Auger [90] that were more recently joined also by the ANTARES and Telescope Array (TA) collaborations [91,92]. In a classical cross-correlation method, the number of observed UHECR-neutrino pairs is counted as a function of the maximal angular separation and then compared to the averaged background, which is calculated once for an isotropic neutrino flux and once for an isotropic UHECR flux. Applying this analysis to the latest combined Auger+TA UHECR and IceCube neutrino data sets, it was found that a maximum departure from the expectation for an isotropic CR flux, keeping the arrival directions of the neutrinos fixed, occurs at an angular distance of 1 • for tracks and 22 • for cascades, with post-trial p-values of 0.48 and 5.4 · 10 −3 , respectively [92]. In parallel, an unbinned likelihood method was used. It is based on a stacking point-source analysis applied to the UHECRs, where the neutrino arrival directions are smeared with a symmetric 2D-Gaussian and stacked as the signal template. The width of the Gaussian is the quadratic sum of the UHECR reconstruction uncertainty and the magnitude of the magnetic deflection, which is assumed to be described by σ MD = 6 • /E CR [100 EeV]. Again, the analysis is applied to track-and cascade-like high-energy neutrino candidates separately. Assuming an isotropic flux of neutrinos, the smallest post-trial p-values are found to be 2.2 · 10 −2 and 1.7 · 10 −2 for the track-and cascade-like events, respectively.
These results are interesting, but the level of correlation is still too small to arrive at a conclusion. It will be interesting to see how the signal evolves as more data are accumulated by the neutrino and UHECR observatories involved. In addition, new improved analysis strategies are being developed making, for example, better use of the measured UHECR composition and energy so that deflections in the EGMF and Galactic magnetic field can be more appropriately accounted for.

Association of UHECR with Source Populations
A related multi-messenger approach aimed at identifying the sources of UHECR from their directional information is to search for correlations with catalogue based astronomical objects. The basic concept is related to the studies presented in the previous section, only that the tracers are specific source populations rather than samples of high energy neutrinos. A set of source candidates that was studied in [14] includes nearby radio-loud AGNs and SBGs 5 , selected as possible sources to accelerate CRs to the highest energies [95,96]. Individual objects found within a given catalogue are expected to contribute differently to the total UHECR flux observed at Earth. Since the CR luminosity of individual sources is not known, we have chosen their observed electromagnetic emission as a proxy. In case of the 17 radio loud γAGN found within 250 Mpc in the 2FHL catalogue from Fermi-LAT [97], the measured integral γ-flux from 50 GeV to 2 TeV was used as a proxy for the UHECR flux. Similarly, 23 (out of a total of 63) SBGs within 250 Mpc and with a flux larger than 0.3 Jy were taken from [98] and weighted with their continuum radio emission at 1.4 GHz. In both cases, attenuation of UHECRs from distant objects in the photon background fields was accounted for using a data-driven scenario that reproduces the composition and spectral constraints obtained by the Observatory [99]. Within this Ansatz, the relative contributions of NGC 4945, NGC 253, and M83 are expected to be about 38, 33, 13 % of the total UHECR flux from SBGs observed at the Auger Observatory. The relative contributions seen at the Observatory from the γAGN would be dominated by the Cen A Core (75 %), followed by M87 (15 %) and Mkn 421 (8 %).
To test these models against the observed UHECR sky maps, a maximum-likelihood analysis was performed with two free parameters aimed at maximizing the degree of correlations of the model maps with UHECR events: a) the fraction of an isotropic component contributing to the total UHECR flux in addition to the source population being tested, and b) the width of a 2D Gaussian smearing around the position of the source candidates. The first free parameter can be interpreted as the contribution from more distant or less luminous sources, and the second parameter effectively accounts for the random scattering of UHECR in the EGMF. This analysis is repeated for a sequence of energy thresholds applied to the UHECR events and the test statistic (TS) as a function of the two parameters is analyzed. The results are shown in Fig. 7. It was found that the SBG-Model rejects the hypothesis of an isotropic sky best at E CR > 39 EeV and yields a post-trial significance of 4.0σ. Here, the SBGs contribute about 10 % to the total flux, and the smearing angle is 12.9 • . In the case of the γAGN model, the best rejection of isotropy is reached at E CR > 60 EeV at a level of 2.7σ with a γAGN fraction to the total flux of about 7 % and a smearing angle 6.9 • . Interestingly, the smearing angle decreases with increasing energy as expected from the larger rigidities of the UHECR in the γAGN model.
The results are very promising and demonstrate departures from a pure isotropic sky with more and more structure becoming visible, starting at above 8 EeV with the observation of a large-scale dipole [12], indications of higher order moments above 16 EeV, [13], blurry spots above 40 EeV, and hints of sharper spots above 60 EeV. There are indications that both SGB and AGN, as well as a mixture of both or similar populations could play a significant role in the production of UHECR. However, it is still too early to conclude about this long standing problem. Future analyses, possibly profiting from full-sky observations with TA, may involve a better modelling of the Galactic magnetic field, account for the change of the composition as a function of energy, and possibly introduce better proxies for the unknown UHECR luminosities, so that 5σ detections are well within reach.

CONCLUSIONS AND PROSPECTS
The Pierre Auger Observatory was designed as a multi-purpose observatory for the high-energy Universe with multi-messenger observations being foreseen from the beginning. In this review, we have summarized a few prominent examples that demonstrate the unprecedented sensitivities to UHE photons, neutrinos, and neutrons. Observations of photons in the so far dark Universe at E γ 10 17 eV or of a neutrino at these energies would be a major breakthrough by itself. However, the non-observation of point sources and diffuse fluxes of photons and neutrinos allowed the derivation of upper bounds that constrain models very effectively. The bounds to diffuse photon and neutrino fluxes have ruled out the top-down models of UHECRs origin motivated by particle physics and also started to constrain the GZK-effect as a dominant process for explaining the observed flux suppression of the most energetic UHECR [33]. Many TeV γ-sources are observed at energy fluxes of the order of 1 eV cm −2 s −1 . Such sources would be visible to the Auger Observatory as strong photon and Galactic neutron sources if their energy spectrum would continue with a Fermi-like energy distribution up to about 10 17 eV. Again, their absence suggests that their maximum source energy does not reach out to the threshold energy of the Observatory and/or that their spectrum is significantly softer than dJ/dE ∝ E −2 .
Point source searches now routinely include also mergers of compact binaries alerted by gravitational wave interferometers. The most spectacular event so far was the neutron star merger GW170817 at a distance of about 40 Mpc. Within the predefined ±500 s search window, the Auger Observatory reached a neutrino flux sensitivity above 100 PeV that was over an order of magnitude higher than of any other neutrino observatory presently operated. Again, the absence of neutrinos at Auger, IceCube and ANTARES allowed constraining the jet properties of the neutron star merger. The third observation run O3 is about to start with many more events being expected in the near future. To accommodate for this, mechanisms are set up to automatically react to GW or other astrophysical alerts and search for neutrinos and photons.
While receiving alerts from a worldwide network of observatories, it is also possible to send alerts from the Observatory. The Auger Observatory is both a triggering and a follow-up partner in the Astrophysical Multimessenger Observatory Network (AMON) [100], which establishes and distributes alerts for immediate follow-up by subscribed observatories with private or Gamma-Ray Coordinates Network 6 notices. AMON registers all vertical (θ ≤ 60 • ), high-quality, Auger SD events, with energy above 3 EeV with a cadence of at most a few minutes. AMON establishes clusters of two or more Auger events received correlated in time and arrival direction as alerts. These are excesses of events expected from neutral (photon or neutron) UHECRs as discussed in Section 4 and Section 5. Unlike the aforementioned neutron and photon searches, this alert channel is sensitive to ultra-high energy neutral transient emission which is less well constrained in the time-integrated searches. For each Auger alert, AMON also receives and transmits the background rate, in other words, the significance of each of the alerts, which depends on the energy and declination of each event. Thus high-significance alerts can be preferentially followed-up. To improve the overall significance of alerts and to enhance future photon and neutrino searches in general, additional improved hardware triggers have been implemented to the existing surface detector array of the Observatory. Besides enhancing the sensitivity to photons and neutrinos, they allow reducing their detection thresholds. Thus, photon alerts for individual ultra-high energy photon candidates from a selection of events which utilizes observables based on the new triggers will be transmitted in near-real time and available for immediate follow-up using the AMON infrastructure. Finally, the infilled array, covering an area of 23.5 km 2 with a 750 m detector spacing, will provide acceptance for neutrinos and photons below the present detection threshold, though with reduced total exposure. It may also be used to generate specific alerts and there are plans to extend the neutrino search to the FD [101].
Presently, the Pierre Auger Observatory is being upgraded to AugerPrime primarily to improve the mass composition measurements and particle physics capabilities with the surface detector array [102][103][104]. For this purpose, 3.8 m 2 plastic scintillation detectors will be installed on each of the existing surface detector stations and the current readout electronics will be replaced by a faster and more powerful one. The new electronics will facilitate three times faster sampling of the signals (120 MHz instead of 40 MHz), will have more monitoring features implemented, and will have more sensitive triggers installed for low-energy showers and those initiated by photons and neutrinos. Moreover, the infilled array will be furnished with large area underground muon detectors [105,106] and each surface detector station will be amended by a radio antenna to provide improved shower information for inclined air showers [107]. AugerPrime will be operated into 2025 and will improve the statistics for composition information and composition enhanced astronomy by about a factor of 10. Clearly, all of these features will open up many new possibilities for improved searches of photons and neutrinos. This suit of enhancements will further strengthen the prominent role of the Pierre Auger Observatory as a multi-messenger observatory for the next decade.