#### Abstract

The detection of ultrahigh-energy neutrinos, with energies in the PeV range or above, is a topic of great interest in modern astroparticle physics. The importance comes from the fact that these neutrinos point back to the most energetic particle accelerators in the Universe and provide information about their underlying acceleration mechanisms. Atmospheric neutrinos are a background for these challenging measurements, but their rate is expected to be negligible above ≈1 PeV. In this work we describe the feasibility to study ultrahigh-energy neutrinos based on the Earth-skimming technique, by detecting the charged leptons produced in neutrino-nucleon interactions in a high mass target. We propose to detect the charged leptons, or their decay products, with the High Altitude Water Cherenkov (HAWC) observatory and use as a large-mass target for the neutrino interactions the Pico de Orizaba volcano, the highest mountain in Mexico. In this work we develop an estimate of the detection rate using a geometrical model to calculate the effective area of the observatory. Our results show that it may be feasible to perform measurements of the ultrahigh-energy neutrino flux from cosmic origin during the expected lifetime of the HAWC observatory.

#### 1. Introduction

The first evidence of ultrahigh-energy neutrinos (in the PeV energy range) from extraterrestrial origin was recently reported [1]. This opened a new field in astroparticle physics that will allow the identification and characterization of the most powerful particle accelerators in the Universe. Neutrinos are not affected by the electromagnetic or strong interactions and thus point back to the source where they were produced, unlike charged cosmic rays. Gamma-rays are another cosmic probe that provides information about the acceleration mechanisms that occur in astrophysical sources. Due to this, there are several dedicated instruments, both ground or space based, performing a continuous survey of the Universe characterizing gamma-ray sources, for example, the Fermi Gamma-Ray Space Telescope and the imaging atmospheric Cherenkov telescopes HESS, MAGIC, VERITAS, and FACT.

The HAWC observatory was designed to detect and characterize the sources of high-energy gamma-rays in the energy range between 100 GeV and 100 TeV [2] and started full operations in April 2015. It is a ground based instrument that detects atmospheric showers by measuring with high precision the arrival time of the particles that compose the air showers. This is done via the Cherenkov light produced by the air shower particles as they enter the 300 water detector tanks that constitute the observatory, with a total water volume of 54 million litres. At energies above 100 TeV, gamma-rays suffer strong absorption from pair production with photons from the cosmic microwave background radiation that strongly reduces their mean free path [3]. For this reason, ultrahigh-energy neutrinos may be a better tool to study the most energetic extragalactic particle accelerators.

The Earth-skimming technique to detect ultrahigh-energy neutrinos has been proposed before; see, for instance, [4–8]. The method consists of using the interaction between a neutrino and a nucleon via the exchange of a boson to produce a charged lepton of the same flavour as the incoming neutrino. Since the neutrino-nucleon cross section is very small, a large-mass target is needed. The natural candidate to produce such interactions is the Earth crust, either by focusing the searches on quasi-horizontal neutrinos that travel along a chord inside the Earth or by using mountains as targets. The produced charged lepton travels essentially in the same direction as the neutrino. Thus, in the considered scheme, the charged lepton will travel either upward if it was moving through a chord inside the Earth crust or quasi-horizontally if the neutrino passed through a mountain.

From the three families of leptons, the electron neutrinos are unfavourable for this type of study because the produced high-energy electrons initiate electromagnetic showers that are easily absorbed by the target mass shortly after production. For *μ* neutrinos, the relatively large mean life of the produced charged lepton combined with the ultrahigh energy will produce a detectable signal only as an ultraenergetic . The neutrinos are the ones that have attracted more interest from the experimental point of view. The reason is that because of their very short mean life, even if they are very energetic the produced charged leptons will decay into secondary particles that would make the detection of the signal easier. However, the neutrinos have proven to be one of the most elusive particles of the standard model, with less than 15 detections up to now [13]. Ultrahigh-energy neutrinos (with searches up to 72 PeV) have eluded direct detection so far, even after the analysis of three years of IceCube data [14]. Even though neutrinos are disfavoured in production mechanisms at the astrophysical acceleration sites, the neutrino flavour mixing that occurs in cosmological distances is expected to produce approximately equal proportions of all neutrino flavours at the Earth.

There have been already experimental attempts to detect ultrahigh-energy neutrinos using the Earth-skimming technique. For instance, the Pierre Auger observatory in Argentina used their surface detectors to look for the electromagnetic signature of extensive air showers initiated by the decay of charged leptons of EeV energies that develop close to the detector [15–17], without finding candidate signals so far. There are also studies that propose to use the fluorescence detector of the Pierre Auger observatory to detect the decay in the atmosphere of charged leptons produced by ultrahigh-energy neutrinos [18]; however, this idea has not been implemented yet. The Ashra-1 collaboration [19] searched for neutrino emission from a GRB in the PeV–EeV energy range using the Earth-skimming technique. The Ashra-1 experiment, located on the Mauna Loa volcano and facing the Mauna Kea volcano in the Hawaii island, aimed to detect neutrinos that converted into charged leptons inside Mauna Kea. Their method consisted of measuring the Cherenkov light emitted by the particles of the atmospheric shower initiated by the decay products of the charged lepton. Their analysis of the GRB081203A did not find signals associated to neutrinos in the PeV–EeV energy range. There have been also studies about the feasibility to use the MAGIC telescopes for the detection of neutrinos [20, 21], by pointing their telescopes below the horizon towards the sea, or by searching for reflections of Cherenkov light by the nearby ground, the sea, or clouds [22, 23]. However, no experimental results have been published yet.

In this paper we propose to adapt the Earth-Skimming technique to use it with the HAWC gamma-ray observatory, employing the Pico de Orizaba volcano as a target for the neutrino-nucleon interactions. However, we propose not to follow the method explored so far of studying the decay products of a charged lepton on the atmosphere. Instead, we propose to reconstruct directly the trajectory of the charged lepton or their decay products as they travel through the HAWC detectors. In this way we do not only restrict our studies to neutrinos but also include the possibility for the detection of ultrahigh-energy neutrinos (a first suggestion of this method was mentioned in [11]). At ultrahigh energies, the produced charged leptons will have an energy approximately equal to that of the original neutrino [7, 11]. These charged leptons, or their decay products, could travel crossing several HAWC detectors depositing large amounts of Cherenkov light, well beyond the average left from both atmospheric muons (≈30 photoelectrons (PEs) [24]) located far from the shower core used in the gamma/hadron discrimination algorithms of HAWC and the high PE noise, also associated to atmospheric muons, considered by the HAWC collaboration to be in the range of 10–200 PEs [25]. Here it is important to point out that the dynamical range of the HAWC electronics goes from a fraction of a PE up to thousands PEs [26].

In order to demonstrate that our proposal to search for tracks produced by charged leptons or their boosted decay products is feasible, we performed some GEANT4 [27] simulations of the decay of charged leptons with an energy of 1 PeV (that decay approximately 50 m after their creation point) and studied the shower evolution in air. We choose an energy of 1 PeV so the charged leptons decay quickly and are in the energy regime of our studies. After decay, the opening angle of the charged products is smaller than 0.7°, thus, after the decay products have travelled 2050 m (the approximate distance from the edge of the volcano to the HAWC array at ≈4090 m a.s.l. is of two kilometers), they would hit at most three columns of HAWC tanks. Moreover, over 95% of the secondary charged particles are contained within an opening angle smaller than 0.2°, therefore producing large Cherenkov signals only within a single row of HAWC tanks, producing a clearly identifiable track. This simplified exercise was done for 1 PeV charged ’s, making it easy to extrapolate the results to higher energy charged leptons, since the opening angle of the decay products is inversely proportional to the energy of the primary charged lepton [28].

Based on this information, we believe that the tracking method is possible in the search for charged leptons, produced by neutrino-nucleon interactions, in the PeV energy range. The amount of Cherenkov light that could be detected by the HAWC PMTs by such energetic particles is expected to be in the range of thousands to tens of thousands PEs. However, given the dynamic range of the current HAWC electronics one could anticipate that the collected light could be used as a proxy to at least set a lower boundary on the energy of the incoming charged lepton, in a similar manner to what IceCube does for muons that pass through the detector.

The paper is organized as follows: in Section 2 we present the calculation of the flux of charged leptons that could be produced by neutrino-nucleon interactions in the Pico de Orizaba volcano. In Section 3, we describe a method to calculate the HAWC effective area based on purely geometrical considerations and evaluate it for different trigger conditions. Then, in Section 4, we present our results for the possible detection rate and address the issue of the expected background signals. In Section 5, we discuss how the current trigger of HAWC can be useful in selecting data for these studies and discuss further a possible background rejection strategy. Finally, the conclusions are presented in Section 6.

#### 2. Calculation of the Flux of Charged Leptons Produced by Earth-Skimming Neutrinos

In order to obtain an estimate of the number of ultrahigh-energy charged leptons that could be produced via the Earth-skimming technique, we follow the formalism developed in [7, 29]. However we use further simplifications due to the detection method that we propose in this paper and also because of the energy of the neutrinos that we plan to study. The differences between the original formalism of [7, 29] and our implementation are pointed out in the text.

The number of charged leptons () that could be detected by the observatory is given bywhere is the charged lepton flux, is the effective area of the observatory, is the live time of the experiment, and is the duty cycle for observations, which basically describes which fraction of the time the experiment is actually able to take data. In this section we describe the calculation of the flux of charged leptons produced by the neutrino-nucleon interactions that occur while the neutrinos traverse the volcano. We start with the differential flux of ultrahigh-energy neutrinos produced by astrophysical sources, which we consider to be isotropic. This differential flux is given bywith () being, respectively, the energy, polar, and azimuthal angle of the neutrinos. Since we are interested in integrating the flux on only a certain region (the one covered by the Pico de Orizaba volcano) and because we are considering an isotropic flux, we can simplify the differential neutrino flux towhere is the solid angle covered by the volcano that is being used as the target for the neutrino-nucleon interactions. After the proposed isotropic neutrino flux passes through the Pico de Orizaba volcano, the produced differential charged lepton flux is given bywith () being, respectively, the energy, polar, and azimuthal angle of the produced charged leptons. The relation between the differential fluxes of incoming neutrinos and the produced charged leptons is thus given bywhere is a function that physically represents the convolution of the probabilities of the different processes that need to take place in order that an ultrahigh-energy neutrino converts into a charged lepton inside the target material and is able to escape the mountain. Given that we are interested in studying ultrahigh-energy neutrinos (), the first simplification comes from the fact that the produced charged lepton will approximately follow the same direction as the original neutrino. The angle between the original neutrino and the produced charged lepton () has been estimated to be smaller than 1 arcmin for energies above 1 PeV for ’s [19]; so, we expect this angle to be negligible. This makes that the angular dependence of can be approximated by delta functions.

The energy dependent part of the function can be written as the integral along the path of the neutrino and the corresponding charged lepton inside the volcanowhere is the survival probability for a neutrino travelling a certain distance inside the volcano. This probability can be written aswhere is the charged current interaction length. Since the average width of the volcano () is much smaller than the interaction length, almost all of the neutrinos will traverse the whole mountain. That defines the integration limit of (8). The interaction length is a function of the density of the medium that the neutrino travels through. In the case of the formalism developed in [7], where the trajectories of the neutrinos were across chords inside Earth, the authors had to consider variations of the Earth density. In this particular case, since we are interested in neutrinos that travel through volcanoes, it is reasonable to consider a constant density along the path of the neutrinos. Thus, the expression for the charged current interaction length can be written as follows:where is Avogadro’s constant and is the charged current cross section. The second probability () that enters the calculation of is that of a neutrino to convert into a charged lepton through a charged current in an infinitesimal distance . This probability is given by

The third process () that is taken into account in the function is the probability that the produced charged lepton is able to escape the volcano, taking into account the charged lepton energy losses in the medium. This process is described by a system of two coupled differential equations:

Where (11) describes the energy loss processes, with describing the ionization energy loss and the radiative energy loss. According to the literature, for example, [7], the effects of are negligible at the energy regime of interest of this work; so, we consider in (11) that . In (12), is the mass of the charged lepton, is the speed of light, and is the charged lepton lifetime. A Monte Carlo study from [11] shows that, for ’s with energies of 1 PeV, after 6 km water equivalent (km.w.e) (6 km.w.e. ≈ 2.3 km in “standard rock,” the average path length inside the volcano), that is, in the extreme case where the charged lepton has to travel the average width of the volcano, and this value approaches unity as the energy of the charged lepton increases. Thus, we consider that, for the energy regime studied in this work, it is appropriate that for ’s we can take a value of (see Table 1).

