Multi-messenger studies with the Pierre Auger Observatory

. The combination of data from observatories measuring ultra-high energy cosmic rays, photons, neutrinos and gravitational waves has provided new insights into the most extreme phenomena in the Universe. Sharing information within a broad community is the foundation of the multi-messenger approach. The Pierre Auger Observatory, the world’s largest cosmic ray detector, provides sensitivity to photons and neutrinos above 10 17 eV, thus contributing signiﬁcantly to this joint e ↵ ort. The latest results from di ↵ use and targeted searches will be reviewed here, along with results from follow-up analyses and future perspectives. In particular, pre-liminary limits on photon ﬂuence from a selection of gravitational wave sources detected by LIGO / Virgo and results of the search for ultra-high energy neutrinos from binary black hole mergers will be presented.


Introduction
In the era of multi-messenger astronomy, ultra-high energy cosmic rays (UHECRs) o↵er the unique opportunity to investigate the nature of astrophysical sources, and particle interactions, in an energy range far beyond that covered by current particle accelerators. The Pierre Auger Observatory [1], the world's largest cosmic ray detector, combines in a hybrid design the information from fluorescence telescopes, observing the longitudinal profile of extensive air showers, with a surface array, measuring the lateral distributions of secondary particles at the ground.
The Observatory consists of a surface detector (SD) array of about 3000 km 2 surrounded by 27 air fluorescence detector (FD) telescopes grouped in four sites, providing a powerful instrument for extensive air shower detection. The SD comprises ⇠ 1600 water-Cherenkov detectors separated by 1500 m (SD-1500) in a triangular grid, plus a smaller nested array of 61 additional detectors spaced by 750 m (SD-750) covering an area of 23.5 km 2 . The FD provides a nearly calorimetric estimate of the primary energy, almost independent of the assumptions on hadronic interaction models at the highest energies. The operation of the FD is however restricted to a duty cycle of about 15% because it can only operate during clear, moonless nights and it is limited by the atmospheric conditions while the SD measurements are made for about 100% of the time. The hybrid paradigm relies on the fact that it is possible to calibrate the SD signal using the events simultaneously recorded by the FD, the so-called hybrid events, avoiding to a large extent the use of Monte Carlo simulations and allowing us to achieve good control of the systematic uncertainties in the energy scale [2]. ⇤ e-mail: lorenzo.perrone@le.infn.it ⇤⇤ e-mail: spokespersons@auger.org

