Search for Ultra-High-Energy Neutrinos with the Telescope Array Surface Detector

TL;DR

This paper reports an upper limit on ultra-high-energy neutrino flux using nine years of Telescope Array surface detector data, employing a multivariate Boosted Decision Trees method to distinguish neutrino-induced showers.

Contribution

It introduces a novel application of BDTs with 16 observables to identify ultra-high-energy neutrinos in surface detector data.

Findings

Established an upper flux limit for neutrinos above 10^18 eV.

Demonstrated the effectiveness of multivariate analysis in neutrino detection.

Provided constraints for models of ultra-high-energy neutrino sources.

Abstract

We present an upper limit on the flux of ultra-high-energy down-going neutrinos for $E > 10^{18}\ \mbox{eV}$ derived with the nine years of data collected by the Telescope Array surface detector (05-11-2008 -- 05-10-2017). The method is based on the multivariate analysis technique, so-called Boosted Decision Trees (BDT). Proton-neutrino classifier is built upon 16 observables related to both the properties of the shower front and the lateral distribution function.

Figures

Equations

$$t\left(r\right)=t_0+t_{\text{mbox}plane}+a\times {\left(1+r/R_L\right)}^{1.5}LDF{\left(r\right)}^{-0.5},$$

$$S\left(r\right)=S_{800}\times LDF\left(r\right),$$

$$f\left(r\right)={\left(\frac{r}{R_m}\right)}^{-1.2}{\left(1+\frac{r}{R_m}\right)}^{-(\eta -1.2)}{\left(1+\frac{r^2}{R_1^2}\right)}^{-0.6},$$

$$t_{\text{mbox}plane}^i=\frac{1}{c}\vec{n}\left({\vec{R}}_i-{\vec{R}}_{\text{mbox}core}\right),$$

$$R_m=90.0\text{leavevmode}\text{mbox}m,\text{leavevmode}{R}_1=1000\text{leavevmode}\text{mbox}m,\text{leavevmode}{R}_L=30\text{leavevmode}\text{mbox}m,$$

$$\eta =3.97-1.79\left(\sec \left(\theta \right)-1\right),$$

$$r=\sqrt{({\vec{R}}_i-{\vec{R}}_{core}{)}^2-(\vec{n}({\vec{R}}_i-{\vec{R}}_{core}){)}^2},$$

$$AoP\left(r\right)=\alpha -\beta \left(r/r_0-1.0\right),$$

$$S_{b}=\sum_{i=1}^{N}\biggl{[}S_{i}\times{\left(\frac{r_{i}}{r_{0}}\right)}^{b}\biggr{]},$$

$${\xi }_{cut}\left(\theta \right)={\xi }_0+{\xi }_1\ast \theta +{\xi }_2\ast {\theta }^2.$$

$$f_{merit}\left({\xi }_0,{\xi }_1,{\xi }_2,\theta \right)=\frac{⟨n_p{⟩}_{90\%\text{mbox}C.L.}}{n_{\nu }},$$

$$A_{eff}^{\nu }=A_{geom}^{MC}\times \frac{N_{pass}}{N_{thrown}}\times N_{flavor}.$$

$$A_{eff}^{\nu }=4.2\times 10^{-3}{\text{mbox}km}^2\text{mbox}sr\text{mbox}yr.$$

$$F_{\nu }=\frac{{\overline{n}}_{\nu }}{A_{eff}^{\nu }},$$

$$E{F}_{\nu }<1.85\times 10^{-6}\text{mbox}GeV\text{mbox}{cm}^{-2}\text{mbox}{s}^{-1}\text{mbox}{sr}^{-1}\left(90\%\text{mbox}C.L.\right).$$

$${\hat{A}}_{eff}^{\nu }=1.2\times 10^{-2}{\text{mbox}km}^2\text{mbox}sr\text{mbox}yr.$$

$$E{\hat{F}}_{\nu }<4.2\times 10^{-6}\text{mbox}GeV\text{mbox}{cm}^{-2}\text{mbox}{s}^{-1}\text{mbox}{sr}^{-1}\left(90\%\text{mbox}C.L.\right).$$

Full text

I Introduction

The Telescope Array (TA) experiment is the largest ultra-high-energy (UHE) cosmic-ray experiment in the Northern Hemisphere, located near Delta, Utah, USA Tokuno . TA is designed to register the extensive air showers (EAS), cascades of secondary particles produced in the interactions of cosmic rays with energies greater than $10^{18}\ \mbox{eV}$ with the Earth’s atmosphere. In the Telescope Array, air showers are registered in two ways: particle density and the shower timing on the ground are measured with the surface detector (SD) array TASD , while the fluorescence light from gas molecules in the atmosphere excited and ionized by the passage of EAS particles is detected with 38 fluorescence telescopes grouped into three fluorescence detector stations – Middle Drum, Black Rock Messa and Long Ridge Tokuno2 . The simultaneous use of both surface detectors and fluorescence telescopes is known as the hybrid technique.

The surface detector facility is an array of 507 plastic scintillator stations arranged on a square grid with 1.2 km spacing covering an approximate area of 700 ${\mbox{km}}^{2}$. Each detector is composed of two layers of 1.2 cm thick extruded scintillator of the $3\ {\mbox{m}}^{2}$ effective area.