The case of the survival probability for charged leptons is more complicated to evaluate, since most of the research has been done for energies above or at 100 PeV [30–32]. We take as a base for our calculations the results presented in [31]. For , and, for , after 2.3 km in “standard rock,” that is, in the extreme case where ’s are produced just after entering the volcano. However, the values of at lower energies cannot be obtained by a simple extrapolation since the value of below 100 PeV is dominated by the mean life time of . By considering this fact, we obtain approximate values of for the PeV energy regime. Table 1 shows the average values of in the energy range of interest of this work. Finally, the factor makes sure that the charged lepton escapes the volume of the volcano with an energy . Based on (11), taking , then can be written as

Thus, we can approximate the function as

After evaluating the integral we obtain an expression for the function.

Equation (15) is valid in the case where we neglect variations of the density of the volcano. Going back to (5), we can substitute (6) and (15) into it. By doing this, we obtain an equation for the differential flux with respect to the energy of the leptons

The next step is to define a function that describes the differential flux of neutrinos as a function of energy. For this, we use the parametrization of the measured astrophysical neutrino flux made by IceCube [12, 33, 34], which has been found to be well described by an unbroken power law

We choose to use the parametrization presented in [33] for the measurement of the astrophysical flux. For this particular parametrization we shifted the mean values of the normalization and spectral index in order to obtain the highest flux, within the allowed range given by their statistical uncertainties. These values are and . We assume for our calculations an expected equal contribution to the astrophysical flux for , due to the neutrino flavour mixing over cosmological distances that would produce approximately equal proportions of all flavours. Then, we extrapolate the measured flux to the energy range 10 PeV to 100 EeV. Moreover, motivated by the most energetic event found by the diffuse flux muon neutrino search (2009–2015) done by IceCube (a track event that deposited in the sensible volume of the detector [35]), we also include in the flux estimations the models proposed in [9], which account for sources of multi-PeV neutrinos that are constrained by the most recent ultrahigh-energy neutrino upper limits set by IceCube [36] and Pierre Auger [15]. These models, for the sum of all neutrino flavours, have the smoothly broken power law functional forms:

The values of the parameters are , and . take the values of 10^{7} GeV for what we refer in this work as Model A, 10^{8} GeV for Model B, and 10^{9} GeV for Model C. The values of the normalizations for each of these models are presented in Table 2.

According to [9], Models A and B would correspond to the spectra produced by BL Lac AGNs, and combinations of Models A and C would follow the expected shape of GZK neutrinos produced from EBL and CMB interactions.

In order to calculate the flux of charged leptons produced by neutrino-nucleon interactions we need to define the input parameters that enter the calculation of the number of produced charged leptons (see (16)). For the value of the parameter that describes the radiative energy losses of the charged leptons as they travel through the rock (), one can find in the literature parametrizations that may differ by up to a factor of two. For instance, it ranges from 0.26 to at 100 PeV in [10]. For our calculations we decided to use the values obtained with the ASW structure functions calculated in [10] and for the results obtained in [11]. For this latter case the results are available up to an energy of 1 EeV, but these values are enough for our calculations. For the numerical integrations, we used intermediate values of and in different bins of energy as shown in Table 3.

For the density of the volcano, we take that of “standard rock,” [11]. For the calculation of the charged current interaction length we use the results from [37], with the cross section for neutrino-nucleon interactions via charged currents in the energy range given by

This parametrization is a well known result. However, as a cross check, we compared the parametrization given by [37] with more recent calculations [38–41]. The result of this is presented in Figure 1. One can notice that the result from [37] follows the general trend predicted by the work of [38] and Sarkar et al. Due to this and the fact that there is an analytical expression for the cross section of [37], we decided to use (19) in our numerical calculations. The anti neutrino-nucleon cross section is taken to be the same [37]. Using these results we get a value of , much smaller than the light-year interaction length for average neutrinos in lead [42].

For the average width of the Pico de Orizaba, we calculated the typical length of a chord that goes through a cone that follows the geometry of the volcano, as will be shown in the following section of the paper. At 4100 m a.s.l., the volcano has a width of ≈6 km and at the summit at ≈5,500 m a.s.l., a width of ≈0.35 km (see Figures 4 and 5). The average path length for the neutrinos is of 2.33 km inside the volcano. Finally, we assume a duty cycle , which agrees with the reported value observed during actual HAWC operations [43]. The energy range for the charged leptons that exit the volcano is taken in the range (, ). In this way we consider an energy range that contains a good fraction of the charged leptons that exit the volcano. This is because the charged leptons escape with an energy distribution, after loosing some of their energy in the medium. For a study of the propagation of a monoenergetic beam of charged leptons through rock see, for instance, the work from [31].

We integrate numerically (16) using the differential neutrinos fluxes from (17) and (18) to obtain the number of charged leptons () that escape the mountain. The results are presented in Tables 4 and 5.

One can notice that the flux of charged leptons is more than an order of magnitude lower than that of ’s. The reason for this can be seen in (16), where the charged lepton flux depends inversely on the radiative energy losses of the charged leptons, which are more than an order of magnitude higher for ’s with respect to ’s (see Table 3). This does not contradict the fact that, at the highest energies, for both ’s and ’s the survival probability . quantifies the probability that the charged leptons are able to escape the volcano, while the factor in (16) appears because of the form of and is independent of . This result is consistent with the arguments developed in [44], where it is pointed out that charged leptons have a higher probability to escape the Earth crust compared to ’s.

In the following section we present a simple method that allows to approximate the effective area of the HAWC observatory to ultrahigh-energy charged leptons.

#### 3. Effective Area Determination

In general, the effective area or acceptance of an observatory is calculated using a detailed Monte Carlo simulation; see, for instance, in [45, 46] studies of the sensitivity of air Cherenkov and fluorescence observatories to neutrinos using the Earth-skimming technique and in [22] an early proposal about using fluorescence measurements on observatories as Pierre Auger in order to search for charged lepton showers. The full Monte Carlo method of calculating the effective area represents a complex and time consuming task. In the early results obtained by the Pierre Auger collaboration, their effective area calculation neglected the effect of the topography that surrounds the observatory and took into account its effect in the systematic error of their observation limits [17]. However, their most recent results have incorporated the topography into their effective area calculations in two of their three analysis channels [15]. In our case we are specifically interested in the effect of the largest volcano that surrounds the HAWC observatory as a target for the neutrino-nucleon interactions. In order to reproduce the topography that surrounds HAWC, we use data from the* Instituto Nacional de Estadística y Geografía* (INEGI, México) [47]. Figure 2 shows the topography that surrounds the HAWC site (indicated by a small red rectangle). The observatory is located between two volcanoes: in the direction Northing-Easting by the Pico de Orizaba and in the opposite direction by the Sierra Negra volcano. Due to its much larger volume, we focus our attention on the Pico de Orizaba volcano as the target for the neutrino-nucleon interactions.

Figure 3 shows the topography of the Pico de Orizaba volcano as seen in a coordinate system centered at the location of the HAWC array at ≈4090 m a.s.l. The approximate symmetry of the Pico de Orizaba volcano makes it easy to motivate a simplification in the calculation of the effective area. We can assume that the geometry of the mountain can be modelled to be conical. Figures 4 and 5 show the profile of the volcano along the Easting and Northing directions, and how we can approximate this profile using a cone with a height of 1500 m ( taken at the altitude above sea level of the HAWC observatory) and a diameter of 6000 m.

Our method to calculate the effective area to detect ultrahigh-energy charged leptons produced by neutrino-nucleon interaction is as follows:(1)We consider that the ultrahigh-energy charged leptons, or their highly boosted decay products, will travel following a straight trajectory along the direction of the initial neutrino.(2)We approximate the region where the charged leptons are produced as a triangular surface, a 2D projection of the volcano, that faces the HAWC observatory. We can use this simplification since here we are only interested in studying the trajectories of the produced charged leptons. The width of the mountain is considered in the charged lepton flux calculation, as presented in Section 2. On the triangular surface we draw a grid made of lines parallel to the - and -axis in the coordinate system shown in Figure 6. The lines of the grid have a separation of 0.5 m in each direction. The surface is located at a distance of 5.7 km to the center of the HAWC observatory, approximately the distance between the center of the detector array and that of the volcano.(3)The detection volume of HAWC is modelled as a rectangular prism, with dimensions () 140 m × 140 m × 4.5 m. This approximation is motivated by the actual configuration of the HAWC array. The observatory is made of a compact group of 300 Water Cherenkov Detectors (WCDs). Each WCD is a cylinder 7.3 m in diameter and 4.5 m in height. The array covers an area of approximately 22,000 m^{2} [43].(4)From each cell of the grid (~14 million), we generate vectors that point towards our model of HAWC, approximated as a rectangular prism. Figure 7 shows a diagram of the definition of the angles used in this work. The azimuth angle covers the range from to (with steps of 1 degree), pointing towards the direction in the coordinate system shown in Figure 6. This is done for the different orientations of the polar angle . The polar angle is defined such that the vector orientation that lies in the plane corresponds to 90° and increases as the orientation gets closer to the plane (see Figure 7).(5)For each orientation , we calculate a differential element of effective area where is the total number of vectors whose directions points towards the detection volume of HAWC with a minimum length of the trajectory across the plane shown in Figure 6 (inside the volume of the rectangular prism). We will refer to this minimum length as the trigger condition that approximately represents how many WCDs would measure Cherenkov light produced by the incoming lepton or their collimated decay products. is the total number of vectors that point towards the rectangular prism, regardless of the trajectory that they have inside the prism. is the section of the area of the triangular surface over which the vectors that point towards the rectangular prism () were generated. Finally, is the azimuthal angle covered by the area , as seen from the center of the rectangular prism.

By following the procedure described above, we can obtain the differential elements of effective area as a function of the polar angle orientation. Figure 8 shows the results for different orientations of the polar angle.

One can see that the differential effective area increases with increasing values of the polar angle, up to a maximum at ≈93°, and then decreases because the area over which the vectors are generated decreases as well, as it corresponds to the upper sections of the triangular grid where the vectors are generated. The different line styles used in the plot show the results obtained for different trigger conditions. For instance, the curve with 14 m as trigger condition indicates the results obtained when the requirement of the path of the charged lepton through the detection volume corresponds to at least a trajectory in the plane inside the HAWC rectangular prism of 14 m. This would correspond roughly to having the ultrahigh-energy charged lepton (or its boosted decay products) to traverse through at least two WCDs. This may not seem enough to even provide a rough idea of the direction of the incoming lepton. However, one should keep in mind that each WCD is equipped with four photomultiplier tubes (PMTs), recording the amplitude and arrival time of the Cherenkov light. By combining the information from the eight PMTs from two WCDs it may be possible to have enough information to reconstruct a trajectory, although details of the angular reconstruction of the tracks are beyond the scope of this work. In Figure 8, we present the results for trigger conditions that correspond roughly to having the charged lepton (or its collimated decay products) traversing from two up to six WCDs. As expected, the effective area that corresponds to a given trigger condition decreases as one increases the number of detectors required to have signals. We are interested in the total effective area that goes into the calculation of the number of possible detections of (1). This total effective area is given bywhere the integration range in goes from 90° to 105° (the angular height of the Pico de Orizaba volcano as seen from the center of the HAWC observatory is of ≈15°). Our results of the effective area calculation as a function of the different trigger conditions are shown in Table 6.

With the results from Sections 2 and 3 we have the flux of charged leptons given a certain isotropic flux of neutrinos and the effective area of the HAWC observatory obtained with a simple trigger condition. Using this, in the next section we present our results for the flux of ultrahigh-energy charged leptons that could be detected by the HAWC observatory.

#### 4. Results, Possible Background Signals, and Discussion

Table 7 shows the number of detectable charged ’s produced by neutrino-nucleon interactions in the Pico de Orizaba volcano, obtained using the effective area results presented in Table 6 and the ultrahigh-energy lepton fluxes from Table 4. We restrict the results to those from charged leptons since, as it was shown in Section 2, the flux from charged leptons is approximately an order of magnitude lower, making their detection not feasible with the method proposed in this work. The detection estimates are shown as a function of the different trigger conditions and of the different ultrahigh-energy neutrino fluxes discussed in Section 2. The eighth column presents an estimate of the background signals that are expected using this method. We estimated the background considering that the main contribution is that of ultrahigh-energy atmospheric muons, coming from the direction of the Pico de Orizaba. In order to calculate this, we used the characterization of the atmospheric muon flux above 15 TeV measured by IceCube [48], that can be modelled as an unbroken power law:

In this case, we shifted the mean values (within the statistical and systematic uncertainties of the IceCube measurement) of the normalization and spectral index in order to obtain the worst case scenario of the background and then extrapolated the muon flux to the energy range of our interest. The values of the parameters that maximize the background muon flux are and a spectral index of 3.73. Thus, the background flux can be calculated using (1), where in this case the charged lepton flux corresponds to the integral of (22) within the energy bins used in this work. Of course this is an overestimation, since these atmospheric muons would have to be able to survive their travel through the volcano; however, we consider this estimation of the background appropriate because the study of the propagation of atmospheric muons through the volcano is beyond the scope of this work. Moreover, it has already been noticed that the atmospheric muon background is relevant up to energies of ≈100 TeV [49]. The last two columns of Table 7 show the estimated time to get a detection for each trigger condition, first for the case of only having as a source the neutrino flux given by the extrapolation of the measurement done by IceCube [12] and also for the case of including GZK neutrinos from [9].

From Table 7, we find that it would take around eight years of data taking in order to be able to find a track-like signal traversing three WCDs based on the extrapolation of the measured astrophysical flux, for neutrinos in the energy range from 10 PeV to 100 PeV. In order to find a track-like signal that propagates through four WCDs it could take more than nine years of data taking, but such a long track may allow reconstructing with better accuracy the direction of the original neutrino. We have verified the predicted detection rate by doing a back of the envelope comparison with the PeV neutrino measurements reported by IceCube [1] and found that, taking into account the different experimental conditions, our results fall in the expected order of magnitude.

The HAWC observatory is currently planned to operate for at least five years [50], but given its importance as a trigger for future instruments such as CTA, it is not unlikely that it may operate for at least 10 years. Table 7 shows also for completeness the results expected for neutrinos in the energy range from 100 PeV to 100 EeV. Charged leptons that do not decay before arriving or inside the HAWC array (energies of ≈100 PeV or larger) will leave almost the same signal because the values of (see Table 3) change very slowly with energy [8]. This condition will make it very difficult to estimate the energy of the incoming . The results also show that the background from ultrahigh-energy muons at the energy regime of interest of this work is very low and would allow a clean identification of a signal coming from the direction of the Pico de Orizaba volcano.

#### 5. Triggering of Signals and Background Rejection

The current software trigger used by the HAWC observatory consists of requiring a certain number of PMTs to be above a charge threshold (it actually consists of two thresholds used to estimate the charge with the Time-over-Threshold method; the thresholds are of 0.25 and 4 PEs) within a time sliding window of the order of ≈150 ns [24]. For the HAWC results presented in [43], the number of PMTs required to trigger an event was of 15. Thus, the multiplicity trigger used by the HAWC collaboration is already useful for the collection of signals that are needed for the neutrino searches proposed in this work. To show this, we can estimate the average number of PMTs that have signals during any 150 ns time window. From [51] we know that the sum of the signals from 112 HAWC’s PMTs in a 60 s time window was of ≈2.285 × 10^{8}. Then, we can use this number to estimate the number of PMTs fired in any 150 ns triggering window, which is of ≈6. Thus, only 9 additional PMTs would be required to have signals within the specified time window in order to produce a trigger useful for the neutrino search proposed in this work. This condition can be easily fulfilled if an ultrahigh-energy lepton or its decay products pass through 3 WCDs, firing at least nine of its 12 PMTs (with a propagation time of the corresponding signals along 3 WCDs of ≈70 ns, well within the current trigger window). Of course, a specialized trigger can be developed in order to obtain a better trigger efficiency. This could be done, for instance, by selecting events with topologies consistent with tracks propagating through the WCDs. This would also have the consequence of reducing the bandwidth needed for this additional trigger.

A second point of interest to discuss is the ability of the instrument to separate the signals produced by ultrahigh-energy charged leptons from those from the background of lower energy atmospheric muons. The amount of emitted Cherenkov photons per unit length inside a HAWC WCD is given by the Frank-Tamm equation; and the total amount of light will be proportional to the the sum of the path lengths of the primary and all the produced secondary charged particles (from bremsstrahlung and pair production processes) that travel inside a WCD and have a velocity larger than the threshold for Cherenkov light production. According to the study presented in [52] this sum of path lengths increases linearly with the energy of the primary particle. Thus, the way to discriminate background signals could be based on the total light yield detected in the WCDs that belong to a given track (e.g., a 1 TeV primary electron would produce an order of magnitude more photons than those from a 100 GeV primary electron). Based on the results presented in [53], the mean momentum of nearly horizontal tracks () is of ≈200 GeV/c, and in this work we aim to detect the tracks produced by leptons of PeV energies, so we could expect a factor larger than 10^{3} in the deposited light in the WCDs for the ultrahigh-energy neutrino initiated signals relative to the average muon background tracks. Moreover, the results presented in [54] show that the integral intensity of muons of ≈1 TeV/c is larger than that of ≈200 GeV/c muons by a factor of ≈52, making it reasonable to easily handle the data rate if a specialized trigger is set to keep the data of tracks related to large deposits of Cherenkov light.

An additional experimental proof that this analysis could be performed, is the ability to observe the cosmic ray shadow produced by the volcano, as, for instance, pointed out in [22]. However, this requires the use of data and algorithms property of the HAWC collaboration that are not publicly available. Nonetheless we are confident, because of the interest on this result, that the Pico de Orizaba cosmic ray shadow will be publicly available soon.

Due to these facts, we believe that it is feasible to separate the huge atmospheric background from the signals we are looking for. We are aware that the precision that could be achieved in determining the energy of the primary tau lepton will be low; nonetheless the primary interest of this work is the detection technique and we leave the development of a possible energy estimator for future work.

#### 6. Conclusions

We presented an estimate of the detection capabilities of the HAWC observatory to study ultrahigh-energy neutrinos interacting in the Pico de Orizaba volcano. We based our study in the analytic method developed by [7, 29] and modified it to a simpler case where the neutrino conversion takes place inside a medium of constant density. The effective area of the HAWC observatory for the detection of ultrahigh-energy charged leptons was calculated geometrically for different trigger conditions. We used the astrophysical neutrino flux measured by IceCube [12] and extrapolated it to energy range from 10 PeV up to 100 EeV and also considered models for multi-PeV neutrinos [9] that are constrained by the most recent data from both IceCube [36] and the Pierre Auger observatory [15]. With this, we found that finding a signal consistent with the propagation of an ultrahigh-energy charged lepton coming from the direction of the Pico de Orizaba volcano will require approximately nine years of data taking if the signal is required to propagate in four WCDs of HAWC. We estimated the expected background for this analysis using the measured atmospheric muon flux above 15 TeV [48] and found that it should be feasible to perform this study with a reasonable signal to background ratio during the lifetime of the HAWC observatory. We also showed that the current software trigger used by HAWC should be sufficient to acquire data for this analysis and that given the dynamic range of the HAWC electronics it is feasible to be able to discriminate the background signals produced by lower energy muons using the light yield detected by the HAWC PMTs.

As an anonymous referee pointed out to us, the detection rate should be taken with caution since there are several quantities that are uncertain to at least some degree, such as the flux of neutrinos at the ultrahigh-energy regime, the neutrino-nucleon cross section, and the parameters that describes the radiative energy loss processes for the leptons traversing the Earth crust. By selecting a different set of the input parameters for the calculations one could find a detection rate that could decrease by a factor of three or even more. However we also did not consider scenarios where the neutrino-nucleon cross section at ultrahigh energies could have enhanced values relative to the standard model predictions, due to the presence of new physics, as described, for instance, by [55]. Further work should be done using a complete Monte Carlo simulation that incorporates the detailed topography of the volcanoes that surround the HAWC observatory and the complete detector simulation to study the detector response. However, this first step shows encouraging results to pursue more detailed studies. The detection rate is certainly low, not comparable to the one that dedicated neutrino experiment can achieve, but, as pointed out in [9], a single detection of a neutrino with an energy higher than 10 PeV would give evidence of a flux beyond what is firmly established. Our results indicate that, although extremely challenging, it is worth trying to detect ultrahigh-energy neutrinos, interacting in the Pico de Orizaba with the HAWC observatory.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

#### Acknowledgments

This work was supported by the Programme UNAM-DGAPA-PAPIIT IA102715, Consejo Nacional de Ciencia y Tecnología (CONACyT), Mexico Grant 254964, and Coordinación de la Investigación Científica UNAM.