Multi-messenger astronomy
UHECRs are mostly protons and other atomic nuclei; a fraction of photons would be an ideal probe for fundamental questions in particle physics and cosmology, such as the nature of dark matter and the possibility of Lorentz invariance violation (LIV). Indeed, the decay of superheavy DM particles could result in large fluxes of UHE Standard Model particles [3], with photons and neutrinos dominating the final state. Also, UHE photons could result from the decay of hadrons produced when protons with ⇠ 10 times the final photon energy interact with microwave/IR/visible/UV background photons, either within the source environment or during their propagation through the intergalactic medium. In the Standard Model, such UHE photons are expected to be quickly absorbed by background photons, but LIV phenomena could inhibit such interactions, resulting in a significantly higher fraction of photons in the UHECR flux.
Only UHECR observatories have capabilities to detect photons at these energies, through the extensive air showers they produce in the atmosphere. In this context, the Pierre Auger Observatory has the highest sensitivity to photon primaries. Due to the lower multiplicity of electromagnetic interactions (compared to hadronic), photons induced showers develop deeper in the atmosphere and have fewer muons thus yielding a larger depth at shower maximum, X max , and a steeper lateral distribution function. These distinctive features are also reflected in a longer rise time for the signal shape observed with SD. Since the start of data collection, various searches for ultra-high-energy with energy above 10 17 eV photons have been performed, either for a di↵use flux, for point sources or for photons associated with transient events. A very recent review of current results can be found here [4] and details about the detection techniques for specific photon and neutrino searches have been presented at this conference [7]. As an example, the current status of the search for a di↵use flux of photons is shown in Fig. 1, top panel. Upper limits set by the fluorescence detector operating in hybrid mode play an important role at the lowest and intermediate energies, up the EeV range, those set by the surface detector dominate at the highest energies, above 10 19 eV.
Several model predictions in the context of dark matter scenarions are severely constrained by the Auger results and also the expectation of photons originating from the interaction of UHECRs with the interstallar medium start to be in tension with a scenario of a proton-dominated chemical composition. Thus, the search for photons can provide strong bounds on fundamental physics and, at the same time, valuable indications on UHECRs mass composition.
While photons have a limited horizon due to their interaction with the interstellar radiation field (few Mpc at 10 18 eV), neutrinos can travel substantially unabsorbed over cosmological distances. Their fluxes depend on the properties of the sources and on the composition of the primary beam. Neutrinos can be searched for by selecting "horizontal showers" with a large electromagnetic com-ponent close to the detector. This provides a clear identifier of a neutrino-induced shower as, for hadronic cascades, the electromagnetic component would be almost completely absorbed and only muons could be detected in the SD. Two main categories of events are considered: Earth-skimming, induced by ⌫ ⌧ travelling upwards from below the Earth crust in directions just below the horizon and producing a ⌧ lepton which can then generate a shower above the SD, and downward-going events due to neutrinos of any flavour. Fig. 1, bottom panel, shows the most recent upper limits on cosmogenic neutrinos [6], along with model predictions and limits from other experiments.
Given the very good angular resolution of the Observatory (. 1 ) specific follow-up searches were performed with the SD detector looking for point-like sources in spacial and time coincidence with gravitational waves (GW) events. These topics will be reviewed in the following sections.

Search for photons from GW events
A subset of the most interesting GW events from the LIGO/Virgo collaborations was selected based on their localization quality and distance [9]. Time periods of 1000 s around and 1 day after the GW events are analyzed. To maximize the photon-hadron separation power, the observables are combined using a Fisher discriminant analysis taking advantage from the deviation of the expected photon lateral distribution function and of the observed signal risetime from a benchmark reference value measured in data. The distributions of the Fisher discriminant of data events and a set of simulated photon events are shown in Fig. 2. A typical photon-induced air shower is expected to have a significantly larger Fisher discriminant value than the average event found in data. Then, a thresh- old value may be placed to define which events will be accepted as "photon candidate events". The photon dis- crimination method is optimized for air showers with incident zenith angles between 30 and 60 and photon energies above 10 19 eV. Following the approach described in [8], out of all air shower events recorded during a period of 16.5 yr (from January 2004 to June 2020), 16 events passed the photon candidate cut which was placed at the median of the distribution of photon simulations (vertical line in Fig. 2.). This number was found to be consistent with the expected hadronic background.
The GW events considered in the analysis described in [9] were recorded by Advanced LIGO and Virgo during their first three observation runs and published in three gravitational wave transient catalogs: GWTC-1 [10], GWTC-2 [11] with its second revision GWTC-2.1 [12], and GWTC-3 [13]. While the first catalog covers the observations of the first two runs O1 (from 2015 September 12 to 2016 January 19) and O2 (from 2016 November 30 to 2017 August 25), the third observation run has been split in two parts, O3a (from 2019 April 1 to 2019 October 1) and O3b (from 2019 November 1 to 2020 March 27), with a maintenance break of 1 month in between. In order to keep the sensitivity to a possible photon signal as high as possible, GW events are additionally selected by their localization quality and distance. Close and well localized sources are preferred over distant and poorly localized ones. Thus, optimal results can be obtained while keeping the expected background at a reasonable level. Four classes of accepted GW events are defined here for which the 50% localization region is analyzed for coincident air shower events (see Fig. 3). These selection criteria can be summarized as: (D L < 1 and ⌦ 50% < 100 deg 2 ) short (class I), (D L < 1 and ⌦ 50% < 20 deg 2 ) long (class II), (D L < 180 Mpc and ⌦ 50% < 100 deg 2 ) long (class III), (D L < 40 Mpc and ⌦ 50% < 720 deg 2 ) long,short (class IV), where "long" ("short") stands for a time window of 1 day (1000 s) starting after (before) 500 s from the time of a GW event. Since classically no photon signal is expected from very distant sources, the fact of considering large distances (as for class I and II) keeps a window open for potential discoveries of new physics, while the maximum distance chosen for class III corresponds to the most optimistic horizon for photon primaries at the highest energies. A final class of accepted GW events (class IV) allows well identified sources within a distance providing a realistic chance of observing a potential UHE photon flux or at least placing strong physical constraints on the fraction of energy transferred into UHE photons. The value of 50 Mpc is defined by the maximum distance a source like GW170817 may have so that the fraction of energy transferred into UHE photons could still be constrained by a non-observation of photons at the SD array. In total, 23 GW events are classified in at least one of the classes described above. Out of these, 7 of their localization contours were at least partly covered by the Auger SD field of view in the long time window and 3 in the short time window. All 7 events in the long time window are found in class II while one of those, GW170817, also passes the selection criteria for classes I, III and IV. The 3 events in the short time window belong exclusively to class I, i.e. no event falls in both time windows. No coincident air showers with energy above 10 19 eV occured for any source in either of the time windows. This is well in agreement with the expected amount of chance coincidences. Consequently, also none of the 16 photon candidate events from [8] was found to be coincident with any of the selected GWs. Following this non-observation of coincident events, for each GW source an upper limit on the number of photons can be placed using the Feldman and Cousins approach [14]. It can then be converted into photon flux upper limits using the directional exposure calculated in [9] assuming for example a power-law shape with spectral index ↵ for the primary photon flux In addition, an upper limit on the spectral fluence F UL of UHE photons arriving from a given source at the Earth can be derived from the flux upper limit by integrating over the chosen time window [t 0 , t 1 ] and energy range [E 0 , E 1 ]: While no assumption on the time dependence of the flux is made, the extrapolation of the flux limits (derived for the time of transit of the source through the field of view) to the full time window is done under the assumption that the average flux during the period for which the source has been in the field of view is representative for the whole time window. The limits on the spectral fluence depend on the exact direction of the GW source and on the assumed spectral shape of the UHE photon flux. Hence, in Fig. 4 the results for F UL are shown for all possible source directions within each localization contour and a variation of the spectral index ↵ between -2.3 and -1.7, for both time windows. In the long time window, all localization regions have been fully covered by the field of view except for GW170818. Hence, this event could not be constrained for all source directions within the ⌦ 50%contour. All three GW sources in the short time window have contours which partly leak out of the field of view. No coincidences are observed. Upper limits on the UHE photon fluence from a GW event are typically at the level of ⇠ 7 MeV cm 2 for the time period 1000 s and of ⇠ 35 MeV cm 2 for the time period of 1 day. Due to the proximity of the binary neutron star merger GW170817, the energy of the source transferred into UHE photons above 40 EeV is constrained to be less than 20% of its total gravitational wave energy. These are the first limits on UHE photons from GW sources.

Search for neutrinos from GW events
The capability of the Pierre Auger Observatory to detect and identify UHE neutrinos has been shown and discussed in section 2 for the di↵use search. The sensitivity to UHE neutrinos covers a large volume from well above 90% of the Sky, in a declination range from -90 to +60 and an accurate knowledge of the Observatory directionaly sensitivity is required for targeted neutrino search. Fig. 5 shows the upper limits derived for point sources as a function of declination [15], under the assumption of a neutrino flux following a power law / E 2 . These limits are derived at Figure 5. Upper limits on the neutrino flux normalization from steady sources [15]. Limits are calculated for a point source at given declination and single flavour neutrinos. energies higher than the ones covered by ANTARES and IceCube (see the corresponding upper limits also shown in Fig. 5), thus demonstrating the complementarity potential of the Pierre Auger Observatory as neutrino detector in the multi-messenger astronomy joint enterprise.
The UHE neutrino detection capability of the Pierre Auger Observatory makes it also a promising instrument for searches of transient sources such as UHE neutrinos from binary black hole (BBH) mergers that have been observed through the emission of gravitational waves. A dedicated search has been performed for two di↵erent durations after each merger: 24-hour and 60-days [16]. By accounting for the time dependent exposure to neutrinos emitted from 3D localization probability distributions provided by LIGO/Virgo, the UHE neutrino luminosity is derived as a function of time after the merger, assuming a universal emission model.
After the first detection of a direct GW signal from a BBH merger, the LIGO and Virgo Collaborations carried out a total of three observational runs (O1, O2, and O3), as reported in section 3. In O2, the first GW signal from a binary neutron star merger, GW170817, was detected, accompanied by an unprecedented multi-wavelength and multi-messenger observation and search campaign [17] featuring a rich variety of observations that remain unique to date. Part of this was also the follow-up search for high and ultra-high energy neutrinos [18]. No neutrino candidates with significant temporal and directional correlation have been observed, leading to constraints in the astrophysical parameter space of the source. In addition, the detection of GW signals from 83 BBH mergers, the most prevalent type of GW events, have been reported and published [11][12][13] by LIGO/Virgo.
The Pierre Auger Observatory sensitivity across the field of view is not uniform (it is much enhanced close to the horizon), and it depends on the transit and residence time through the field of view of a specific detection channel. An example is again the follow-up search of GW170817, which was located in the most sensitive region directly below the horizon at the time of the merger, allowing to set very stringent UHE neutrino limits in particular for short-term emission. Further objects with follow-up searches for UHE neutrinos include the neutrino-emitting blazar TXS 0506+056 [19].
To quantify the direction-dependent sensitivity to UHE neutrinos for the case of the Pierre Auger Observatory, the e↵ective area A e f f is defined as a function of the local zenith angle and neutrino energy, through the relation between the rate of detected neutrino events and the neutrino energy spectral flux ⌫ : The average of A e f f over the year 2016 is shown in Fig. 6 top panel, as a function of neutrino energy for various zenith angles. It is overall largest for ⌧ neutrinos arriving from directly below the horizon, peaking at a zenith angle of 90.8 . At this angle, the extensive air showers induced by UHE neutrinos are very well distinguishable from those induced by UHECRs, while the Earth is more transparent to UHE neutrinos than at higher zenith angles, especially at the highest energies. This is the reason why the e↵ective area generally decreases quickly for zenith angles above 91 . Here, UHE neutrino energy spectra are assumed to be / E 2 which is a common benchmark assumption for neutrino analyses at the highest energies. The e↵ective area can then be folded with the assumed spectrum, yielding to an e↵ective area per energy A A is proportional to the expected rate of detected neutrinos as it combines the e ciency and the e↵ective geometric area of the SD. A is shown in Fig. 6, bottom as a function of the zenith angle, using the benchmark A e f f from the top panel of the same figure. The non-observation of neutrino candidates with a background low enough to be compatible with zero allows to set an upper limit on the expected number of induced neutrinos that can be converted in terms of the UHE neutrino luminosity of the sources. This can be considered as an astrophysical interpretation of the flux limit at Earth that one could also deduce from the non-observation of neutrino candidates at the given search times and directions. The Pierre Auger Observatory approach is based on the assumption of a time-dependent universal UHE neutrino luminosity L up (t) with which each BBH merger emits neutrinos isotropically, starting at the time of the merger and lasting for the duration of the respective search period [16]. The universal isotropic luminosity upper limit obtained by combining all runs exhibits smaller variations than individual sources or single runs, as the variations in the visibilities of all sources are averaged out. The combination of several sources and di↵erent runs substantially improves the sensitivity to the source class of BBH mergers as compared to single sources or runs. The result for the universal isotropic luminosity for the 24-hour (top) and 60-days (bottom) follow-up period is shown in Fig. 7, along with the contributions from the individual LIGO/Virgo observational runs. The variations in L up are caused by the changing visibilities of the sources during a sidereal day. As an example, the most prominent minima in the curve representing the contribution of O2 occur due to the fact that the localizations of the relatively close-by sources GW170608 and GW170814 coincide with the Earth-skimming band at the corresponding times after the merger. The corresponding limit on the total energy emitted in UHE neutrinos per source is ⇠ 5,6 10 51 ergs. This bound does not depend on the search duration.

Search for upward-going extensive air showers with the FD
The FD is sensitive to upward-going events and can be used to search for upward-going showers without limitation on the zenith angle of the emerging induced showers. The first results from this study, triggered in the context of what reported by the ANITA Collaboration in the energy range above 10 17 eV, have been reported in [20]. Given its operation time and wide field of view, the FD has the potential to support or constrain the "anomalous" observations by the ANITA detector [21], interpreted as upward-going air showers of unexplained nature. We have used 14 years of data collected by the FD to search for upward-going showers using a set of quality selection criteria defined using 10% of the full data sample. To distinguish candidates from false positives, calculate the exposure and obtain the expected background, dedicated simulations for signal (upward-going events) and background (downward-going events) have been performed. Results of the analysis after unblinding the data set are presented. Upward-going protons have been simulated with energies between 10 16.5 and 10 18.5 eV, zenith angle between 90 and 180 and height of first interaction between 4 and 9 km. Protons have been chosen because they can be easily adapted to fit any interesting scenarios such as neutri-nos or Beyond Standard Models (BSM) particles. Inclined downward-going events landing behind an FD site can also mimic an upward-going track in the fluorescence telescope so events with zenith angle between 0 and 90 have also been simulated as background. When simulating an upward-going air shower the height of the first interaction point H fi can significantly change the trigger e ciency of the FD. For this reason a double di↵erential exposure has been calculated as where E cal is the energy released by the shower in the air, S gen is the surface area of generation (a square of 100 ⇥ 100 km 2 ), T is the 14 years of operation of the FD, ⌘ is the fraction of events passing the selection and ✓ is the zenith angle. Fig. 8 shows this exposure based on upwardgoing proton simulations as a function of E cal and H f i . As expected the FD exposure increases with energy being null for energies below 10 17 eV. A blind analysis has been per- Figure 8. Double di↵erential exposure with log10(E cal /eV) on the x-axis and the height of first interaction on the y-axis for upward-going events [20].
formed using 10% of FD data to study the background due to misreconstructed downward-going events [20]. Moreover upward-going lasers, used by the Collaboration for atmospheric monitoring, represent another possible source of background. They are mainly shot by two facilities located in the middle of the array and by 4 LIDARs located at each FD site. Lasers are mostly fired vertically and the large majority of them is rejected based on their known time stamp. However it may happen that some lasers leak into the data sample as genuinely upward-going events. An algorithm has been developed to identify and reject the remaining lasers by exploiting the time of each event and its position inside the array. A profile constrained geometry fit has been applied to the selected sample testing if any possible upward-going geometry can explain the event. In that case, downward geometries have been tested too. A variable X = atan( 2log( L down L all )) · 2/⇡ has been defined to compare the likelihood of the two reconstructions so that an event with X = 0 is more likely a downward-going event, while if X = 1 the upward reconstruction is favoured. Fig. 9 shows the distribution of the X variable for the events in the signal and background simulations as well as for the 10% data sample. According to this distribution a cut value of X = 0.55 has been set with an expected background of ⇠ 0.5 events in the full data sample. Figure 9. Distribution of the X variable (see text for details) for the events from 10% of data (black), background simulations (red), signal simulations (green) [20].
The unblinding procedure leads to only one event passing all the quality cuts, including the final one on the X variable. This result is compatible with the expected background and an integral upper limit to the flux of upward-going air showers has been set at 3.6 10 20 cm 2 sr 1 s 1 and 8.5 10 20 cm 2 sr 1 s 1 after weightening the exposure with E 1 cal and E 2 cal , respectively. A specific scenario assuming upward-going air showers are initiated by ⌧ leptons has been also considered and the corresponding sensitivity of the FD detector calculated in detail [22]. Results have been presented also at this conference [23].