An improved limit to the diffuse flux of ultra-high energy neutrinos from the Pierre Auger Observatory

Neutrinos in the cosmic ray flux with energies near 1 EeV and above are detectable with the Surface Detector array of the Pierre Auger Observatory. We report here on searches through Auger data from 1 January 2004 until 20 June 2013. No neutrino candidates were found, yielding a limit to the diffuse flux of ultra-high energy neutrinos that challenges the Waxman-Bahcall bound predictions. Neutrino identification is attempted using the broad time-structure of the signals expected in the SD stations, and is efficiently done for neutrinos of all flavors interacting in the atmosphere at large zenith angles, as well as for"Earth-skimming"neutrino interactions in the case of tau neutrinos. In this paper the searches for downward-going neutrinos in the zenith angle bins $60^\circ-75^\circ$ and $75^\circ-90^\circ$ as well as for upward-going neutrinos, are combined to give a single limit. The $90\%$ C.L. single-flavor limit to the diffuse flux of ultra-high energy neutrinos with an $E^{-2}$ spectrum in the energy range $1.0 \times 10^{17}$ eV - $2.5 \times 10^{19}$ eV is $E_\nu^2 dN_\nu/dE_\nu<6.4 \times 10^{-9}~ {\rm GeV~ cm^{-2}~ s^{-1}~ sr^{-1}}$.

Neutrinos in the cosmic ray flux with energies near 1 EeV and above are detectable with the Surface Detector array of the Pierre Auger Observatory. We report here on searches through Auger data from 1 January 2004 until 20 June 2013. No neutrino candidates were found, yielding a limit to the diffuse flux of ultra-high energy neutrinos that challenges the Waxman-Bahcall bound predictions. Neutrino identification is attempted using the broad time-structure of the signals expected in the SD stations, and is efficiently done for neutrinos of all flavors interacting in the atmosphere at large zenith angles, as well as for "Earth-skimming" neutrino interactions in the case of tau neutrinos. In this paper the searches for downward-going neutrinos in the zenith angle bins 60 • −75 • and 75 • −90 • as well as for upward-going neutrinos, are combined to give a single limit. The

I. INTRODUCTION
The flux of ultra-high energy cosmic rays (UHECRs) above ∼ 5 × 10 19 eV is known to be suppressed with respect to that extrapolated from lower energies. This feature has been seen in the UHECR spectrum [1,2], with the position of the break being compatible with the Greisen-Zatsepin-Kuzmin (GZK) effect [3], i.e. the interaction of UHECRs with the cosmic microwave background (CMB) radiation. However, other explanations are possible, most prominently a scenario where the limiting energy of the UHECR sources is being observed [4]. Key to distinguishing between these two scenarios is the determination of the composition of the UHECRs [5,6], with the second scenario predicting increasing fractions of primaries heavier than protons as energy increases [4].
Above ∼ 5 × 10 19 eV cosmic-ray protons interact with CMB photons and produce ultra-high energy cosmogenic neutrinos of energies typically 1/20 of the proton energy [7]. Their fluxes are uncertain and at EeV energies they depend mostly on the evolution with redshift z of the unknown sources of UHECRs, and on their spectral features at injection. Protons typically produce more neutrinos than heavier primaries do [8,9], so measurement of the neutrino flux gives information on the nature of the primaries. In this respect the observation of UHE neutrinos can provide further hints on the dominant scenario of UHECR production [9], as well as on the evolution with z of their sources which can help in their identification [9,10].
UHE neutrinos are also expected to be produced in the decay of charged pions created in the interactions of cosmic rays with matter and/or radiation at their potential sources, such as Gamma-Ray Bursts or Active Galactic Nuclei among others [11]. In fact, at tens of EeV, neutrinos may be the only direct probe of the sources of UHECRs at distances farther than ∼ 100 Mpc.
A breakthrough in the field was the recent detection with the IceCube experiment of three neutrinos of energies just above 1 PeV, including a 2 PeV event which is the highest-energy neutrino interaction ever observed, followed by tens of others above ∼ 30 TeV representing a ∼ 5.7 σ excess above atmospheric neutrino background [12]. The measured flux is close to the Waxman-Bahcall * Now at Fermilab, Batavia, IL, USA † Deceased ‡ Now at CERN, Geneva, Switzerland § Also at Vrije Universiteit Brussels, Belgium ¶ auger spokespersons@fnal.gov upper bound to the UHE neutrino flux [13], although with a steeper spectrum [37].
In the EeV energy range, i.e. about three orders of magnitude above the most energetic neutrinos detected in IceCube, neutrinos have so far escaped detection by existing experiments. These can be detected with a variety of techniques [14], among them with arrays of particle detectors at ground.
In this work we report on the search for EeV neutrinos in data taken with the Surface Detector array (SD) of the Pierre Auger Observatory [15]. A blind scan of data from 1 January 2004 up to 20 June 2013 has yielded no neutrino candidates and an updated and stringent limit to the diffuse flux of UHE neutrino flux has been obtained.