Ultra-high-energy cosmic rays pose a number of long-standing questions in astrophysics: neither their origin, nor the acceleration mechanisms responsible for the observed UHECR flux have been identified so far. Recent observation of dipole anisotropy by the Pierre Auger Observatory Aab:2017tyv indicates that highest-energy cosmic rays are born in the extragalactic sources. Because of their properties, cosmic-ray neutrinos are thought to be one of the efficient tools to study UHECR models Aloisio:2015ega ; vanVliet:2017obm ; Heinze:2015hhp ; Giacinti:2015pya . Neutrinos and photons are produced in hadron-hadron and photon-hadron interactions through the production and decay of pions and kaons. In contrast to photons, a negligible fraction of neutrinos is absorbed during their propagation to Earth due to the small cross-sections of its interactions with other particles. Neutrinos are not deflected by magnetic field which means they point directly to their sources and may give a hint on the type of UHECR sources and on the mechanisms responsible for the acceleration of a parent particle Waxman:1998yy ; Waxman:1997ti ; Engel:2001hd ; Murase:2010gj ; Aartsen:2016ngq . Moreover, if the primary particles are not purely protons, the neutrino flux will be suppressed considerably, so they may also be a probe of a mass composition of the UHECR Aloisio:2017waw ; vanVliet:2019nse .

UHE neutrinos may be born in the Universe in three different types of sources and processes:

Astrophysical neutrinos are born in the hadron interactions of the UHECR with radiation or matter near its astrophysical sources. Most promising among the sources are active galactic nuclei (AGN) and its sub-class, blazars, where neutrinos are born in the proton interactions with the AGN photon fields Stecker:1991vm ; galaxy clusters Berezinsky:1996wx and starburst galaxies Thompson:2006np , whose intergalactic and interstellar medium serves as a target for UHECR; supernovae Budnik:2007yh and hypernovae Wang:2007ya , suggested for a long time as the cosmic-ray accelerators; and gamma-ray bursts (GRBs) Waxman:1997ti . 2. 2.

Primary particles and nuclei may interact with the cosmic microwave background (CMB) and the extragalactic background light (EBL) during their propagation in the Universe. Primary protons undergo photo-pion production ($p+\gamma\rightarrow p/n+\pi^{0}/\pi^{+}$) process, while the relevant process for nuclei is the photodisintegration ($\left(A,Z\right)+\gamma\rightarrow\left(A^{\prime},Z^{\prime}\right)+\left(Z-Z^{\prime}\right)p+\left(A-A^{\prime}+Z-Z^{\prime}\right)n$). The photo-hadronic interactions at the highest energies should result in the flux of the so-called cosmogenic neutrinos Berezinsky:1970xj ; Stecker:1973sy ; Hill:1983xs . 3. 3.

Lastly, UHE neutrinos may be born in top-down models in decays of massive objects, such as heavy dark matter particles and topological defects Berezinsky:1991aa ; Bhattacharjee:1991zm ; AlvarezMuniz:2000es ; Gelmini:1999ds ; Kusenko:2000fk in the processes such as $D\rightarrow\nu+all$, or a possible rare decay $D\rightarrow 3\nu$.

Multiple methods have been used to search for UHE neutrinos Baret:2011zz ; Anchordoqui:2009nf . Neutrinos may interact with the Earth’s atmosphere through charged ($\nu+N\rightarrow\mbox{lepton}+X$) and neutral ($\nu+N\rightarrow\nu+X$) current interactons (CC and NC from now on). Electrons and tau-leptons produced in NC interactions as well as final-state products $X$ of hadronic interactions generate extensive air showers, which, in turn, may be detected with ground arrays, fluorescence telescopes and radio antennas. Such neutrino events are usually called down-going.

However, neutrino-air interaction cross-sections are much smaller than the ones for protons or other nuclei Kusenko:2001gj thus the probability to invoke an EAS is also suppressed. To overcome this obstacle, it was suggested in Beresinsky:1969qj ; Berezinsky:1975zz , that air showers initiated by the UHE neutrinos may be observed at large zenith angles, which increases the slant depth neutrinos travel in the atmosphere. The probability of neutrino-air interaction is constant at any point along their trajectories, and this allows to distinguish neutrino-induced showers as the ones which develop deep in the atmosphere, unlike the protons and nuclei showers which are initiated in the upper layers.

Another possibility for neutrino to invoke an EAS is the CC interactions with the minerals in the Earth’s crust – so-called Earth-skimming events. They usually occur close to the exit point below the surface of the Earth and the EAS develops upwards in the atmosphere. Whilst muon and electron neutrinos produce well-contained tracks and cascades after one CC-interaction, ${\nu}_{\tau}$ generates a $\tau$-lepton, and a characteristic interaction length of such process is about $20\ \mbox{km}$. The decay length of a $\tau$-lepton is then sufficient enough for it to escape the Earth’s crust and generate an EAS Bottai:2002nn ; Feng:2001ue .

Thirdly, one can observe the radio-emission from neutrino passing through dense matter, such as ice or lunar regolith, caused by the Askaryan effect Askaryan . In the Askaryan effect, the coherent radiowave Cherenkov radiation is initiated by the particles moving through the medium with speeds close to the speed of light.

Previously, cosmic-ray neutrino searches have been performed by various experiments. At the ultra-high energies, neutrino studies were performed by the Fly’s Eye Baltrusaitis:1985mt , the HiRes Abbasi:2008hr , the Westerbork Synthesis Radio Telescope Scholten:2009ad , the RESUN project Jaeger:2009zb , the LUNASKA James:2009sf ; Bray:2015lda and RICE Kravchenko:2011im experiments, the Pierre Auger Observatory Aab:2015kma ; Aab:2019auo , the ARIANNA array Barwick:2014pca and the ANITA Gorham:2008yk , ANITAII Gorham:2010 and ANITAIII balloon payloads Allison:2018cxu . No neutrinos with energies more than $10^{16}\ \mbox{eV}$ have been detected so far and only the upper limits on the differential flux were derived.

At the TeV to PeV energies, results were obtained with the Askaryan Radio Array (ARA) Allison:2015eky ; Allison:2019xtn , the ANTARES telescope Albert:2017oba and the IceCube observatory Albert:2018vxw . An astrophysical neutrino flux above the atmospheric neutrino background was discovered by IceCube Aartsen:2013jdh in the three years of its data. The latest measured astrophysical neutrino flux based on the IceCube six-year data is comprised of 82 events Aartsen:2017mau , including two neutrinos with energies above 1 PeV, and a 2 PeV event Aartsen:2014gkd .

The present Paper is dedicated to the ultra-high-energy down-going neutrino search using solely the data from the Telescope Array Surface Detector array TASD . The data from the nine years of the Telescope Array surface detector operation from May 11, 2008 to May 10, 2017 are used in the present analysis. Fluorescence detectors operate only on clear moonless nights, which results in approximately 10 % duty cycle, while the SD array duty cycle reaches 95 % which make the use of its data much more beneficial in terms of acquired statistics.

While there are no observables, which can be obtained from the SD data and are as sensitive to the composition of primary particles, as the depth of the shower maximum, $X_{\mbox{max}}$, an alternative approach was proposed. The data analysis employs the multivariate analysis techniques to obtain the proton-neutrino Boosted Decision Trees (BDTs) classifier Breiman ; Schapire built upon 16 observables. For each event, a classifier returns a single value $\xi$, which is then available to one-dimensional anaylsis and allows one to distinguish between different hypotheses.

The Paper is organized as follows: in the Section II data and Monte-Carlo sets are described. In the Section III multivariate analysis method is addressed together with its implementation to the neutrino search. Finally, results and discussions are provided in Section IV.

II Data set and simulations

II.1 Data set

The data from the nine years of the Telescope Array surface detector operation from May 11, 2008 to May 10, 2017 are used in the present Paper. Each event is comprised of the time-dependent signals (waveforms) from both upper and lower layers of each triggered detector. The waveforms are recorded by the 12-bit flash analog-to-digital converters (FADC) with the 50 MHz sampling rate and are converted to Minimum Ionizing Particles (MIP) TASD at the calibration stage. The station is marked as triggered if the signal exceeds $0.3\ \mbox{MIP}$.

II.2 Event reconstruction and cuts

Each surface detector array event is comprised of a set of observables $t_{i}\left(r_{i}\right)$ and $S_{i}\left(r_{i}\right)$, where $t_{i}\left(r_{i}\right)$ is the time of arrival of an EAS to the $i$-th detector of the event and $S_{i}\left(r_{i}\right)$ is the signal from the station. These observable quantities are used to perform a joint fit of the geometry and lateral distribution function (LDF). The shower front is approximated with empirical functions proposed by Linsley Linsley and later modified in the AGASA experiment Teshima . The $S_{800}$ parameter, the signal value at $800\ \mbox{m}$ from the shower core used as an energy estimator, is then obtained using pulse heights in the counters together with the event geometry information Takeda .

For comparison, a number of Monte-Carlo events are generated with the CORSIKA package Heck and GEANT4-based detector simulator Agostinelli . MC events are also reconstructed with the same code and procedure for the data.

The Linsley front curvature parameter is determined through the following fit to the shower front using the LDF with 7 free parameters: $x_{\mbox{core}}$, $y_{\mbox{core}}$, $\theta$, $\phi$, $S_{800}$, $t_{0}$, $a$ TAgammalim :

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\vec{x}_{\mbox{core}}$ and $\vec{y}_{\mbox{core}}$ are the locations of the shower core, $\vec{x}_{i}$ and $\vec{y}_{i}$ are the locations of each station of an event, obtained from the pre-defined coordinate system of the array centered at the Central Laser Facility (CLF) Takahashi:2011zzd , $t_{plane}$ is the arrival timing of the shower plane at the distance $r$, $\vec{n}$ – unit vector towards the direction of arrival of a primary particle, $c$ is a speed of light and $a$ is the Linsley front curvature parameter.

The following cuts are applied:

the event includes 5 or more triggered stations; 2. 2.

zenith angle $\theta\in[45^{\circ};90^{\circ}]$; 3. 3.

reconstructed core position inside the array with the distance of at least $1200\ \mbox{m}$ from the edge of the array; 4. 4.

Effective ${\chi}^{2}/\mbox{d.o.f.}$ doesn’t exceed 5 for the joint geometry and LDF fit.

No energy cut is applied, and after the cuts, the SD data set contains 197250 events.

II.3 Simulations

For the Monte-Carlo simulations, the CORSIKA software package Heck is used along with the QGSJETII-03 model for high-energy hadronic interactions Ostapchenko , FLUKA FLUKA ; FLUKA2 for low energy hadronic interaction and EGS4 EGS4 for electromagnetic processes. Interactions of primary neutrinos are handled with HERWIG Corcella:2000bw integrated with CORSIKA. While a newer version of the QGSJET model, QGSJETII-04 QGSJETII:04 , is also widely incorporated in the related analyses, the previous version is found to be adequate to describe the data Abbasi:2018wlq ; and an upper limit on the neutrino flux is expected to be conservative in respect to the hadronic interaction models as discussed in Section IV.

The neutrino MC set is comprised of 3000 CORSIKA showers generated in the energy range $3\times 10^{17}-3\times 10^{20}\ \mbox{eV}$. Due to the low neutrino-interaction cross-section a primary neutrino thrown with CORSIKA will most likely travel through the Earth’s atmosphere without interaction. To overcome this obstacle the height of the primary interaction is simulated with a separate routine and fixed manually within CORSIKA.

The procedure consists of the following steps: at first, a uniform neutrino flavour flux ${\nu}_{e}:\bar{{\nu}}_{e}:{\nu}_{\mu}:\bar{{\nu}}_{\mu}:{\nu}_{\tau}:\bar{{\nu}}_{\tau}=1:1:1:1:1:1$ is assumed and a randomly chosen flavour is assigned for the event. Then, zenith angle $\theta\in[0;90^{\circ}]$ is thrown assuming an isotropic flux. For a given zenith angle, the slant depth of the atmosphere $T_{atm}$ is calculated.

Then, we continue with an event if the interaction happens in the atmosphere. The probability of the latter is given by $T_{\mbox{atm}}/T_{\mbox{int}}$, where $T_{\mbox{int}}=\left(M/N_{A}\right)/\sigma_{CC+NC}$ is the interaction slant depth in low interaction probability approximation $T_{atm}/T_{int}\ll 1$, where $M$ is the average molar mass of air, $N_{A}$ is the Avogadro number. Neutrino-nucleus interaction cross-sections $\sigma_{CC+NC}$ are taken from CooperSarkar:2007cv .

Finally, the interaction slant depth is thrown in the range $[0;T_{atm}]$ and the derived value is fixed within CORSIKA input file. Simulated neutrino showers were used to throw 80 millions of events on the detector by randomly distributing the core location and azimuthal angle.

The highly-inclined proton MC set consists of 2400 proton showers with $45^{\circ}<\theta<90^{\circ}$ in the energy range $3\times 10^{17}-3\times 10^{20}\ \mbox{eV}$ used to throw 100 millions of events of each type on the SD array.

For both neutrino and highly-inclined proton MC sets, the usual CORSIKA treatment of atmosphere as a flat disk becomes inefficient at high zenith angles, and one needs to take the curvature of the atmosphere into account. This is done within CORSIKA with a special CURVED atmosphere option Heck .

CPU-time-saving thinning procedure is utilized for both of the MC sets Hillas . Instead of following each of millions of particles born in an EAS, the method follows in detail only those of the particles which have an energy above a certain threshold, which is specified by a fraction of the primary energy $\epsilon_{th}$. Below this limit, only one particle out of the secondaries produced in a certain interaction is randomly selected and assigned with a weight to ensure energy conservation. In the present work, the thinning level of $\epsilon_{th}=10^{-6}$ is used.

Statistical properties of a shower are restored with a dethinning procedure Stokes . The detector response is simulated by the GEANT4 package Agostinelli . Real-time array status and detector calibration information for nine years of observations are used for each simulated event TAdataMC .

II.4 Composition-sensitive observables

During the reconstruction procedure, the following composition-sensitive observables are obtained for each data and MC event:

Linsley front curvature parameter, as described in section II.2. 2. 2–3.

Area-over-peak (AoP) of the signal at 1200 m and AoP slope parameter Abraham :

• $AoP\left(r\right)$ is fitted with a linear fit:

[TABLE]

where $r_{0}=1200$ m, $\alpha$ is $AoP\left(r\right)$ value at 1200 m and $\beta$ is its slope parameter.

Number of detectors hit. 2. 5.

Number of detectors excluded from the fit of the shower front by the reconstruction procedure AbuZayyad3 . 3. 6.

Effectuve $\chi^{2}/\mbox{d.o.f.}$ 4. 7–8.

$S_{b}$ parameter for $b=3$ and $b=4.5$ Ros . The definition of the parameter is the following:

[TABLE]

where $S_{i}$ is the signal of i-th detector, $r_{i}$ is the distance from the shower core to this station in meters and $r_{0}=1200\ \mbox{m}$ – reference distance. The value $b=3$ and $b=4.5$ are used as they provide the best separation. 5. 9.

The sum of the signals of all the detectors of the event. 6. 10.

Asymmetry of the signal at the upper and lower layers of detectors. 7. 11.

Total number of peaks within all FADC (flash analog-to-digital converter) traces. Number of peaks is averaged over both upper and lower layers of all stations participating in the event. FADC peak is defined as a time bin with a signal above 0.2 Vertical equivalent muons (VEM) higher than signal in the 3 preceding and 3 following time bins. 8. 12.

Number of peaks for the detector with the largest signal. 9. 13.

Number of peaks present in the upper layer and not in the lower. 10. 14.

Number of peaks present in the lower layer and not in the upper.

Together with the zenith angle and $S_{800}$, the aforementioned variables comprise a set of 16 composition-sensitive observales used to build a multivariate classifier. They are discussed in more detail in Abbasi:2018wlq .

III Method

III.1 Boosted Decision Trees

The Boosted Decision Tree (BDT) concept belongs to a number of so-called multivariate analysis techniques which allow to efficiently treat the data described by vector-type variables. BDTs are used in a variety of problems, where it is necessary to define, whether a particular event belongs to the signal or to the background. Boosted Decision Trees were previously succesfully used in a number of problems related to cosmic-ray data analysys Abbasi:2018wlq ; Abbasi:2018ywn ; Abbasi:2019mwi .

In general, a single classifier, a tree is built as follows:

For each variable a splitting value with best separation is found. This value divides the full range of the values of the variable into two ranges, which are called branches. It will be mostly signal in one branch, and mostly background in another; 2. 2.

The algorithm is repeated recursively on each branch. It may use a new variable or reuse the same one; 3. 3.

The decision tree will iterate until the stopping criterion is reached (for example, number of events in a branch). The terminal node is called a leaf.

A single classifier may not be efficient enough to provide a good separation between signal and background events, in this case it is called weak. Instead, one may build a number of weak classifiers – a forest to create a strong one. This is the idea which is the basis of the concept of boosting. The adaptive boosting (AdaBoost) algorithm is employed in the present work Schapire ; Freund with the number of trees $N_{\rm Trees}=1000$. In AdaBoost, a weak classifier is run multiple times on the training data, and each event is weighted by how incorrectly it was classified. An improved tree with reweighted events may now be built, and as a result, averaging over all trees allows the creation of a better classifier.

The method finally results in a single value $\xi$ for each event. $\xi$ resides in the range $\xi\in[-1;1]$, where $\xi=1$ corresponds to a pure signal event , $\xi=-1$ – pure background event.

In this Paper, BDT method available as a part of the ROOT root Toolkit for Multivariate Data Analysis (TMVA) package Hocker is used. The neutrino and highly-inclined proton Monte-Carlo sets are split into three and two parts with equal statistics, correspondingly. The first part of both sets is used to build and train the BDT classifier based on 16 variables listed in II.4. Highly-inclined-proton-induced MC showers are used as a background and neutrino-induced ones as signal events. The second part is used for the cut optimization. The third part of the neutrino MC set is used for the exposure calculation. The classifier is applied to the data set as well to the second and third parts of the Monte-Carlo sets.

The zenith angle distribution for MC events plotted against the $\xi$ parameter is shown in Figure 1, where neutrino MC is shown with green pluses, while highly-inclined proton MC is shown with red crosses. The $\left(\xi,\ \theta\right)$ points are used in the cut optimization.

III.2 Cut optimization

The $\xi$ parameter derived in the Section III.1 is then used for the neutrino candidate search. For this purpose, one needs to derive a cut on $\xi$ which depends on the zenith angle. The cut is derived by using the highly-inclined proton MC set as the null-hypothesis, a set without any neutrino events.

In this case, a neutrino candidate should have a $\xi$ parameter larger, than $\xi_{cut}\left(\theta\right)$, where the latter is defined with a quadratic polynomial function:

[TABLE]

$\xi_{cut}\left(\theta\right)$ is derived with the use of highly-inclined proton and neutrino MC sets by minimizing the merit function, which corresponds to the mean expected value of the upper limit on neutrino flux. The merit function is calculated in the following way: some initial values of $\xi_{0}$, $\xi_{1}$ and $\xi_{2}$ are assumed as a starting point for the optimization. Applying the ‘‘starting’’ cut $\xi_{cut}^{0}$ to the highly-inclined proton and neutrino MC sets one obtains the number $n_{p}$ and $n_{\nu}$ of events, correspondingly, which pass the applied cut.

$n_{\nu}$ is proportional to the neutrino exposure of the SD. The number of neutrino candidates is derived from $n_{p}$. Among the events from the ‘‘null hypothesis’’ set which pass the cut, there are both proton events which pass the cut due to the incomplete separation between proton and neutrino events by the $\xi$ variable as well as actual neutrino events. Thus, we define the number of neutrino candidates as the random Poisson variable with the number of observed events being the $n_{p}$, normalized to the statistics of the experiment.

Since $n_{p}$ is always a small value, we apply the Feldman-Cousins method FC to estimate the number of neutrino candidates as the 90% C.L. upper limit on the mean value $\langle n_{p}\rangle$ of the Poisson variable with the number of observed events $n_{p}$ and the expected number of background events $b=0$.

Finally, the merit function is estimated as:

[TABLE]

and it is minimized with the use of simplex algorithm Murty .

The obtained optimized $\xi_{cut}=0.3161+0.0612\times\theta+0.0076\times{\theta}^{2}$ is shown in Figure 1, where zenith angle is assumed to be in radians.

III.3 Exposure calculation

The geometrical exposure for the 9 year SD observational period for the $0^{\circ}<\theta<90^{\circ}$ is $A_{geom}^{MC}=55500\ {\mbox{km}}^{2}\ \mbox{sr}\ \mbox{yr}$.

The exposure is then obtained as the ratio of number of MC neutrino events which have passed all the quality cuts and the $\xi$ cut to the number of neutrinos thrown in the atmosphere, multiplied by the geometrical exposure and the number of neutrino flavors:

[TABLE]

It corresponds to the number of neutrino MC events which pass all the cuts $N_{pass}=8460$ (21.3 % of the number of neutrino set events which have passed the cuts). The number of events thrown in the atmosphere is $N_{thrown}=2.81\times 10^{11}$, which finally leads to the effective exposure for the down-going neutrino:

[TABLE]

IV Results and discussion

Figure 2 shows $\xi$ parameter distribution histogram in the energy range $E>10^{18.0}\ \mbox{eV}$, where proton MC is shown with red line, neutrino MC is shown with green line, while black dots represent the data. One may see that the highly-inclined proton MC set does not follow the data distribution in detail due to the heavier composition as obtained in Abbasi:2018wlq as well as due to the possible systematic effects related to hadronic interaction models, including the “muon excess” problem Abbasi:2018fkz , which is expected to manifest itself more on the higher zenith angles.

The zenith angle-$\xi$ scatter plot for data is shown in Figure 3 with the obtained $\xi_{cut}$ function.

There are no neutrino candidates in the data set and hence using FC , we derive the upper limit on the number of neutrino events of all flavors in the data: $\overline{n}_{\nu}=2.44$ at 90% C.L.

By definition, the integral neutrino flux depends on the number of neutrino events and on the effective neutrino exposure:

[TABLE]

thus we arrive at the following upper limit on the single flavor diffuse neutrino flux for $E>10^{18}\ \mbox{eV}$:

[TABLE]

The obtained limit in comparison with the Pierre Auger Observatory Aab:2015kma , ANITA-II ANITAII , RICE Kravchenko:2011im , IceCube Aartsen:2013dsm , ARA Allison:2019xtn and HiRes Abbasi:2008hr results is shown in Figure 4 as well as the the Waxman-Bahcall limit Waxman:1998yy .

One may be also interested in the results of the neutrino search without the cut optimization procedure explained in III.2. For this purpose, following the line of Homola , we employ the cut based on the median of the $\xi$ distribution for neutrino.

From the distribution of the $\xi$ parameter for neutrino set, see Fig. 2, one obtains the median value of $\xi_{med}=0.205$, which corresponds to 50% of MC neutrino passing the $\xi_{med}$ cut ($N_{pass}=19856$).

The effective exposure for the down-going neutrino in this case:

[TABLE]

One obtains $35$ neutrino candidates in the data, which corresponds to the upper limit on neutrino events $\hat{\overline{n}}_{\nu}=46.4$ at 90% C.L. It leads to the upper limit on the single flavor diffuse neutrino flux for $E>10^{18}\ \mbox{eV}$, shown in Fig. 4:

[TABLE]

One should note note that given any fixed cut on $\xi$ one arrives at valid neutrino flux limit, while as expected, the optimization procedure results in the stronger limits.

Let us finally discuss the effect of the hadronic interaction model systematics and the composition uncertainty on the result of the Paper. Both these effects infer the analysis through the Monte-Carlo set used as a background, i.e. the highly-inclined proton MC set. This set is used for the construction of the classifier, including the BDT algorithm and the optimization of the cut. An optimal sensitivity is by construction achieved in the case of a perfect agreement of the data and Monte-Carlo. Incomplete desciption of the observed events by a Monte Carlo set leads either to reduced exposure or to an increased number of false candidates compared to an optimal classifier. Therefore the neutrino flux limit derived is a conseravtive estimate which is close to optimal limit given the reasonable agreement of data and Monte-Carlo demonstrated in the present Paper. The sensitivity of the method can be further improved by producing hadronic Monte-Carlo which better reproduces the bulk of data events.

Acknowledgment

The Telescope Array experiment is supported by the Japan Society for the Promotion of Science (JSPS) through Grants-in-Aid for Priority Area 431, for Specially Promoted Research JP21000002, for Scientific Research (S) JP19104006, for Specially Promoted Research JP15H05693, for Scientific Research (S) JP15H05741, for Science Research (A) JP18H03705 and for Young Scientists (A) JPH26707011; by the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo; by the U.S. National Science Foundation awards PHY-0601915, PHY-1404495, PHY-1404502, and PHY-1607727; by the National Research Foundation of Korea (2016R1A2B4014967, 2016R1A5A1013277, 2017K1A4A3015188, 2017R1A2A1A05071429); IISN project No. 4.4502.13, and Belgian Science Policy under IUAP VII/37 (ULB). The development and application of the multivariate analysis method is supported by the Russian Science Foundation grant No. 17-72-20291 (INR). The foundations of Dr. Ezekiel R. and Edna Wattis Dumke, Willard L. Eccles, and George S. and Dolores Doré Eccles all helped with generous donations. The State of Utah supported the project through its Economic Development Board, and the University of Utah through the Office of the Vice President for Research. The experimental site became available through the cooperation of the Utah School and Institutional Trust Lands Administration (SITLA), U.S. Bureau of Land Management (BLM), and the U.S. Air Force. We appreciate the assistance of the State of Utah and Fillmore offices of the BLM in crafting the Plan of Development for the site. Patrick Shea assisted the collaboration with valuable advice on a variety of topics. The people and the officials of Millard County, Utah have been a source of steadfast and warm support for our work which we greatly appreciate. We are indebted to the Millard County Road Department for their efforts to maintain and clear the roads which get us to our sites. We gratefully acknowledge the contribution from the technical staffs of our home institutions. An allocation of computer time from the Center for High Performance Computing at the University of Utah is gratefully acknowledged. The cluster of the Theoretical Division of INR RAS was used for the numerical part of the work. This work was partially supported by the Collaborative research program of the Institute for Cosmic Ray Research (ICRR), the University of Tokyo.

Список литературы

• (1)

H. Tokuno et al. [Telescope Array Collaboration], J. Phys. Conf. Ser.  293, 012035 (2011).

• (2)

T. Abu-Zayyad et al. [Telescope Array Collaboration], Nucl. Instrum. Meth. A 689, 87 (2013) [arXiv:1201.4964 [astro-ph.IM]].

• (3)

H. Tokuno et al. [Telescope Array Collaboration], Nucl. Instrum. Meth. A 676, 54 (2012) [arXiv:1201.0002 [astro-ph.IM]].

• (4)

A. Aab et al. [Pierre Auger Collaboration], Science 357, no. 6537, 1266 (2017) [arXiv:1709.07321 [astro-ph.HE]].

• (5)

R. Aloisio, D. Boncioli, A. di Matteo, A. F. Grillo, S. Petrera and F. Salamida, JCAP 1510, no. 10, 006 (2015) [arXiv:1505.04020 [astro-ph.HE]].

• (6)

A. van Vliet, J. R. Hörandel and R. Alves Batista, arXiv:1707.04511 [astro-ph.HE].

• (7)

J. Heinze, D. Boncioli, M. Bustamante and W. Winter, Astrophys. J.  825, no. 2, 122 (2016) [arXiv:1512.05988 [astro-ph.HE]].

• (8)

G. Giacinti, M. Kachelriess, O. Kalashev, A. Neronov and D. V. Semikoz, Phys. Rev. D 92, no. 8, 083016 (2015) [arXiv:1507.07534 [astro-ph.HE]].

• (9)

E. Waxman and J. N. Bahcall, Phys. Rev. D 59, 023002 (1999) [arXiv:hep-ph/9807282].

• (10)

E. Waxman and J. N. Bahcall, Phys. Rev. Lett.  78, 2292 (1997) [arXiv:astro-ph/9701231].

• (11)

R. Engel, D. Seckel and T. Stanev, Phys. Rev. D 64, 093010 (2001) [arXiv:astro-ph/0101216].

• (12)

K. Murase and J. F. Beacom, Phys. Rev. D 81, 123001 (2010) [arXiv:1003.4959 [astro-ph.HE]].

• (13)

M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. Lett.  117, no. 24, 241101 (2016) [arXiv:1607.05886 [astro-ph.HE]].

• (14)

R. Aloisio, Nuovo Cim. C 40, no. 3, 142 (2017) [arXiv:1707.06188 [astro-ph.HE]].

• (15)

A. van Vliet, R. Alves Batista and J. R. Horandel, arXiv:1901.01899 [astro-ph.HE].

• (16)

F. W. Stecker, C. Done, M. H. Salamon and P. Sommers, Phys. Rev. Lett.  66, 2697 (1991).

• (17)

V. S. Berezinsky, P. Blasi and V. S. Ptuskin, Astrophys. J.  487, 529 (1997) [arXiv:astro-ph/9609048].

• (18)

T. A. Thompson, E. Quataert, E. Waxman and A. Loeb, [arXiv: astro-ph/0608699].

• (19)

R. Budnik, B. Katz, A. MacFadyen and E. Waxman, Astrophys. J.  673, 928 (2008) [arXiv:0705.0041 [astro-ph]].

• (20)

X. Y. Wang, S. Razzaque, P. Meszaros and Z. G. Dai, Phys. Rev. D 76, 083009 (2007) [arXiv:0705.0027 [astro-ph]].

• (21)

V. S. Berezinsky and G. T. Zatsepin, Yad. Fiz.  11, 200 (1970).

• (22)

F. W. Stecker, Astrophys. Space Sci.  20, 47 (1973).

• (23)

C. T. Hill and D. N. Schramm, Phys. Lett. B 131, 247 (1983).

• (24)

V. S. Berezinsky, Nucl. Phys. B 380, 478 (1992).

• (25)

P. Bhattacharjee, C. T. Hill and D. N. Schramm, Phys. Rev. Lett.  69, 567 (1992).

• (26)

J. Alvarez-Muniz and F. Halzen, Phys. Rev. D 63, 037302 (2001) [arXiv:astro-ph/0007329].

• (27)

G. Gelmini and A. Kusenko, Phys. Rev. Lett.  84, 1378 (2000) [arXiv:hep-ph/9908276].

• (28)

A. Kusenko and M. Postma, Phys. Rev. Lett.  86, 1430 (2001) [arXiv:hep-ph/0007246].

• (29)

B. Baret and V. Van Elewyck, Rept. Prog. Phys.  74 (2011) 046902.

• (30)

L. A. Anchordoqui and T. Montaruli, Ann. Rev. Nucl. Part. Sci.  60, 129 (2010) [arXiv:0912.1035 [astro-ph.HE]].

• (31)

A. Kusenko and T. J. Weiler, Phys. Rev. Lett.  88, 161101 (2002) [arXiv:hep-ph/0106071].

• (32)

V. S. Berezinsky and G. T. Zatsepin, Phys. Lett.  28B, 423 (1969).

• (33)

V. S. Berezinsky and A. Y. Smirnov, Astrophys. Space Sci.  32, 461 (1975).

• (34)

S. Bottai and S. Giurgola, Astropart. Phys.  18, 539 (2003) [arXiv:astro-ph/0205325].

• (35)

J. L. Feng, P. Fisher, F. Wilczek and T. M. Yu, Phys. Rev. Lett.  88, 161102 (2002) [arXiv:hep-ph/0105067].

• (36)

G. Askaryan, Sov.Phys. JETP 14, 441 (1962).

• (37)

R. M. Baltrusaitis, G. L. Cassiday, J. W. Elbert, P. R. Gerhardy, E. C. Loh, Y. Mizumoto, P. Sokolsky and D. Steck, Phys. Rev. D 31, 2192 (1985).

• (38)

R. U. Abbasi et al., Astrophys. J.  684, 790 (2008) [arXiv:0803.0554 [astro-ph]].

• (39)

O. Scholten et al., Phys. Rev. Lett.  103, 191301 (2009) [arXiv:0910.4745 [astro-ph.HE]].

• (40)

T. R. Jaeger, R. L. Mutel and K. G. Gayley, Astropart. Phys.  34, 293 (2010) [arXiv:0910.5949 [astro-ph.IM]].

• (41)

C. W. James, R. D. Ekers, J. Alvarez-Muniz, J. D. Bray, R. A. McFadden, C. J. Phillips, R. J. Protheroe and P. Roberts, Phys. Rev. D 81, 042003 (2010) [arXiv:0911.3009 [astro-ph.HE]].

• (42)

J. D. Bray et al., Phys. Rev. D 91, no. 6, 063002 (2015) [arXiv:1502.03313 [astro-ph.HE]].

• (43)

I. Kravchenko et al., Phys. Rev. D 85, 062004 (2012) [arXiv:1106.1164 [astro-ph.HE]].

• (44)

A. Aab et al. [Pierre Auger Collaboration], Phys. Rev. D 91, no. 9, 092008 (2015) [arXiv:1504.05397 [astro-ph.HE]].

• (45)

A. Aab et al. [Pierre Auger Collaboration], arXiv:1906.07422 [astro-ph.HE].

• (46)

S. W. Barwick et al. [ARIANNA Collaboration], Astropart. Phys.  70, 12 (2015) [arXiv:1410.7352 [astro-ph.HE]].

• (47)

P. W. Gorham et al. [ANITA Collaboration], Phys. Rev. Lett.  103, 051103 (2009) [arXiv:0812.2715 [astro-ph]].

• (48)

P. W. Gorham et al. [ANITA Collaboration], Phys. Rev. D  82, 022004 (2010) [arXiv:1003.2961 [astro-ph.HE]].

• (49)

P. W. Gorham et al. [ANITA Collaboration], Phys. Rev. D 98, no. 2, 022001 (2018) [arXiv:1803.02719 [astro-ph.HE]].

• (50)

P. Allison et al. [ARA Collaboration], Phys. Rev. D 93, no. 8, 082003 (2016) [arXiv:1507.08991 [astro-ph.HE]].

• (51)

P. Allison et al. [ARA Collaboration], arXiv:1912.00987 [astro-ph.HE].

• (52)

A. Albert et al. [ANTARES Collaboration], Phys. Rev. D 96, no. 6, 062001 (2017) [arXiv:1705.00497 [astro-ph.HE]].

• (53)

A. Albert et al., [arXiv:1808.03531 [astro-ph.HE]].

• (54)

M. G. Aartsen et al. [IceCube Collaboration], Science 342, 1242856 (2013) [arXiv:1311.5238 [astro-ph.HE]].

• (55)

M. G. Aartsen et al. [IceCube Collaboration], [arXiv:1710.01191 [astro-ph.HE]].

• (56)

M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. Lett.  113, 101101 (2014) [arXiv:1405.5303 [astro-ph.HE]].

• (57)

L. Breiman et al., Wadsworth International Group (1984).

• (58)

R.E. Schapire, Mach. Learn. 5 (1990) 197.

• (59)

J. Linsley, L. Scarsi, Phys. Rev. 128 (1962) 2384.

• (60)

M. Teshima et al., J. Phys. G 12, 1097 (1986).

• (61)

M. Takeda et al., Astropart. Phys.  19, 447 (2003) [arXiv:astro-ph/0209422].

• (62)

D. Heck et al., Report FZKA-6019 (1998), Forschungszentrum Karlsruhe.

• (63)

S. Agostinelli et al. [GEANT4 Collaboration], Nucl. Instrum. Meth. A 506, 250 (2003).

• (64)

T. Abu-Zayyad et al. [Telescope Array Collaboration], Phys. Rev. D 88, no. 11, 112005 (2013) [arXiv:1304.5614 [astro-ph.HE]].

• (65)

Y. Takahashi et al. [Telescope Array Collaboration], AIP Conf. Proc.  1367, 157 (2011).

• (66)

S. Ostapchenko, Nucl. Phys. Proc. Suppl.  151, 143 (2006) [arXiv:hep-ph/0412332].

• (67)

T. T. Böhlen et al., Nucl. Data Sheets 120, 211 (2014)

• (68)

A. Ferrari, P. R. Sala, A. Fasso and J. Ranft, CERN-2005-010, SLAC-R-773, INFN-TC-05-11.

• (69)

W. R. Nelson, H. Hirayama, D.W.O. Rogers, SLAC-0265 (permanently updated since 1985).

• (70)

G. Corcella, I. G. Knowles, G. Marchesini, S. Moretti, K. Odagiri, P. Richardson, M. H. Seymour and B. R. Webber, JHEP 0101, 010 (2001) [arXiv:hep-ph/0011363].

• (71)

S. Ostapchenko, Phys. Rev. D 83, 014018 (2011).

• (72)

R. U. Abbasi et al. [Telescope Array Collaboration], arXiv:1808.03680 [arXiv:astro-ph.HE].

• (73)

A. Cooper-Sarkar and S. Sarkar, JHEP 0801, 075 (2008) [arXiv:0710.5303 [hep-ph]].

• (74)

A. M. Hillas, Nucl. Phys. Proc. Suppl.  52B, 29 (1997).

• (75)

B. T. Stokes et al., Astropart. Phys. 35, 759 (2012).

• (76)

T. Abu-Zayyad et al. [Telescope Array Collaboration], [arXiv:1403.0644 [astro-ph.IM]].

• (77)

J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett.  100, 211101 (2008) [arXiv:0712.1909 [astro-ph]].

• (78)

T. Abu-Zayyad et al. [Telescope Array Collaboration], ApJL 768 (2013) L1.

• (79)

G. Ros et al., Astropart. Phys.  35, 140 (2011) [arXiv:1104.3399 [astro-ph.HE]].

• (80)

R. U. Abbasi et al. [Telescope Array Collaboration], Astropart. Phys.  110, 8 (2019), [arXiv:1811.03920 [astro-ph.HE]].

• (81)

R. U. Abbasi et al. [Telescope Array Collaboration], arXiv:1904.00300 [astro-ph.HE].

• (82)

Y. Freund, R.E. Schapire, Proc. ICML (1996) 148.

• (83)

R. Brun and F. Rademakers, Proceedings AIHENP’96 Workshop, Lausanne, Sep. 1996, Nucl. Inst. & Meth. in Phys. Res. A 389 (1997) 81-86. See also http://root.cern.ch/.

• (84)

A. Hocker et al., PoS ACAT (2007) 040 [physics/0703039 [PHYSICS]].

• (85)

G. J. Feldman and R. D. Cousins, Phys. Rev. D 57 (1998) 3873.

• (86)

R. U. Abbasi et al. [Telescope Array Collaboration], Phys. Rev. D 98, no. 2, 022002 (2018) [arXiv:1804.03877 [astro-ph.HE]].

• (87)

K. G. Murty, Wiley (1983).

• (88)

P. W. Gorham et al. [ANITA Collab.], Phys. Rev. D 85, 049901(E) (2012) [arXiv:1003.2961 [astro-ph.HE]].

• (89)

M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. D 88, 112008 (2013) [arXiv:1310.5477 [astro-ph.HE]].

• (90)

A. Aab et al. [Pierre Auger Collaboration], JCAP 1704, 009 (2017) [arXiv:1612.01517 [astro-ph.HE]].