II. SEARCHING FOR UHE NEUTRINOS IN AUGER
The concept for identification of neutrinos is rather simple. While protons, heavier nuclei, and even photons interact shortly after entering the atmosphere, neutrinos can initiate showers quite deep in the atmosphere. At large zenith angles the atmosphere is thick enough so that the electromagnetic component of nucleonic cosmic rays gets absorbed and the shower front at ground level is dominated by muons ("old" shower front). On the other hand, showers induced by neutrinos deep in the atmosphere have a considerable amount of electromagnetic component at the ground ("young" shower front). The Surface Detector array (SD) of the Pierre Auger Observatory is not directly sensitive to the muonic and electromagnetic components of the shower separately, nor to the depth at which the shower is initiated. In the ∼1600 water-Cherenkov stations of the SD of the Pierre Auger Observatory, spread over an area of ∼3000 km 2 , separated by 1.5 km and arranged in a triangular grid, the signals produced by the passage of shower particles are digitised with Flash Analog to Digital Converters (FADC) with 25 ns resolution. This allows us to distinguish narrow signals in time induced by inclined showers initiated high in the atmosphere, from the broad signals expected in inclined showers initiated close to the ground.
Applying this simple idea, with the SD of the Pierre Auger Observatory [15] we can efficiently detect inclined showers and search for two types of neutrino-induced showers at energies above about 1 EeV: 1. Earth-skimming (ES) showers induced by tau neutrinos (ν τ ) that travel in a slightly upward direction with respect to ground. ν τ can skim the Earth's crust and interact relatively close to the surface inducing a tau lepton which escapes the Earth and decays in flight in the atmosphere, close to the SD.
2. Showers initiated by any neutrino flavor moving down at large angles with respect to the vertical that interact in the atmosphere close to the surface detector array through charged-current (CC) or neutral-current (NC) interactions. We include here showers induced by ν τ interacting in the mountains surrounding the Pierre Auger Observatory. Although this latter process is exactly equivalent to the "Earth-skimming" mechanism, it is included in this class because such showers are also going downwards. In the following we will refer to all these types of showers as "downward-going" (DG) ν-induced showers.
With the aid of Monte Carlo simulations we have established that this search can be performed efficiently as long as it is restricted to showers with zenith angles θ > 60 • . Due to the characteristics of these showers depending on the zenith angle, the search in this channel was performed in two angular subranges: (a) "low" zenith angle (DGL) corresponding to 60 • < θ < 75 • and (b) "high" zenith angle (DGH) with 75 • < θ < 90 • .
A. General procedure The identification of potential neutrino-induced showers is based on first selecting those events that arrive in rather inclined directions, and then selecting among them those with FADC traces that are spread in time, indicative of the early stage of development of the shower and a clear signature of a deeply interacting neutrino triggering the SD.
First of all, events occurring during periods of data acquisition instabilities [16] are excluded. For the remaining events the FADC traces of the triggered stations are first "cleaned" to remove accidental signals [17] induced mainly by random atmospheric muons arriving closely before or after the shower front. These muons are typically produced in lower energy showers (below the energy threshold of the SD of the Auger Observatory) that arrive by chance in coincidence with the triggering shower. A procedure to select the stations participating in the event described in [17,18] is then applied, with the event accepted if the number of accepted stations N st is at least three (four) in the Earth-skimming (downward-going) selections.
From the pattern (footprint) of stations at ground a length L along the arrival direction of the event and a width W perpendicular to it characterizing the shape of the footprint are extracted [17]. The ratio L/W ∼ 1 in vertical events, increasing gradually as the zenith angle increases. Very inclined events typically have elongated patterns on the ground along the direction of arrival and hence large values of L/W . A cut in L/W is therefore a good discriminator of inclined events. Another indication of inclined events is given by the apparent speed V of the trigger from a station i to a station j, averaged over all pairs (i, j) of stations in the event. This observable denoted as V is obtained from the distance between the stations after projection along L and from the difference in trigger times of the stations. In vertical showers V exceeds the speed of light since all triggers occur at roughly the same time, while in very inclined events V is concentrated around the speed of light. Moreover its Root-Mean-Square (RMS(V )) value is small. For downward-going events only, a cut on the reconstructed zenith angle θ rec is applied [18].
Once inclined showers are selected the next step is to identify young showers. A Time-over-Threshold (ToT) trigger 1 is usually present in SD stations with signals extended in time, while narrow signals induce other local triggers. Also the Area-over-Peak ratio (AoP), defined as the ratio of the integral of the FADC trace to its peak value, normalized to the average signal produced by a single muon, provides an estimate of the spread-in-time of the traces, and serves as an observable to discriminate broad from narrow shower fronts. In particular, a cut on AoP allows the rejection of background signals induced by inclined hadronic showers, in which the muons and their electromagnetic products are concentrated within a short time interval, exhibiting AoP values close to the one measured in signals induced by isolated muons. These observables are used by themselves in the search for ν candidates, or combined in a linear Fisher-discriminant polynomial depending on the selection as described later in this work.
As a general procedure and to optimize the numerical values of the cuts and tune the algorithms needed to separate neutrino-induced showers from the much larger background of hadronic showers, we divided the whole data sample (1 January 2004 -20 June 2013) into two parts (excluding periods of array instability). A selection dependent fraction of the data ∼ 20%, along with Monte Carlo simulations of UHE neutrinos, is dedicated to define the selection algorithm, the most efficient observables and the value of the cuts on them. These data are assumed to be overwhelmingly constituted of background showers. The applied procedure is conservative because the presence of neutrinos in the training data would result in a more severe definition of the selection 1 This trigger is intended to select sequences of small signals in the FADC traces spread in time. It requires at least 13 bins in 120 FADC bins of a sliding window of 3 µs above a threshold of 0.2 I peak V EM (the peak value of the signal expected for a vertical muon crossing the station), in coincidence in 2 out of 3 PMTs [16].
criteria. The remaining fraction of data is not used until the selection procedure is established, and then it is "unblinded" to search for neutrino candidates. We used real data to train the selections instead of Monte Carlo simulations of hadronic showers, the primary reason being that the detector simulation may not account for all possible detector fluctuations that may induce events that constitute a background to UHE neutrinos, while they are contained in data. It is important to remark that this is the same selection procedure and training period as in previous publications [17,18], which is applied in this work to a larger data set.
Regarding the Monte Carlo simulations, the phase space of the neutrino showers reduces to three variables: the neutrino energy E ν , the incidence zenith angle θ and the interaction depth D in the atmosphere for downwardgoing neutrinos, or the altitude h c of the τ decay above ground in the case of Earth-skimming neutrinos. Showers were simulated with energies from log(E τ /eV) = 17 to 20.5 in steps of 0.5, zenith angles from 90.1 • to 95.9 • in steps of 0.01 rad (ES) and from 60 • to 90 • in steps of 0.05 rad (DG). The values of h c range from 0 to 2500 m (in steps of 100 m) whereas D is uniformly distributed along the shower axis in steps of 100 g cm −2 .
We have described the general procedure to search for Earth-skimming ν τ and downward-going ν-induced showers. However the two searches (ES and DG) differ in several aspects that we describe in the following sections.

B. Earth-skimming (ES) neutrinos
With Monte Carlo simulations of UHE ν τ propagating inside the Earth, we have established that τ leptons above the energy threshold of the SD are efficiently produced only at zenith angles between 90 • and 95 • . For this reason, in the Earth-skimming analysis we place very restrictive cuts to select only quasi-horizontal showers with largely elongated footprints: L/W > 5 and V ∈ [0.29, 0.31] m ns −1 with RMS(V )< 0.08 m ns −1 (see Table I In the ES selection, the neutrino identification variables include the fraction of stations with ToT trigger and having AoP> 1.4 for data prior to 31 May 2010 [17]. This fraction is required to be above 60% of the triggered stations in the event. The final choice of the values of these cuts was made by requiring zero background events in the training data sample, corresponding to 1% of the events recorded up to that date. For data beyond  Table I. Grayfilled histogram: the data in the training period. Black histogram: data in the search period. These two distributions are normalised to the same number of events for comparison purposes. Blue histogram: simulated ES ντ events. The dashed vertical line represents the cut on AoP > 1.83 above which a data event is regarded as a neutrino candidate. An exponential fit to the tail of the distribution of training data is also shown as a red dashed line (see text for explanation).
1 June 2010 a new methodology and a new set of efficient selection criteria was established based on an improved and enlarged library of ES simulated ν τ events and on a larger period of training data. In particular, we used the average value of AoP ( AoP ) over all the triggered stations in the event as the main observable to discriminate between hadronic showers and ES neutrinos. The new methodology allows us to place the value of the cut on AoP using the tail of its distribution as obtained in real data (which was seen to be consistent with an exponential shape as shown in Fig. 1). This tail was fitted and extrapolated to find the value of the cut corresponding to less than 1 expected event per 50 yr on the full SD array. As a result, an event is tagged as a neutrino candidate if AoP > 1.83 (see Table I and Fig. 1). The new methodology is not applied to the data prior to 31 May 2010 since that data period was already unblinded to search for UHE neutrinos under the older cuts [17]. Roughly ∼ 95% of the simulated inclined ν τ events producing τ leptons above the energy threshold of the SD are kept after the cut on AoP . The search for neutrinos is clearly not limited by background in this channel.  0.08 plus a further requirement that the reconstructed zenith angle θ rec > 75 • (see [18] and Table I for full details).
In the low zenith angle range (DGL) corresponding to 60 • < θ < 75 • , L/W , V and RMS(V )/ V are less efficient in selecting inclined events than the reconstructed zenith angle θ rec , and for this reason only a cut on θ rec is applied, namely 58.5 • < θ rec < 76.5 • , which includes some allowance to account for the resolution in the angular reconstruction of the simulated neutrino events.
After the inclined shower selection is peformed, the discrimination power is optimized with the aid of the multivariate Fisher discriminant method [19]. A linear combination of observables is constructed which optimizes the separation between background hadronic inclined showers occuring during the downward-going training period, and Monte Carlo simulated ν-induced showers. The method requires as input a set of observables. For that purpose we use variables depending on the dimensionless Area-over-Peak (AoP) observable -as defined above -of the FADC traces.
In the DGH channel, due to the inclination of the shower the electromagnetic component is less attenuated at the locations of the stations that are first hit by a deep inclined shower (early stations) than in the stations that are hit last (late stations). From Monte Carlo simulations of ν−induced showers with θ > 75 • we have established that in the first few early stations the typical AoP values range between 3 and 5, while AoP tends to be closer to 1 in the late stations. Based on this simple observation and as already reported in [18], we have found a good discrimination when the following ten variables are used to construct the linear Fisher discriminant variable F: the AoP and (AoP) 2 of the four stations that trigger first in each event, the product of the four AoPs, and a global parameter that measures the asymmetry between the average AoP of the early stations and those triggering last in the event (see [18] for further details and Table I).
The selection of neutrino candidates in the zenith angle range 60 • < θ < 75 • (DGL) is more challenging since the electromagnetic component of background hadronic showers at ground increases as the zenith angle decreases because the shower crosses less atmosphere before reaching the detector level. Out of all triggered stations of an event in this angular range, the ones closest to the shower core exhibit the highest discrimination power in terms of AoP. In fact it has been observed in Monte Carlo simulations that the first triggered stations can still contain some electromagnetic component for background events and, for this reason, it is not desirable to use them for discrimination purposes. The last ones, even if they are triggered only by muons from a background hadronic shower, can exhibit large values of AoP because they are far from the core where muons are known to arrive with a larger spread in time. Based on the information from Monte Carlo simulations, the variables used in the Fisher discriminant analysis are the individual AoP of the four or five stations (depending on the zenith angle) closest to the core, and their product [20]. In the DGL analysis it is also required that at least 75% of the triggered stations closest to the core have a ToT local trigger [20].
Once the Fisher discriminant F is defined, the next step is to define a numerical value F cut that efficiently separates neutrino candidates from regular hadronic showers. As was done for the variable AoP in the Earth-skimming analysis, F cut was fixed using the tail of the distribution of F in real data, which is consistent with an exponential shape in all cases. An example is shown in Fig. 2. The tail was fitted and extrapolated to find the value of F cut corresponding to less than 1 expected event per 50 yr on the full SD array [18,20]. Roughly ∼ 85% (∼ 60%) of the simulated inclined ν events are kept after the cut on the Fisher variable in the DGH (DGL) selections. The smaller efficiencies for the identification of neutrinos in the DGL selection are due to the more stringent criteria in the angular bin θ ∈ (60 • , 75 • ) needed to reject the larger contamination from cosmic-ray induced showers.

A. Data unblinding
No events survived when the Earth-skimming and downward-going selection criteria explained above and summarized in Table I are applied blindly to the data collected between 1 January 2004 and 20 June 2013. For each selection the corresponding training periods are excluded from the search. After the unblinding we tested the compatibility of the distributions of discriminating observables in the search and training samples. Examples are shown in Fig. 1 for the AoP variable in the Earth-skimming analysis, and in Fig. 2 for the Fisher variable in the DGH analysis. In particular fitting the tails of the corresponding distributions to an exponential, we obtained compatible parameters within 1 σ statistical uncertainties.
B. Exposure calculation

Neutrino identification efficiencies
The selection criteria in Table I, were also applied to neutrino-induced showers simulated with Monte Carlo, and the identification efficiencies ES , DGH , DGL for each channel -defined as the fraction of simulated events passing the cuts -were obtained.
A large set of Monte Carlo simulations of neutrinoinduced showers was performed for this purpose, covering the whole parameter space where the efficiency is expected to be sizeable. In the case of Earth-skimming ν τ induced showers, the efficiency depends on the energy of the emerging τ leptons E τ , on the zenith angle θ and on the altitude of the decay point of the τ above ground. These efficiencies are averaged over azimuthal angle and the τ decay channels. The maximum efficiency that can be reached is 82.6%, the 17.4% remaining corresponds to the channel in which the τ decays into a µ which is un-likely to produce a detectable shower close to ground. In the case of downward-going neutrinos the identification efficiency depends on neutrino flavor, type of interaction (CC or NC), neutrino energy E ν , zenith angle θ, and distance D measured from ground along the shower axis at which the neutrino is forced to interact in the simulations.
The identification efficiencies depend also on time, through the changing configuration of the SD array that was growing steadily since 2004 up to 2008, and because the fraction of working stations -although typically above 95% -is changing continuously with time. Also the continuous monitoring of the array reveals a slight evolution with time of the optical properties of the water-Cherenkov stations (see below). Although the number of working stations and their status are monitored every second and as a consequence the SD configuration is known with very good accuracy at any instant of time, in practice, to avoid having to cope with an impractically large number of configurations, different strategies were devised to calculate in an accurate and less timeconsuming manner the actual identification efficiencies (as explained in [17,18,20]).
The evolution of the optical properties of the water-Cherenkov stations was taken into account in an effective way in the calculation of the exposure. The main effect of this evolution is a decrease with time of the decay-time of the light as obtained from the monitoring data that revealed a continuous decrease of ∼ 10% from 2004 until the end of the data period used in this work (20 June 2013). This induces a reduction of the AoP and, as a consequence, the trigger efficiency changes with time. These changes were accounted for in the calculation of the exposure by dividing the whole data set into three separate periods and assuming that in each of them the decay-time of the light in the tank remained approximately constant as seen in data. A conservative approach was adopted by choosing constant values of the light decay-time below the actual curve in the three periods.

Combination of selections
In previous publications [17,18,20] the fraction of νinduced Monte Carlo events identified as neutrino candidates was obtained by applying each particular set of selection criteria (ES, DGH, DGL) only to its corresponding set of simulated showers (ES, DGH or DGL). In this work the fraction of selected events is further increased by applying the three sets of criteria to each sample of simulated showers (ES, DGH, DGL) regardless of channel. With this procedure the fraction of identified Monte Carlo events is enhanced as, for instance, an ES simulated shower induced by a ν τ might not fulfill the requirements of the ES selection, but might still pass the DGH or DGL criteria, and hence contribute to the fraction of identified events. The enhancement in the fraction of events when applying this "combined" analysis depends on the particular set of Monte Carlo simulations. For instance applying the three criteria to the DGH Monte Carlo sample identifies a fraction of neutrino events ∼ 1.25 larger than when the DGH criteria are applied alone, the enhancement coming mainly from events with 3 stations rejected by the DGH criteria but accepted by ES. The application of the three criteria to the ES Monte Carlo sample however results in a smaller enhancement ∼ 1.04.

Exposure calculation
For downward-going neutrinos, once the efficiencies DG (E ν , θ, D, t) are obtained, the calculation of the exposure involves folding them with the SD array aperture and the ν interaction probability at a depth D for a neutrino energy E ν . This calculation also includes the possibility that downward-going ν τ interact with the mountains surrounding the Observatory. Integrating over the parameter space except for E ν and in time over the search periods and summing over all the interaction channels yields the exposure [18,20].
In the Earth-skimming channel, ES (E τ , θ, X d ) are also folded with the aperture, with the probability density function of a tau emerging from the Earth with energy E τ (given a neutrino with energy E ν crossing an amount of Earth determined by the zenith angle θ), as well as with the probability that the τ decays at an altitude h c [17]. An integration over the whole parameter space except for E ν and time gives the exposure [17].
The exposures E ES , E DGH and E DGL obtained for the search periods of each selection are plotted in Fig. 3 along with their sum E tot . The exposure to Earth-skimming neutrinos is higher than that to downward-going neu-trinos, partially due to the longer search period in the Earth-skimming analysis, and partially due to the much larger neutrino conversion probability in the denser target of the Earth's crust compared to the atmosphere. The larger number of neutrino flavors and interaction channels that can be identified in the DGH and DGL analysis, as well as the broader angular range 60 • < θ < 90 • partly compensates the dominance of the ES channel. The ES exposure flattens and then falls above ∼ 10 19 eV as there is an increasing probability that the τ decays high in the atmosphere producing a shower not triggering the array, or even that the τ escapes the atmosphere before decaying. At the highest energies the DGH exposure dominates. The DGL exposure is the smallest of the three, mainly due to the more stringent criteria needed to apply to get rid of the larger background nucleonic showers in the zenith angle bin 60 • < θ < 75 • .
The relative contributions of the three channels to the total expected event rate for a differential flux behaving with energy as dN ν (E ν )/dE ν ∝ E −2 ν are ES:DGH:DGL∼0.84:0.14:0.02 respectively, where the event rate is obtained as:

C. Systematic uncertainties
Several sources of systematic uncertainty have been considered. Some of them are directly related to the Monte Carlo simulation of the showers, i.e., generator of the neutrino interaction either in the Earth or in the atmosphere, parton distribution function, air shower development, and hadronic model.
Other uncertainties have to do with the limitations on the theoretical models needed to obtain the interaction cross-section or the τ energy loss at high energies. In the Earth-skimming analysis the model of energy loss for the τ is the dominant source of uncertainty, since it determines the energy of the emerging τ s after propagation in the Earth; the impact of this on the downward-going analysis is much smaller since τ energy losses are only relevant for ν τ interacting in the mountains, a channel that is estimated to contribute only ∼ 15% to the DGH exposure [18].
The uncertainty on the shower simulation, that stems mainly from the different shower propagation codes and hadronic interaction models that can be used to model the high energy collisions in the shower, contributes significantly in the ES and DG channels.
The presence of mountains around the Observatorywhich would increase the target for neutrino interactions in both cases -is explicitly simulated and accounted for when obtaining the exposure of the SD to downwardgoing neutrino-induced showers, and as a consequence does not contribute directly to the systematic uncertainties. However, it is not accounted for in the Earthskimming channel and instead we take the topography around the Observatory as a source of systematic uncertainty.
In the three channels the procedure to incorporate the systematic uncertainties is the same. Different combinations of the various sources of systematic uncertainty render different values of the exposure and a systematic uncertainty band of relative deviation from a reference exposure (see below) can be constructed for each channel and for each source of systematic uncertainty. For a given source of uncertainty the edges of the ES, DGH and DGL bands are weighted by the relative importance of each channel as given before and added linearly or quadratically depending on the source of uncertainty. In Table II we give the dominant sources of systematic uncertainty and their corresponding combined uncertainty bands obtained in this way. The combined uncertainty band is then incorporated in the value of the limit itself through a semi-Bayesian extension [21] of the Feldman-Cousins approach [22].
In the calculation of the reference exposure the νnucleon interaction in the atmosphere for DG neutrinos (including CC and NC channels) is simulated with HER-WIG [23]. In the case of ν τ CC interactions, a dedicated, fast and flexible code is used to simulate the τ lepton propagation in the Earth and/or in the atmosphere. The τ decay is performed with the TAUOLA package [24]. In all cases we adopted the ν-nucleon cross-section in [25]. In a second step, the AIRES code [26] is used to simulate the propagation in the atmosphere of the particles produced in the high energy ν interaction or in the τ lepton decay. The types, energies, momenta and times of the particles reaching the SD level are obtained. The last stage is the simulation of the SD response (PMT signals and FADC traces). This involves a modification of the "standard" sampling procedure in [27] to regenerate particles in the SD stations from the "thinned" air shower simulation output, that was tailored to the highly inclined showers involved in the search for neutrinos. Light production and propagation inside the station is based on GEANT4 [28] with the modifications to account for the evolution of the light decay-time explained above. These two latter changes roughly compensate each other, with the net result being a few percent decrease of the exposure with respect to that obtained with the standard thinning procedure and a constant average value of the light decay-time.

IV. RESULTS
Using the combined exposure in Fig. 3 and assuming a differential neutrino flux dN (E ν )/dE ν = k · E −2 ν as well as a ν e : ν µ : ν τ = 1 : 1 : 1 flavor ratio, an upper limit on the value of k can be obtained as: Total ∼ +37%, -28% Table II. Main sources of systematic uncertainties and their corresponding combined uncertainty bands (see text for details) representing the effect on the event rate defined in Eq. (1). The uncertainty due to "Simulations" includes: interaction generator, shower simulation, hadronic model, thinning and detector simulator. The uncertainty due to "τ energyloss" affects the ES channel and also the DGH but only to ντ with θ 88 • going through the mountains surrounding the Pierre Auger Observatory. However it does not affect the DGL channel. The topography around the Observatory is not accounted for in the ES channel and is taken as a systematic uncertainty that would increase the event rate.   Table I to Monte Carlo simulations of UHE neutrinos (see text for explanation). Also shown are the individual exposures corresponding to each of the three selections. For the downward-going channels the exposure represents the sum over the three neutrino flavors as well as CC and NC interactions. For the Earth-Skimming channel, only ντ CC interactions are relevant.
The actual value of the upper limit on the signal events (N up ) depends on the number of observed events (0 in our case) and expected background events (conservatively assumed to be 0), as well as on the confidence level required (90% C.L. in the following). Using a semi-Bayesian extension [21] of the Feldman-Cousins approach [22] to include the uncertainties in the exposure we obtain 3 N up = 2.39. The single-flavor 90% C.L. limit is: k 90 < 6.4 × 10 −9 GeV cm −2 s −1 sr −1 . (3) The limit applies in the energy interval ∼ 1.0 × 10 17 eV − 2.5 × 10 19 eV where the cumulative number of events as a function of neutrino energy increases from 5% to 95% of the total number, i.e. where ∼ 90% of the total event rate is expected. It is important to remark that this is the most stringent limit obtained so far with Auger data, and it represents a single limit combining the three channels where we have searched for UHE neutrinos. The limit to the flux normalization in Eq. (3) is obtained integrating the denominator of Eq. (2) in the whole energy range where Auger is sensitive to UHE neutrinos. This is shown in Fig. 4 , along with the 90% C.L. limits from other experiments as well as several models of neutrino flux production (see caption for references). The denominator of Eq.
(2) can also be integrated in bins of energy, and a limit on k can also be obtained in each energy bin [35]. This is displayed in Fig. 5 where the energy bins have a width of 0.5 in log 10 E ν , and where we also show the whole energy range where there is sensitivity to neutrinos. The limit as displayed in Fig. 5 allows us to show at which energies the sensitivity of the SD of the Pierre Auger Observatory peaks. The search period corresponds to an equivalent of 6.4 years of a complete Auger SD array working continuously. The inclusion of the data from 1 June 2010 until 20 June 2013 in the search represents an increase of a factor ∼ 1.8 in total time quantified in terms of equivalent full Auger years with respect to previous searches [17,18]. Further improvements in the limit come from the combination of the three analysis into a single one, using the procedure explained before that enhances the fraction of identified neutrinos especially in the DGH channel.
In Table III we give the expected total event rates for several models of neutrino flux production.
Several important conclusions and remarks can be stated after inspecting Figs. 4 and 5 and Table III: 1. The maximum sensitivity of the SD of the Auger Observatory is achieved at neutrino energies around EeV, where most cosmogenic models of ν production also peak (in a E 2 ν × dN/dE ν plot). 2. The current Auger limit is a factor ∼ 4 below the Waxman-Bahcall bound on neutrino production in optically thin sources [13]. The SD of the Auger Observatory is the first air shower array to reach that level of sensitivity.
3. Some models of neutrino production in astrophysical sources such as Active Galactic Nuclei (AGN) are excluded at more than 90% C.L. For the model yields a value of Nup = 2.39 slightly smaller than the nominal 2.44 of the Feldman-Cousins approach.
#2 shown in Fig. 14 of [32] we expect ∼7 neutrino events while none was observed.
4. Cosmogenic ν models that assume a pure primary proton composition injected at the sources and strong (FRII-type) evolution of the sources are strongly disfavored by Auger data. An example is the upper line of the shaded band in Fig. 17 in [31] (also depicted in Figs. 4 and 5), for which ∼4 events are expected and as consequence that flux is excluded at ∼98% C.L. Models that assume a pure primary proton composition and use the GeV γ-ray flux observations by the Fermi-LAT satellite detector as an additional constraint, are also disfavored. For instance for the model shown as a     Waxman-Bahcall '01 [13] Figure 5. Upper limit to the normalization of the diffuse flux of UHE neutrinos (at 90% C.L. and in bins of width 0.5 in log 10 Eν -see text for details) from the Pierre Auger Observatory (straight steps). We also show the corresponding limits from ANITAII [29] (dot-dashed line) and IceCube [30] (dashed line) experiments (with appropriate normalizations to take into account the energy bin width, and to convert to single flavor), along with expected fluxes for several cosmogenic neutrino models [9,31,33] as well as the Waxman-Bahcall bound [13] (all converted to single flavor).
solid line in the bottom right panel of Fig. 5 in [33] (also depicted in Figs. 4 and 5 in this work), corresponding to the best-fit to the cosmic-ray spectrum as measured by HiRes, we expect ∼3.2 events. As a consequence that model is excluded at more than 90% C.L. For this particular model we also show in Figs. 4 and 5 the 99% C.L. band resulting from the fitting to the HiRes spectrum down to E min = 10 19 eV. The Auger limit is also approaching the solid line in the upper left panel of Fig. 5 in [33], a model that assumes extragalactic protons above E min = 10 17.5 eV [36], for which ∼ 1.6 events are expected (see Table III). The Auger direct limits on cosmogenic neutrinos are also constraining part of the region indirectly bounded by Fermi-LAT observations. 5. The current Auger limit is less restrictive with the cosmogenic neutrino models represented by the gray shaded area in the bottom panel of Fig. 4 (∼0.5 to ∼1.4 events are expected as shown in Table III) which brackets the lower fluxes predicted under a range of assumptions for the composition of the primary flux (protons or mixed), source evolution and model for the transition from Galactic to extragalactic cosmic-rays [9] The same remark applies to models that assume pure-iron composition at the sources. A 10-fold increase in the current exposure will be needed to reach the most optimistic predictions of cosmogenic neutrino fluxes if the primaries are pure iron, clearly out of the range of the current configuration of the Auger Observatory. 6. A large range of exotic models of neutrino production [34] are excluded with C.L. larger than 99%.

V. ACKNOWLEDGMENTS
The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort