Probing New Physics at Future Tau Neutrino Telescopes

We systematically investigate new physics scenarios that can modify the interactions between neutrinos and matter at upcoming tau neutrino telescopes, which will test neutrino-proton collisions with energies $ \gtrsim 45~{\rm TeV}$, and can provide unique insights to the elusive tau neutrino. At such high energy scales, the impact of parton distribution functions of second and third generations of quarks (usually suppressed) can be comparable to the contribution of first generation with small momentum fraction, hence making tau neutrino telescopes an excellent facility to probe new physics associated with second and third families. Among an inclusive set of particle physics models, we identify new physics scenarios at tree level that can give competitive contributions to the neutrino cross sections while staying within laboratory constraints: charged/neutral Higgs and leptoquarks. Our analysis is close to the actual experimental configurations of the telescopes, and we perform a $\chi^2$-analysis on the energy and angular distributions of the tau events. By numerically solving the propagation equations of neutrino and tau fluxes in matter, we obtain the sensitivities of representative upcoming tau neutrino telescopes, GRAND, POEMMA and Trinity, to the charged Higgs and leptoquark models. While each of the experiments can achieve a sensitivity better than the current collider reaches for certain models, their combination is remarkably complementary in probing the new physics. In particular, the new physics will affect the energy and angular distributions in different ways at those telescopes.


Introduction
The IceCube observatory has made significant progresses in measuring the ultra-high-energy (UHE) astrophysical neutrino flux [1][2][3][4][5]. As an elusive messenger, neutrinos have been utilized along with cosmic rays, gamma rays and gravitational waves to understand the nature of cosmic accelerators [6][7][8][9][10][11][12]. There is a guaranteed flux of UHE neutrinos produced by the scattering of cosmic rays with the cosmic photon background, i.e., the cosmogenic neutrinos [13] with typical energy around EeV ≡ 10 9 GeV, associated with the Greisen-Zatsepin-Kuzmin (GZK) cutoff structure [14][15][16] in the cosmic ray spectrum. After the successful observations of extraterrestrial UHE neutrinos up to PeV energies at IceCube, a campaign of experimental programs has been launched to measure the cosmogenic neutrinos at extreme EeV energy scales. A promising class of such observatories is the tau neutrino telescope, which is sensitive to the ν τ component in the cosmogenic neutrino flux.
The recent anomalies, arising in the measurements of the muon anomalous magnetic moment [79,80] and B-meson decays [81][82][83], have indicated the existence of new physics that has a preferable coupling to the second (c, s and µ) and third (b and τ ) families over the first one. While the second family is still accessible, the third family is difficult to probe in laboratory. In this respect, tau neutrino telescopes are naturally sensitive to the new physics lying in the second and third families. With an EeV incoming neutrino, the COM energy of the neutrino-proton scattering is as high as 45 TeV. This is much higher than what can be achieved in laboratory 1 . At very high energy scales, the parton distribution functions (PDFs) of heavy quarks increase rapidly with small momentum fraction [85], making processes associated with second and third generations of quark partons, which are suppressed by orders of magnitude at LHC, increasingly important for processes at tau neutrino telescopes.
To enhance the effective volume of neutrino interactions, tau neutrino telescopes will be deployed at a high altitude [86][87][88][89][90][91][92][93][94][95][96][97], e.g., located on mountains, balloon-, or satellite-borne. They will monitor the extensive air shower events emerged from the surface of the Earth target. Since the cosmic rays are shielded by the thick Earth medium, if there is such a shower event, most Earth POEMMA Trinity GRAND Figure 1: A cartoon of the detection principles of tau neutrino telescopes. Three representative programs with different geography are shown: the mountain-valley telescope GRAND (in purple), the mountain-top telescope Trinity (in red), and the space-borne project POEMMA (in blue). A tau neutrino enters the Earth or mountain, and scatters with matter to produce a tau lepton (the first bang, invisible). The tau lepton eventually decays in the air (the second bang), generating extensive air showers which induce radio signals for GRAND, or Cherenkov light for Trinity and POEMMA. For GRAND, both mountain-penetrating and Earth-skimming neutrinos are observable, but Trinity and POEMMA (in Limb mode) mainly look for Earth-skimming events near the horizon.
probably it is produced by a tau decaying in the air 2 , which in turn is generated from a tau neutrino interacting with the Earth. The extensive air shower events can then be detected in the form of radio waves, Cherenkov light, or fluorescence [98][99][100][101][102]. A schematic diagram for the working principles of tau neutrino telescopes is given in Fig. 1, where we have chosen three representative tau neutrino telescopes GRAND [103,104], Trinity [105][106][107][108] and POEMMA [109].
Given the numerous upcoming programs which search for neutrinos at EeV energies, we intend to systematically investigate the sensitivities of tau neutrino telescopes to the particle physics models modifying neutrino-matter interactions [110][111][112][113], with an emphasis on the tau sector. The new physics scenarios we consider enter into neutrino scattering processes as in Fig. 2, which directly modifies the interactions between neutrinos and normal matter. Note that in this work, we do not assume the existence of other long lived particles in the new physics sector, e.g., light 2 The decay length of tau is approximately cτ τ ≈ 50 km · (E τ /EeV), while that of muon reads cτ µ ≈ 4 × 10 8 km · (E µ /EeV), which is much larger than the Earth diameter for a typical energy of EeV. Electrons, on the other hand, will immediately lead to a cascade in the medium after being produced. sterile neutrinos. The structure of the rest of the work is organized as follows. In Sec. 2, we discuss the strategy to propagate neutrinos and taus in matter at extreme energies. In Sec. 3, we introduce the future projects aiming for the detection of GZK neutrinos, and set up the frame for three representative tau neutrino telescopes, i.e., GRAND, Trinity and POEMMA. Sec. 4 summarizes the typical particle physics models that can modify the interactions between neutrinos and matter. Charged Higgs and leptoquark (LQ) models are found to have the largest contributions. In Sec. 5, we illustrate how the new physics modifies the neutrino interaction and study its consequence at tau neutrino telescopes. The sensitivities of GRAND, Trinity and POEMMA telescopes are explored. We find these three telescopes are complementary to each other in probing the new physics modifying neutrino interactions, due to their characteristic experimental configurations and in the way the new physics changes the energy and angular distributions of the tau events. Finally, we make our conclusion in Sec. 6.

Neutrino and Tau Propagation
Neutrinos will be severely attenuated while propagating in matter at EeV energies, around which we expect a bump of the cosmogenic neutrino flux [114][115][116][117][118][119][120][121]. To set the scale, the charged-current (CC) cross section of neutrinos in matter is roughly σ CC ≈ 10 −32 cm 2 · (E ν /EeV) 0.363 [122][123][124], corresponding to a mean free path of L CC ν ≈ 560 km · (E ν /EeV) −0.363 in standard rock with a density of 2.65 g · cm −3 . In comparison, the Earth diameter is around 12742 km. There are two competing effects of having a large neutrino cross section at tau neutrino telescopes: (i) for CC interactions, one can have a larger conversion rate from ν τ to tau and thus collect more events at the detector; (ii) during propagation the neutrino flux will be more attenuated, which implies an opposite consequence. These effects allow us to extract the information of neutrino cross sections at EeV energies, and probe possible deviations from the SM caused by new physics.
As for tau at or above EeV energies, after being produced from CC interactions they will suffer from significant energy loss in medium before they decay [125][126][127][128][129][130][131]. Four different processes contribute to the energy loss of tau, i.e., pair production, photonuclear reaction, ionization and bremsstrahlung. Pair production and photonuclear processes dominate the energy loss at the EeV energy scale. There have been dedicated codes developed to solve the propagation of neutrinos and taus in matter, to our knowledge e.g., nuSQuIDS [132], NuTauSim [128], NuPropEarth [129], nuPyProp [130], PROPOSAL [125] and TauRunner [131,133]. In this work we have developed our own program to propagate neutrinos and taus. The task is to solve a coupled set of integro-differential equations describing the neutrino and tau propagation. We summarize the basic strategy below. The propagation equations can be greatly simplified by taking into account that the neutrino oscillation effect is negligible at energies above a few TeV. Including possible new physics contributions to the cross section and assuming isoscalar nucleon targets, the propagation of neutrino and tau fluxes is governed by the following equation set [126]: tau decay and hard energy loss Here, t is time or equivalently the distance traveled by neutrinos and taus, dΦ ν /dE ν and dΦ τ /dE τ are the differential fluxes of neutrinos and tau, respectively, with the derivative with respect to the solid angle Ω not explicitly shown for convenience, i.e., dΦ ν,τ /dE ν ≡ d 2 Φ 0 ν,τ /(dE ν dΩ), the factor N A is the Avogadro constant, ρ is the mass density of matter and A is the mass number of the target atom. Furthermore, σ CC SM (or σ NC SM ) and σ CC NP (or σ NC NP ) are the CC (or NC) cross sections of SM and new physics contributions, respectively, z ≡ 1 − y represents the fraction of energy going into the final-state lepton with the inelasticity parameter y being the fraction of energy losses, and Γ τ is the tau decay rate. Note that the new physics contributions include possible interference with the SM.
For other unexplained terms in Eqs. (1) and (2), σ pair , σ photo , σ brem , σ ion stand for the energy loss of tau in matter via pair production, photonuclear, bremsstrahlung and ionization,  Figure 3: The outputs of neutrino and tau fluxes as a function of the energy, by injecting a cosmogenic neutrino flux (gray curve) into a 50 km (left-panel) or 2000 km (right-panel) thick standard rock. The black and red curves stand for the neutrino and tau fluxes, respectively. The complete results including all terms are shown as dashed curves, while those neglecting the tau energy loss are given as dotted curves. Note on the vertical axis the differential flux d 2 Φ 0 /(dEdΩ) has been simplified to a label 'Φ' for short. respectively, and β pair , β photo , β brem , β ion are the corresponding parameters for continuous energy loss. For those processes, the cross section of scattering dσ/dy increases rapidly and can even be divergent as the inelasticity y goes to zero, while the energy loss rate ∝ y ·dσ/dy is always under control. For numerical reasons, the tau energy loss should be separated into stochastic and continuous contributions [126]. The stochastic part corresponds to the hard scattering process, where the interaction between tau and medium is dealt with at the cross-section level. The continuous part collects those effects with very small energy loss per scattering (soft). Taking pair production for example, the rate of soft energy loss with dE τ /dt = −ρβ pair E τ reads Here, y cut is a cutoff parameter, above which tau experiences hard scatterings, and below which we can integrate over the inelasticity (0, y cut ) to obtain the continuous energy loss. Note that y cut should not be too large, otherwise the stochastic process may not be fully captured. A reasonable cutoff value would be y cut ∈ (0.001, 0.01) [125].
To solve the integro-differential equations, we can discretize the momentum space into a number of bins, and translate Eqs. (1) and (2) into a large set of ordinary differential equations. In this work, we confine the energy range of interest to be E ∈ (10 6 , 10 12 ) GeV, to which the tau neutrino telescopes are sensitive. The number of bins should be large enough in order to converge to the actual solution of the differential equations. We find that a bin number of 300 is sufficient to have an accurate output of results. Fig. 3 demonstrates the solution of neutrino and tau fluxes after a typical cosmogenic neutrino Table 1: Future neutrino telescopes aiming for the detection of cosmogenic neutrinos. In the fourth (fifth) column, we list the sensitive energy (neutrino flavor) of the telescope. The sensitivity to the neutrino flux strength is given in the sixth column, assuming the flux follows a power law spectrum i.e., E 2 ν Φ ν = constant. Note that for simplicity the label 'Φ ν ' here stands for the differential spectrum d 2 Φ 0 ν /(dEdΩ). Abbreviations for 'mountain-valley' (Mtn-val), 'mountaintop' (Mtn-top), 'atmospheric Cherenkov' (Atm-Cher), 'fluorescence' (Fluo), 'atmospheric radio' (Atm-radio) and 'Askaryan effect' (Aska), have been used. Experiments written in boldface are considered in this work. We should further emphasize that the sensitivity given here is subject to change by the final experimental design, and the ultimate sensitivity can be obtained by properly rescaling from their current proposed exposures. flux from the active galactic nuclei [118] (the gray curve) is injected into a mountain or Earth. The left (right) panel stands for the case that the flux transverses a 50 km (2000 km) thick standard rock with density ρ = 2.65 g · cm −3 . The black and red curves represent the neutrino and tau fluxes, respectively. The dashed curves are the complete results, while the dotted curves are generated without the tau energy loss. In the left panel, for L = 50 km which corresponds to the typical thickness of a mountain, the neutrino flux experiences a negligible attenuation effect. In comparison, the effect of tau energy losses dominates over the tau decay rate, and becomes significant above E τ ∼ 10 8 GeV. In the right panel, the thickness of the standard rock L = 2000 km roughly corresponds to a flux emerging from the Earth with an elevation angle of 10 • . In this case, the attenuation effect of neutrinos becomes considerable at our concerned energies. The output of tau flux is also reduced compared to the left panel owing to the significantly attenuated neutrino flux in medium. An ideal choice of the traveled length is roughly the mean free path of neutrinos, where neutrinos have one scattering on average but are not severely depleted. Most of the information about neutrino interactions is contained in these events.

Tau Neutrino Telescopes
The water-and ice-based Cherenkov techniques [152,153] are not optimized for the detection of cosmogenic neutrinos at EeV energies. These include the past programs DUMAND [154], BAIKAL [155] and AMANDA [156], the running observatories IceCube [1] and ANTARES [157], as well as the future proposals such as IceCube-Gen2 [143], Baikal-GVD [158], KM3NeT [159] and P-ONE [160]. Even though the neutrino cross section increases with energies, i.e., σ ν ∝ E 0.363 , the UHE neutrino flux per decade in energy usually drops faster with a power law spectrum dΦ/d log 10 E ν ∝ E −1 ν . For those cosmogenic neutrinos, a much larger detection volume is thus required, which can be achieved by placing the detectors at a high altitude. It is the hadronic decay of tau in the air from ν τ CC interaction with matter that induces the extensive air shower, which subsequently creates detectable signals like radio waves, Cherenkov light or fluorescence. The detector and the Earth target together form a huge telescope, with an unprecedented effective volume (or area). Depending on the geography of the telescope, they can be typically classified into the following categories: • Balloon-borne telescopes. The ANITA experiment [161,162] has tested the feasibility of neutrino detection via radio waves from extensive air showers as well as the Askaryan effect [163,164] (sensitive to all neutrino flavors). To increase the detection volume, the radio antennas are attached with a balloon, e.g., floating with a height of 30 − 40 km for ANITA. Though no neutrino signal from the Askaryan effect has been detected so far at ANITA [165], there are several extensive air shower events [161,162,166], which are suspected to be tau neutrino candidates. Two of them from ANITA's first and third flights are anomalous [161,166], as their steep incoming angles are in tension with the Standard Model expectation. New physics explanations typically require introducing long-lived degrees of freedom. The fourth flight does not see any anomalous events, but instead four neutrinolike events from the near-horizon direction (Earth-skimming) have been detected [162]. Whether or not these events are truly of neutrino origin remains to be clarified, which is beyond the scope of the present work.
• Space-borne telescopes. The detector can also be attached to a satellite in orbit [87,89], and a much larger detection volume can be achieved compared to the balloon-borne experiment. The space-borne telescope is mainly sensitive to extensive air showers from Earth-skimming neutrino events. To compensate the light attenuation over large baseline, the detection of atmospheric Cherenkov light or fluorescence is usually preferred.
• Mountain telescopes. For the mountain-valley experiment, as proposed by Refs. [90,91], the detector can be placed on one side of a mountain and monitor another mountain over a valley. Furthermore, there are also proposals to place the detector on the top of the mountain, overwatching a thin strip over the horizon, which we will refer to as mountaintop telescope [142] specifically. Compared to balloon-and space-borne telescopes, the ground detector arrays are more extensible and easier to maintain. Almost all detection techniques (radio waves, Cherenkov light and fluorescence) can be utilized for the mountain telescope.    Table 1. The sensitivity curves are reproduced or converted from the original proposals [103-109, 134-141, 143-148, 150, 151] for 90% confidence level, i.e., to collect 2.44 events within a decade of neutrino energy. The existing limits from ANITA [165], Auger [167] and IceCube [168] experiments are given as dotted, dashed and solid gray curves, respectively. The cosmogenic neutrino flux is taken from Ref. [118] for comparison. Right-panel: Our simulated sensitivities to the all-flavor cosmogenic neutrino flux for GRAND200k with 10-year exposure (purple curves), for POEMMA360 with 5-year exposure (blue curve), and for Trinity with 10year exposure (red curve). The thicker (thinner) curve for GRAND is generated assuming the elevation angle of the mountain which hosts the antenna to be β = 3 • (5 • ). The flavor ratio has been set to ν e : ν µ : ν τ = 1 : 1 : 1 with equal fraction of neutrinos and antineutrinos, and one can simply rescale the curve if a different flavor ratio is chosen.
In Table 1, we list to the best of our knowledge the existing and proposed telescopes aiming for the detection of cosmogenic neutrinos. Among these telescopes, the GRAND, Trinity and POEMMA (in Limb mode) experiments represent the mountain-valley, mountain-top and satellite setups, respectively, and have the most outstanding sensitivities to the diffuse tau neutrino flux among similar proposals. We also note that the GRAND proposal has the best diffuse flux sensitivity among these three, and Trinity's sensitivity is slightly better than that of POEMMA in limb mode. Hence we will take them as three representative prototypes in our later analysis. Note that we do not consider to explore the potential of all-flavor neutrino telescopes [112], e.g., those with the Askaryan effect, in this work.
To illustrate the detection possibility, in the left panel of Fig. 4 we summarize the all-flavor sensitivities of these telescopes to the isotropic diffuse flux of cosmogenic neutrinos, along with the prediction of these neutrinos from the active galactic nuclei [118]. Those sensitivity curves are reproduced from the corresponding references in Table 1 by requiring the event number over a decade of neutrino energy interval to be ln(10) · E ν · dN /dE ν = 2.44, which corresponds to 90% confidence level to observe a positive signal [169]. For comparison, the observation time has been unified for similar proposals. To see the current observational status, we recast (with proper rescaling) the existing limits of ANITA [165], Auger [167] (ν τ search) and IceCube [168] as dotted, dashed and solid gray curves, respectively. None of them are able to provide enough sensitivity to the predicted cosmogenic neutrino flux, yet. But in the future, a larger accumulation time or an experimental upgrade for Auger and IceCube might lead to a discovery of cosmogenic neutrinos. In particular, Auger with a scaling of current exposure by a factor of three will have the potential to observe one event at 90% confidence level, given the cosmogenic flux shown in Fig. 4. However, as we will see in later discussions, to study neutrino interactions a sufficiently larger event number is required, which will be challenging for the current running experiments. In the right panel, we give the sensitivity results of our simulations, to be discussed in the following. The notable improvement of POEMMA compared to the left panel should be ascribed to a better experimental configuration [109] than the one adopted in the previous estimate [169]. We will comment on how accurate our simulation is with respect to the published results when we discuss the experiments in detail.

GRAND
One of the major targets of the GRAND experiment [103,104] is to detect the shower of tau decay initiated by ν τ interacting inside the mountain or underneath the Earth horizon. As the extensive air shower propagates in the geomagnetic field a net electric dipole will be developed, which results in strong radio emissions. The radio antenna array of GRAND will be placed on a slope of the mountain which acts as a large projection screen 3 , facing towards another mountain which acts as interaction target. The distance between adjacent antennas will be 1 km, such that with an array of 10000 antennas GRAND10k can cover an inclined surface of 10 4 km 2 . The ultimate stage GRAND200k is to have 20 separate replicates of GRAND10k sites, which can greatly enhance the sensitivity compared to a solo GRAND10k array. With 10-year exposure GRAND200k is able to achieve a world-leading sensitivity E 2 d 2 Φ 0 /(dEdΩ) ∼ 10 −10 GeV · cm −2 · sr −1 · s −1 . This is nearly two orders of magnitude beyond the typical floor of cosmogenic neutrino flux, which indicates a capability of collecting O(100) cosmogenic neutrino events.
An ideal GRAND site should feature both a radio-quiet environment as well as a suitable topography. A possible location of GRAND10k array is at the Tian Shan Mountain, China, where the simulation of neutrino events has been performed [103]. In order to strictly calculate the neutrino event number, the detailed geographical profile near the site is in principle required. However, in the present work we consider the following toy setup to capture the main conclusion of GRAND while keeping a phenomenological simplicity. We assume the antenna array to be uniformly deployed on a 150 × 66 km 2 mountain slope (inclined by β = 3 • from the horizontal), opposing to another mountain with the shortest distance in the valley being 10 km. The target mountain has a height of 2.5 km and a width of 40 km. The length along the valley is set to 150 km, identical to the length of the antenna array. We find that with this toy setup we can closely reproduce the event distributions as well as the sensitivity to tau neutrinos by the official simulation of GRAND.
The antenna array of GRAND10k is wide enough to contain almost all radio pulses from the shower. For instance, the area of a Cherenkov ring imprinted in the antenna screen is around π(l × 1 • ) 2 / sin β ∼ (7 − 46) km 2 , where l = O(50) km is the distance from the shower to the antenna, 1 • is the typical angle of the Cherenkov ring, and β ∼ (3 • − 20 • ) is the yet-tobe-determined elevation angle of the mountain slope with antennas. This area is much smaller than the screen size ∼ 10 4 km 2 . Thus, the simplified methodology assuming a point-like detector, which is much smaller than the scale of a Cherenkov ring, cannot be applied to GRAND. However, it is a good approximation that an event can be accepted if the neutrino's line of sight intersects with the antenna screen.
The total event number can be obtained by integrating over the neutrino flux coming from different angles: where S is the area of antenna screen, θ tr is the angle between the neutrino trajectory and the normal vector of the mountain slope with antennas, and T = 10 year is the data collection time.
The detection probability reads where s is the distance traveled by the tau from the mountain or Earth surface before its decay, and the probability density of the decay reads p decay = Γ(E τ )e −Γ(E τ )s with Γ being the tau decay rate. The probability p det is determined by two criteria. First, the tau converted from neutrino should decay and complete the shower development before reaching the antenna array. Second, the induced signal strength at the antenna should exceed the voltage threshold, which has a conservative value 75 µV or a more aggressive one 30 µV [103]. Namely, with H(x) being the Heaviside-function. The aggressive threshold V min radio = 30 µV will be chosen corresponding to the sensitivity curve put by GRAND. The voltage V radio depending on the detailed antenna response is proportional to the radio flux density which scales as ∝ E τ /l 2 . Instead of a simulation of the radio wave production and the response at antennas, we obtain the voltage by using a scaling relation V radio ∼ V 0 · (E τ /10 8 GeV) · (100 km/l) 2 . The GRAND results can be approximately reproduced from our toy setup if we take V 0 ∼ 2 µV. Note that the dependence of antenna's response on the incoming direction is neglected by using this relation. In practice, all these effects should be taken into account with a dedicated simulation of radio waves, which is however beyond the scope of this work.
In Fig. 4, we have shown all-flavor sensitivity of GRAND200k to the diffuse neutrino flux. The thicker (thinner) purple curve stands for the ten-year sensitivity if the elevation angle of the screen mountain slope is chosen to be β = 3 • (5 • ) in our toy setup. With β = 3 • we are able to reproduce the sensitivity curve in Fig. 4 of Ref. [103] (rescaled to ten years) at a reasonably close level. The event rate benefits from a steeper screen slope, which will possess a larger field of view (FOV) and a better sensitivity to Earth-skimming neutrinos. Furthermore, the event rate at small neutrino energies, e.g. E ν 10 8 GeV, is extremely sensitive to the voltage threshold to trigger the antenna, which is subject to the final experimental design.

POEMMA
The POEMMA experiment [109,169,170] consists of two identical satellites in orbit with an altitude of 525 km. It can operate in two different observation modes: (i) POEMMA-Stereo, aiming for the detection of cosmic rays or neutrinos above 20 EeV via the isotropic fluorescence emission of extensive air shower in the atmosphere; (ii) POEMMA-Limb, for the tau neutrino observation via the Cherenkov light emission of tau decays. We will focus on the POEMMA-Limb mode, which has a much lower energy threshold, E ν 10 PeV, and a better sensitivity at EeV energies than POEMMA-Stereo.
On each POEMMA satellites, there is an optical system which collects and focuses Cherenkov light to the camera. The FOV of the system is 45 • , which can be extended to 360 • in azimuth for the POEMMA360 design. The detection band for Cherenkov light is 200−900 nm for POEMMA. To avoid overwhelming backgrounds from the Sun and moonlight, the Cherenkov camera on the satellite can only operate with a ∼ 20% duty cycle, which is much smaller than for an experiment based on radio waves. This will limit the effective exposure of such telescopes. With five years of observation, POEMMA is able to push the sensitivity to cosmogenic neutrino fluxes down to We will closely follow POEMMA's configuration to generate our events [109]. For a given initial neutrino fluxes, the event number of POEMMA should be obtained with the formula N P = Φ τ P det T cos θ tr dSdΩ tr , which can be written more explicitly as where θ ⊕ is the zenith angle of the tau emergence point with the z axis pointing from the Earth center to the satellite, θ tr is the zenith angle of the tau trajectory but here with z being the vector perpendicular to the Earth surface, and φ tr is the corresponding azimuth angle. The output tau flux only depends on the emergence angle β tr = 90 • − θ tr . The exposure is taken to be T = 5 year × 20%. The probability P det that the Cherenkov light from tau decays can be captured by the detector, is determined by with s being the distance traveled in the atmosphere before the tau decays. There are three requirements that a tau decay event will be accepted by the telescope: (i) the satellite camera (approximated as a point) should be within the Cherenkov angle of the shower event; (ii) the tau decay event takes place within the FOV of the telescope camera; (iii) a sufficient number of photoelectrons in the photomultiplier tube (PMT) should be collected. Hence, the distancedependent probability p det has the following expression where we fix the Cherenkov angle as θ Ch = 1.5 • , and s FOV is the maximal distance within the FOV for a given trajectory. Here, N min PE is the minimum acceptable number of photoelectrons generated by Cherenkov photons registered in the camera, and we fix it as N min PE = 10 following POEMMA. The photoelectron number N PE generated by an extensive air shower can be estimated via the relation [171] with λ being the photon wavelength in units of nm, (λ) the quantum efficiency for the photoelectron conversion, A opt the optical area depending on the angle of the incoming photon, R s the distance from the tau decay point to the telescope. The overall magnitude is fixed by comparing the distributions to Fig. 17 of Ref. [109]. For POEMMA, the integration should be performed over the frequency band (200 − 900) nm. The photon detection efficiency is frequency-dependent, and we will adopt the one for S14520 SiPM array from Fig. 17 of Ref. [109]. The optical area for each POEMMA satellite is 5.71 m 2 on-axis, while for off-axis angles it follows the relation in Fig. 28 of Ref. [109]. The optical depths τ atm and τ aer , by Rayleigh scattering off the atmosphere and by Mie scattering off aerosols, are estimated following Ref. [171], assuming the atmosphere height to be 8 km and the aerosol layer height 1 km. Now we are ready to compute the ν τ events at POEMMA with the tau flux emerged from the Earth surface obtained from the last section. The sensitivity of POEMMA360 with the 360 • FOV in azimuth to all-flavor diffuse neutrino flux is shown in Fig. 4. Five years of observation with a duty cycle of 20% is able to push the sensitivity curve far beyond the threshold of cosmogenic neutrino flux. Our derived sensitivity of POEMMA360 is better than an early evaluation [169], which should be ascribed to the improved configuration considered in Ref. [109], e.g., a larger optical area and a better quantum efficiency of PMTs in the concerned frequency band. However, the potential of POEMMA to the diffuse neutrino flux is weaker than the GRAND200k projection around the cosmogenic flux peak 10 8 − 10 9 GeV.

Trinity
Trinity will deploy the imaging atmospheric Cherenkov telescope on the mountaintop with an altitude of ∼ 2 km [105][106][107][108]. One may think of Trinity as a ground analogue of POEMMA, so most of the considerations to derive neutrino events at POEMMA can apply directly to Trinity after modifying some of the experimental parameters. The complete configuration of Trinity will be made up of three stations, each of which has an imaging system with 360 • field of view in azimuth and 5 • in zenith. The horizon of the Trinity telescope with 2 km height is around 88.56 • from the vertical. A primary investigation of Trinity has found that with 2 • FOV above and 3 • below the horizon, one can achieve an excellent sensitivity for such setups [105]. The effective collecting area of each Trinity mirror for Cherenkov light is 10 m 2 , which can be triggered with a minimum number of N min PE = 24 photoelectrons to sufficiently reject the background. Similar to POEMMA, the duty cycle is limited to 20% each year.
Because of the low altitude of Trinity compared to that of POEMMA, Cherenkov photons can be seen by the mirror even far outside the 1.5 • Cherenkov cone. For instance, a 10 8 GeV shower that is 135 km away from the telescope can have sufficient photoelectrons to trigger the mirror, even if the telescope is 10 • away from the shower axis. To derive the Trinity events, the major difference from POEMMA is in Eq. (9): That is, we do not require the telescope to be within the 1.5 • Cherenkov cone. For the photoelectron number N PE , we adopt fitted functions in Ref. [105] for a height of 2 km. Our calculation of the 10-year sensitivity of complete Trinity setup is given as the solid red curve in the right panel of Fig. 4, which is reasonably close to the Trinity official result [107] in the left panel.

Minimal New Physics Scenarios
In this section we summarize possible new physics contributions that directly modify neutrinomatter interactions at the tree level, which can be probed at tau neutrino telescopes. We maximize the effect to see how much the new physics can contribute to the deviations from the SM. As has been mentioned, the new physics is involved via the diagrams in Fig. 2, where we must have a neutrino and a matter constituent (e, q or g) in the initial state and a lepton (ν e,µ,τ or τ ) in the final state. The CC production of e and µ final states is not relevant for tau neutrino telescopes considered here. Based on the above consideration, we give in Fig. 5 the inclusive new physics cases, where the gauge symmetry and baryon number are always conserved but the lepton number can be possibly violated by two units. We have shown all possible contributing diagrams that modify the neutrino-matter interactions at tree level, which include both neutrino-nucleon and neutrinoelectron collisions. One can generalize these scattering processes to a model-independent framework of effective operators by integrating out the heavy degrees of freedom. However, we point out that the leading contribution for some processes involves resonance production, which is not a trivial task to include in the effective field theory [172].

Neutrino-nucleon collision
The center-of-mass energy of neutrino-nucleon collision reads √ s ≈ 2E ν M ≈ 43 TeV E ν /EeV for an incoming neutrino with energy E ν and nucleon with mass M ≈ 0.94 GeV. The collision at very high energy scales benefits from the increasing number of sea quarks and gluons in the proton. For processes relevant for the neutrino-nucleon scattering, we find the following models: • Charged Higgs (H ) -There are many SM extensions such as the two Higgs doublet model [173], supersymmetric models, left-right symmetric model [174], the type-II seesaw Figure 5: The inclusive contributing diagrams that modify the neutrino-matter interactions as in Fig. 2 at tree level. Loop-suppressed processes are ignored in this work. The symbol 'l' in final states stands for the charged lepton or the neutrino. Note that we assume that the baryon number is preserved, ∆B = 0 up to a very high energy scale, and the lepton number can be preserved or violated by two units, i.e. ∆L = 0, 2. The new physics diagrams that can induce significant signals at tau neutrino telescopes are highlighted in boldface. model [175][176][177][178][179], radiative neutrino mass models [180][181][182][183], axion models [184,185] and dark matter models [186] possessing charged Higgs bosons in their particle spectra. A prototypical example can be taken as the Zee model [187], which is one of the most wellknown neutrino mass models for generating neutrino masses and mixings radiatively at the one-loop level. We consider the following interaction forms where y q ij and y αβ are the Yukawa coupling constants for quarks and leptons, respectively, {i, j} are the quark flavor indices with U ≡ {u, c, t} and D ≡ {d, s, b} including both left-and right-handed fields, and {α, β} are the lepton flavor indices with E ≡ {e, µ, τ }. In this work we focus on the couplings y q cs and y τ τ , which are not strictly constrained by laboratory searches.
The combination of y q cs and y τ τ will switch on the CC conversion from ν τ to tau, i.e., ν τ + s(c) → τ + c(s). The modification to neutrino cross section with M H = 90 GeV and y q cs = y τ τ = 1 is shown in Fig. 6. Other quark couplings, e.g., y q ud , are more severely constrained by the experimental searches. In fact, at the energy scale relevant for tau neutrino telescopes, y q cs can lead to comparable effects as y q ud because of the increasing number of cc pairs in the nucleon. As for other leptonic couplings, y αβ with α, β = e, µ are not relevant for tau neutrino telescopes. Furthermore, y τ α with α = e, µ will contribute to the process ν τ → e, µ during neutrino propagation, resulting in only a depletion in the tau neutrino flux. The effect of y ατ with α = e, µ is very similar to y τ τ but with the conversion ν e,µ → τ . To be definite, we only switch on y τ τ for the later numerical analysis.
Here we mainly focus on leptoquarks that will induce observable signatures at tau neutrino telescopes. Since it involves a neutrino in the initial state, there are only four leptoquark possibilities which are denoted as S 1 (3, 1, −1/3), S 3 (3, 3, 1/3), R 2 (3, 2, 7/6),R 2 (3, 2, 1/6) in the literature. Here we are specifying the leptoquarks by their SM quantum numbers {SU(3) C , SU(2) L , U(1) Y } and the electric charge is defined as Q = T 3 + Y. The possible interaction forms for the initial neutrino are qS 1 , qS 3 , D cR 2 , U c R 2 at the renormalizable level, with denoting the lepton doublet and q the quark doublet. Note again that we concentrate on scenarios where there is no extension in the fermionic spectrum in addition to SM fermions. Taking the S 1 leptoquark as an example, the relevant Yukawa Lagrangian can be expressed as: where y LL and y RR are Yukawa couplings, and V represents the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix. For the S 1 leptoquark, there are two channels contributing to the neutrino scattering: ν +D → LQ → τ +j / ν +j and ν +g → LQ+j → τ +2j / ν +2j. Considering collider bounds, we assume the couplings are dominant among the q = s and = τ in order to maximize the possible contribution. See the later discussion in this section for more details. The cross section of S 1 LQ production with M LQ = 400 GeV is given in Fig. 6.
• Charged gauge boson (W ) -Similar to the charged Higgs, a heavy new charged gauge boson W can induce similar signatures at tau neutrino telescopes. However, due to the tight constraint from existing experiments [201], i.e., M W > 5 TeV, such a heavy W does not lead to any significant effects at tau neutrino telescopes. There are a few models such as Ref. [202], where the W mass can be much lighter (∼ 1.8 TeV) under some broad-width assumption, however we find that this is still not enough to give observable imprint at tau neutrino telescopes; see the left panel of Fig. 6.
• Neutral gauge boson (Z ) -A new neutral gauge boson Z exists in many extensions of SM in the gauge sector, such as gauged L µ − L τ [203][204][205][206], B − L models [189,[207][208][209][210], dark sector models [211][212][213], and left-right symmetric model [174]. To have an observable neutrino-proton scattering process, a Z should be coupled to both quarks and leptons, for example with the interaction This will contribute to the NC scattering process: ν α + q → ν β + q. With the same cross section, the effect of NC process is weaker than the CC one, because for the t-channel exchange of vector boson the final-state neutrino takes away most of the energy of the incoming neutrino. This results in an ineffective scattering other than possible mixture among different neutrino flavors.
In addition, such Z coupled to quarks are tightly constrained by collider searches. If Z talks to quarks of the first family as well as leptons, the lower bound of its mass reads M Z 5 TeV [214]. A minimal scenario is that Z only talks to the third family, for instance the U (1) B−L model [215]. The searches of bbτ τ final states can place a limit M Z 100 GeV for order one couplings [216]. However, for an anomaly-free gauge extension consistent with SM [215], Z will also couple to the first and second generations of quarks via mixing, and the mono-jet searches at LHC set a stringent constraint M Z 500 GeV [217] with order one couplings. In the left panel of Fig. 6, we show the case for a Z coupled primarily to bottom quark and tau neutrino with M Z = 500 GeV and an O(1) coupling.
• Neutral scalar (φ) -Various BSM theories predict neutral scalars in their particle spectra, such as the two Higgs doublet model [173]. Similar to the Z scenario, a neutral BSM scalar φ should be coupled to both quarks and leptons to have an observable neutrinonucleon scattering process. However, similar to Z , if φ couples to quarks of the first two generations, the coupling is tightly constrained. In addition, if φ inherits flavor violating couplings, it will give rise to flavor changing neutral-current processes, which is stringently bounded [218]. This leads us to consider a third-generation-philic scalar scenario, where φ mostly couples to bottom quarks and tau neutrinos. This scenario can be constrained by pp → bb + / E T searches [219]. We find that this bound is not very significant. In the left panel of Fig. 6, we show the case of φ coupled to b and ν τ with M φ = 90 GeV and an O(1) coupling. Because at large Q 2 there are more and more sea quarks (including heavy quarks like c and b) with small momentum fraction, the contribution of a neutral scalar coupled to b quarks is only a factor of two smaller than charged Higgs coupled to s and c quarks. After considering the above new physics scenarios, we find only the charged/neutral Higgs and leptoquark models can lead to detectable signals of neutrino-nucleon scatterings at tau neutrino telescopes, if the laboratory constraints are taken into account. In the left panel of Fig. 6, we depict neutrino-nucleon cross sections for the pure SM case (gray curve) and inclusive new physics scenarios. The uncertainties of SM cross section induced by the 1σ PDF errors computed from the CT18 set [85] are also shown for comparison, as the gray band. In the following, we shall elaborate on the charged Higgs (the neutral Higgs scenario is similar but weaker) and leptoquark models, which will be further explored in the rest of the work. First we briefly discuss the limits on charged Higgs from collider searches. The s-channel Drell-Yan process, mediated by either γ or Z boson, can pair-produce charged scalars at the LEP experiment. Here we are concentrating on the scenario where the charged scalar mostly couples to τ ν and cs. Hence, the branching ratios of H decaying to ντ and cs modes should add up to one. In the left panel of Fig. 7, we show the limits from charged Higgs searches [220] at LEP by looking at final state signature cscs or τ ντ ν. At LHC, charged Higgs can similarly be produced in pairs via Drell-Yan processes. After being produced in pairs, if the charged Higgs decays back to τ ν final state, it will be further constrained from supersymmetric stau searches [221] because both of them give rise to the same final-state signature in the zero neutralino mass limit. However, we find that in our scenario, where the charged Higgs is emerged from the Zee model similar to Ref. [202], the observed cross-section limit is still larger than the theory prediction (for detailed collider analysis, see Ref. [202]). Note that there are several dedicated charged Higgs searches [222][223][224], however in all these scenarios a large H tb coupling has been considered. Since in our scenario, the charged Higgs only talks to cs, these bounds are not directly applicable. On the other hand, since the luminosity of cs quarks is two orders of magnitude lower than that of ud quarks at LHC [225], other search limits are very weak in our scenario.
We continue with the experimental constraints on leptoquark models. Leptoquarks can be pair produced copiously via gg and qq fusion processes at the LHC, and the pair-production rate is uniquely determined by the LQ mass, irrespective of their Yukawa couplings. Leptoquarks can also be singly produced in association with leptons through s-and t-channel qg fusion processes, and the production rate depends on both LQ mass and LQ Yukawa couplings. However, LHC constraints from single-production of LQ are not so severe compared to the pair-production limits [226,227] unless the Yukawa couplings to the first and second-generation quarks are too large (> 1) [180,228]. Other than these direct limits from pair and single production of LQ, there are indirect limits on the mass and Yukawa couplings of LQs from the dilepton searches at the LHC [229] because the process pp → l + l − gets significantly modified due to the t-channel LQ exchange. However, due to our judiciary choice of Yukawa texture relevant for tau neutrino telescopes, only the process pp → τ + τ − gets modified, keeping the pp → e + e − /µ + µ − processes almost unaltered. Whereas, the LHC limits are not so strict for pp → τ + τ − signature [230]. The most stringent limit comes from LQ pair production searches. As we previously mentioned, we are interested in the scenario where LQs mostly decay into τ j or/and νj final state.
Hence, the searches for the final-state signatures containing two neutrinos and two jets (ννjj) and two tau leptons and two jets (τ jτ j) impose the best constraints. Although there is a dedicated search for νjνj final-state signatures [231], there are no such dedicated searches for τ jτ j final states. We recast the limits for τ jτ j final-state signatures [198] from τ bτb searches [227] considering a b-jet misidentification rate of 1.5% as light jets with a b-tagging efficiency of 70% [232]. The green and purple shaded regions in the right panel of Fig. 7 depict the bounds from τ jτ j and νjνj final states. We can see that the former is less stringent than the latter one, and the Yukawa couplings can be chosen in such a way that our LQ can be as light as 400 GeV, making it consistent with all the experimental searches for a branching ratio to τ j of 80% and that to νj of 20%. The branching ratio arrangement requires the LQ Yukawa coupling of τ j to be stronger than that of νj, and we see this is natural for S 1 and R 2 LQ models. Whereas for S 3 andR 2 LQ models, the Yukawa couplings to neutrino and charged leptons are the same, leading to 50% νj and τ j and hence a lower bound M LQ > 700 GeV; see Appendix for details. The lower limit of LQ mass can be further reduced by allowing other channels to partly share the branching ratios. The consequence of leptoquarks on km 3 -scale water-or ice-based Cherenkov detectors has been explored in previous literature [233][234][235][236][237][238][239][240][241], but the bound and future sensitivity of IceCube have been found to be difficult to exceed the current LHC bounds [241].

Neutrino-electron collision
There might be leptophilic forces which can evade hadron collider limits while affecting the neutrino-electron scattering at tau neutrino telescopes. The impact can be enhanced by the resonant production of exotic charged particles from neutrino-electron scatterings, e.g., ν + e − → ∆ − in type-II seesaw scenario or ν + e − → H − , W − [242,243], similar to the Glashow resonance in the SM [28]. For EeV incoming neutrinos, the resonance is at M H = √ s ≈ 1 TeV. The cross section for the process ν α + e − → H − reads [242] where Γ H = Y 2 M H /(16π) is the decay width with Y being the coupling strength. The resonance will be smeared due to the spread of initial neutrino flux and the energy distribution of final states.
Integrating over a decade of neutrino energy around the resonance, we obtain the flux-averaged cross section as 2 × 10 −34 cm 2 for M H = 1 TeV and Y = 1. This should be compared to the CC cross section σ CC (1 EeV) ≈ 10 −32 cm 2 . The smallness of resonance enhancement can be ascribed to the relation σ res ∝ M −2 H . Hence, we cannot gain much from the resonance at the energy scale of cosmogenic neutrinos, which is consistent with the result in Ref. [110].
For lower mediator masses, the effect of resonance becomes more important and may dominate over the CC process similar to the Glashow resonance at E ν = 6.3 PeV. However, on the one hand, both the flux and sensitivity of tau neutrino telescopes drops at lower neutrino energies. On the other hand, the limits of non-standard neutrino interactions (NSIs) will come into play if there is a large exotic coupling between electron and neutrino [242]. These limits include the LEP experiment, BOREXINO, IceCube, etc. Among them, the strongest one is found to be from IceCube atmospheric neutrino data [244], which set y/M H < 0.2 (100 GeV) −1 [242]. For a lower incoming neutrino energy, say E ν = 0.1 EeV, the resonance is at M H = √ s ≈ 319 GeV with y < 0.6, the cross section averaged over a decade in energy turns out to be 8 × 10 −34 cm 2 , negligible compared to the SM CC one σ CC (0.1 EeV) ≈ 4 × 10 −33 cm 2 .
In the meantime, the t-channel exchange of a vector mediator appreciates the forward divergence. Taking the charged mediators H − and W − for example, the cross sections for ν τ + e − → e − + ν β are found to be For small mediator masses, the vector case (W − ) features a forward enhancement, i.e., σ ∼ g 2 τ e g 2 βe /(4πM 2 W ) with M W √ s. This, however, does not apply to the spin-flipping scalar mediator (H − ). The NSI limits will also apply to the t-channel processes. Recalling y/M < 0.2 (100 GeV) −1 , the optimal cross sections can be estimated by taking g τ e = g βe = 1 and M H = M W = 500 GeV, which yields σ W ∼ 8 × 10 −35 cm 2 and σ H ∼ 2 × 10 −36 cm 2 for E ν = 1 EeV. Similar conclusions can be made for neutral mediators, e.g., Z . But, in this case the forward enhancement of vector force (i.e., small momentum transfer) does not noticeably alter the final-state neutrino energy, which weakens the effect significantly. Hence, we conclude that both s-and t-channel neutrino-electron scatterings do not play significant roles at tau neutrino telescopes. We give in the right panel of Fig. 6 cross sections of s-and t-channel processes with optimal mediator choice saturating the experimental limits. This visually demonstrates the negligible effect compared to hadronic processes, which are irreducible backgrounds for the neutrino-electron scattering.

Numerical Results
After implementing the neutrino and tau propagation algorithm into GRAND, POEMMA and Trinity setups, we are ready to investigate the modified neutrino interaction induced by new physics effects. For the GRAND experiment, we assume the matter in the mountain and Earth is composed of uniform standard rock with density ρ = 2.65 g · cm −3 . This is a good approximation for GRAND, as the largest depth in Earth reachable by Earth-skimming neutrinos is around 24 km with an elevation angle of 5 • . For POEMMA, we need to solve the propagation equations with a realistic Earth matter profile, for which we will adopt the PREM model.

New physics effects
Because of the decay and energy loss of taus, the largest length scale L τ that a tau can travel in medium is around 50 km for all energies. Thus, the tau events registered in the telescope are essentially produced within a layer of thickness 50 km close to the matter surface. Beyond this thickness, taus converted from ν τ can only play a role in the regeneration of ν τ via tau decays. In contrast, neutrinos at EeV energies possess a much longer mean free path of L CC ν ≈ 560 km · (E ν /EeV) −0.363 L τ ≈ 50 km in the standard rock. Hence the final event number is basically dependent on the nearly-constant neutrino flux strength within L τ underneath the matter.
The neutrino event number can be roughly estimated with a relation N ≈ Φ out ν · σ CC τ · L τ , where Φ out ν is the neutrino flux near the surface, σ CC τ is the tau production cross section and L τ is the tau existing length scale. Considering two well-separated length scales of taus and neutrinos, the presence of new physics can alter the tau neutrino event number through two possible effects, as has been previously noted. First, the new physics enhancement of CC cross section will increase the production rate of taus in the layer close to the surface, e.g., by ∼ δσ CC τ ·L τ with δσ CC τ denoting the possible deviation induced by new physics. Second, the neutrino flux will experience attenuation before reaching the matter layer near the Earth surface, by a factor Φ out ν ∝ exp (−σ tot · L ν ), where L ν is the distance traveled by neutrinos in matter and σ tot is the total cross which depletes the neutrino flux. Note that σ tot contains both the CC interacting producing taus as well as processes not related to the tau production, e.g., the NC interaction. The modification in the final event number due to new physics can be estimated by where Φ in ν is the initial neutrino flux and δσ tot is the modification of total neutrino cross section due to new physics. Whether we have an increased or decreased event number depends on the sign of δσ CC τ −σ CC τ L ν δσ tot , which can be simplified to 1−σ CC τ L ν in the assumption of σ tot = σ CC τ . In practice, the telescope will monitor a wide field of view corresponding to a range of neutrino baselines L ν in matter. We can therefore resolve the presence of new physics by exploring the variance of events as a function of L ν .
The energy and angular distributions of event number in GRAND, POEMMA and Trinity for two new physics scenarios are illustrated in Fig. 8. In the top two panels, we give the event distributions for GRAND200k assuming the elevation angle of antenna screen is 3 • . The middle and bottom panels show the distributions of POEMMA360 and Trinity (three stations), respectively, and both of them have the full FOV in azimuth. In all panels, the black curves stand for the SM scenario, while the blue and red ones are for the charged Higgs with M H = 90 GeV and y cs = y τ τ = 1 and leptoquark with M S 1 = 400 GeV and y LL 1,sτ (cτ ) = y RR 1,cτ = 1, respectively. In each panel, both the differential distribution (upper one) with respect to the tau energy or the elevation angle and the binned event number with error bars (lower one) are given. In the following, we make some observations on Fig. 8.
• With the given setup, GRAND200k can collect much more events than POEMMA360, as one can expect from their sensitivities to diffuse neutrino flux in Fig. 4. For the angular distribution in the right panel, the elevation angle α < 0 • corresponds to the event coming from above the horizon, which transverses only the mountain target. Because the typical mountain thickness ∼ 100 km is much smaller compared to the neutrino attenuation length, events with α < 0 • are very useful to normalize the initial neutrino flux. On the other hand, an event with α > 0 • (Earth-skimming neutrinos) corresponds typically to a longer traveling distance. The attenuation effect is what we will use to extract the information of the neutrino cross section. On average, the additional contributions from leptoquark and charged Higgs models will increase the event rate at GRAND, and distort the angular  distribution of events. Compared to GRAND, Trinity can collect more Earth-skimming neutrinos.
• The angular distribution of POEMMA events reaches the maximum around the elevation angle α = 4 • , and then it drops rapidly as we go to larger elevation angles due to the attenuation effect. Note that for POEMMA there is a one-to-one correlation between the elevation angle and the distance traveled by neutrinos in Earth. In comparison to GRAND, the new physics contribution reduces the event number significantly at POEMMA. This is due to that the majority of FOV from the POEMMA satellite corresponds to very long chord lengths in the Earth. The event modification due to new physics is negative in Eq. (18) with L ν 1/σ CC τ .
• The energy distribution of events at a telescope is determined by two factors: (i) the input of initial neutrino flux, which has been taken from Ref. [118]; (ii) the average distance traveled by neutrinos. Because neutrinos interact more strongly at higher energies, one can observe a shift of event to lower energies as the neutrino flux has experienced multiple scatterings over a very long distance. This effect can be reflected by the comparison between energy distributions of GRAND and POEMMA, the latter of which has a longer average baseline.
It should be stressed that the GRAND, POEMMA and Trinity experiments are remarkably complementary, because the new physics enhances the event rate at one experiment but reduce the rate at another. Their combination will be very useful to maximize the sensitivity to new physics. Next, we shall investigate the new physics effect in a statistically quantitative approach.

Sensitivities
The largest systematic uncertainty stems from the unknown priors of the diffuse neutrino flux. The magnitude and shape of the cosmogenic neutrino flux is partly model-dependent, resting on, e.g., the evolution model of cosmic sources and the initial chemical component of cosmic rays. To be conservative, we can assume that the initial tau neutrino flux is completely unknown, which is to be derived from the same data set at tau neutrino telescopes that we use to probe the new physics models. This induces a degeneracy between the new physics effect and the initial flux input. For instance, a larger or smaller input of initial neutrino flux can mimic the effect of event excess or depletion due to new physics. The key to resolve the degeneracy is relying on the angular distribution of neutrino events at tau neutrino telescopes. In Fig. 9, we show as an example the event distributions at GRAND by varying the initial neutrino flux such that the deviation between the SM and new physics is minimized. We find that the new physics effect in energy distributions can be completely compensated by varying the input of the unknown initial neutrino flux. In comparison, the difference in angular distributions is stable against the change in the initial flux. In order to quantify the flux uncertainty, we split the initial neutrino flux into 24 bins in the energy range of (10 6 , 10 12 ) GeV. Four bins form a decade in energy, which should be fine enough to control the variation in the cosmogenic flux 4 . We then let the magnitude of neutrino flux in each energy bin vary freely while fitting the simulated data. The further observation of cosmic rays will no doubt enhance our knowledge of the nature of cosmic rays, and hence reduce the systematical uncertainty of neutrino flux. As another benchmark, we shall also explore the sensitivity by assuming the initial neutrino flux to be completely fixed by cosmic ray observations. The true sensitivity in the future experiment should be in between this ideal case and the conservative treatment above.
Another systematic error originates in the PDF uncertainties. This can be conveniently evaluated by using the available PDF set with errors [245]. To accommodate this uncertainty, we repeat the whole computation for all the 59 PDF sets from the CT18 PDFs [85], and take the difference in each computation as the error induced by the corresponding PDF. The uncertainty of SM CC cross section in Fig. 6 is calculated in this way, and a similar procedure will be performed for the final events.
The energy resolution of POEMMA is subject to the long distance from the emerging point of tau to the satellite camera. Owing to the rapid scattering of photons off atmospheric molecules and aerosols, POEMMA has a very limited energy resolution to the extensive air shower induced by tau decay. The energy resolution of tau for 1 EeV is not provided by the POEMMA collaboration, but we can estimate it according to the empirical relation ∆E/E ≈ 1/ N PE 30% [171]. The fluctuation of air shower development might worsen the energy reconstruction. A conservative choice of the resolution ∆E/E = 100% will be adopted in this work. The angular resolution of POEMMA is limited by the Cherenkov cone, typically 1.5 • in the view of the satellite camera. This can be translated into the resolution on the tau elevation angle according give the results without considering the uncertainty of initial neutrino flux. For comparison, we also give the perturbative limit of Yukawa couplings, the LEP [220] and LHC [230] constraints on the lower masses of charged Higgs and leptoquark, as well as the leptoquark constraint from LHC high-p T dilepton tail searches.
to the relation α = arcsin (R ⊕ + h s )/R ⊕ sin θ s − 90 • , with θ s being the viewing angle of the satellite with respect to the vertical. For ∆θ s = 1.5 • near the horizon, the difference in tau elevation angle is as large as ∆α ∼ 8 • , due to the very high altitude of POEMMA satellite. The ground-based GRAND and Trinity experiments outperform POEMMA in both energy and angular resolutions. For instance, the energy resolution of GRAND can be as good as 15%, and the angular reconstruction can achieve a level of sub-degree (less than 0.5 • ) on average [103].
By summing over the two-dimensional grid of energy and elevation angle, we obtain the minimum of χ 2 via where n exp i and n th i are the nominal experimental event number and the theoretical prediction from new physics models, respectively, in the i-th bin. The total bin number in the energy and angle grid is N bins . Note that for our current sensitivity analysis, we do not include higher order processes which in our analysis will appear as part of the theoretical systematics [113]. To generate the experimental data n exp i , we use the cosmogenic neutrino flux in Ref. [118] as the input, and adopt the SM cross section without new physics, i.e., assuming the SM as the true model. For each given parameter choice of new physics, the theoretical expectation n th i can be calculated with a randomly given diffuse neutrino flux. The minimum of chisquare, χ 2 min , will be obtained by scanning over initial neutrino fluxes, which is expected to reduce the statistical significance due to the degeneracy between flux and new physics effects. The value of χ 2 min can then be used to constrain new physics models regardless of priors of the diffuse neutrino flux. In Fig. 10, the purple, blue and red solid curves represent the sensitivities of GRAND, PO-EMMA and Trinity, respectively, to the charged Higgs (left panel) and leptoquark (right panel) models. The combined sensitivity of those three experiments is given as the black curves. The existing laboratory constraints on the scenario we are considering are also presented for comparison. We observe that these three telescopes can have sensitivities surpassing the current collider limits. Some further comments on the results are given below.
• The mountain-based telescopes GRAND and Trinity have comparable sensitivities to the charged Higgs and leptoquark models, if the initial diffuse neutrino flux is given as a wellknown theoretical prior. However, as for the case with unknown initial flux, Trinity has the best sensitivity. Though the effective exposure of GRAND is the largest among three telescopes, many events of GRAND are coming from the direction above the horizon (penetrating only the mountain), where the degeneracy between the new physics contribution and flux uncertainty is difficult to resolve. Whereas, Trinity collects exclusively the Earthskimming neutrinos which contain more information about neutrino interactions. Since the final sites of GRAND are yet to be determined, our result here only represents a special scenario with the antenna screen being deployed on a 3 • inclined mountain. A steeper host slope for GRAND, which has a wider FOV for Earth-skimming neutrinos, will certainly improve the results here.
• The POEMMA experiment is subject to the high altitude of satellite, i.e., 525 km, in comparison to the mountain-based telescope ∼ 2 km. A higher altitude will worsen the resolution of neutrino elevation angle, while the new physics mostly manifest itself by altering the angular distribution of events. Therefore, the ground-based Trinity setup seems to be more optimized in probing the neutrino interaction of our concern.
The combination of GRAND, POEMMA and Trinity greatly enhances the sensitivity. Instead of a simple sum of χ 2 min , these three telescopes are complementary to each other in resolving the diffuse flux uncertainty and probing the new physics effect. The above discussions are made for the case with the initial neutrino flux unknown. As a comparison, we also give the ideal case where the initial neutrino flux is completely known and not minimized over as the dotted curves in Fig. 10.

Conclusions
Other than supplementing the multimessenger astronomy, tau neutrino telescopes can also play the role of a particle collider which collides a high-energy neutrino beam with proton and electron. We have systematically investigated the new physics scenarios that modify the neutrino-matter interactions relevant for tau neutrino telescopes including: charged and neutral Higgs, leptoquark, as well as neutral and charged gauge bosons. These extended scenarios have already been under tight constraints from ground-based colliders like LEP and LHC. Among them, we find the charged/neutral Higgs and leptoquark can have significant imprint on the neutrino-proton scattering, if the existing experimental limits are considered. In particular, this will require the new particles to exclusively have large couplings with the second or third family, for which the related processes can be suppressed at LEP or LHC. Tau neutrino telescopes probe the neutrino-proton COM energy as high as 45 TeV, where the PDFs of second and third generations of quarks are not that suppressed and in some case can be comparable to that of u and d. The gauge bosons are mainly subject to feasible model construction which constrains the lower mass together with collider searches.
By solving the neutrino and tau propagation equations, we have generated the events at GRAND, POEMMA and Trinity with configurations close to their realistic experimental setups. Their sensitivities to the new physics scenarios are limited by the unknown prior of the cosmogenic neutrino flux, which might be improved with future observations of cosmic rays. Under the assumption of completely unknown flux priors, we find that the angular distribution of events carries most information about the absolute cross section of neutrinos. With a χ 2 -analysis, we have obtained the sensitivities of three representative tau neutrino telescopes to the parameter space of charged Higgs and leptoquark models, i.e., Fig. 10. We should keep in mind that those sensitivity curves are expected to shift according to the final design of these three telescopes.
Tau neutrino telescopes, as a particle collider, feature highest neutrino-matter colliding energies. However, we need to point out that their sensitivity to the new physics scenario (in particular, those which have been explored in the work) is more or less restricted for a number of reasons. First, the initial neutrino flux is not under control and spreads over a wide range, compared to the nearly monoenergetic collider beam, which will partly erase the potential new physics signatures. Second, even though we also have an electron target in matter, neutrino-proton collision will produce a large irreducible background. This restriction limits the new physics searches through the clean pure leptonic portal. Third, the event topology at tau neutrino telescopes is not as broad as in the collider spectrometer, and the detailed product of neutrino-matter collision is not a direct observable.

A Neutrino Cross Section: Standard Model
The expressions of neutrino cross sections in the SM are partly available in the literature [122][123][124]255,256], but we present them here for completeness. For a general process ν +N → l+q f , the ultimate differential cross section with respect to the Bjorken scaling variables can be obtained by converting from the parton-level cross section in the center-of-mass frame. In the following expressions, we assume the most general cases which do not require the final states l and q to be massless. For q = u, d, s, c and b, it is a good approximation to set m q = 0 for our concerned neutrino energy. But for top quark, the production threshold corresponds to a Bjorkenscaling variable m 2 t /(2E ν M ) ∼ 1.6 × 10 −5 for E ν = 1 EeV, which cannot be neglected at energy scales of our interest. To incorporate the effect of final-state hadron mass, the usually adopted Bjorken-scaling variable x in PDFs should be replaced by x in the slow-rescaling prescription [122,123,257]:ŝ where x = x + m 2 q /(2M E ν y) is the modified Bjorken variable, y = (1 − E l /E ν ) is the inelasticity, s is the square of total energy in the parton COM frame, and Q 2 = −t is the square of momentum transfer. In the limit of m q → 0, we recover the massless case x = x.
At the parton level, the charged-current cross section between neutrino and quark parton in the COM frame reads for the down-and up-type quarks D and U , respectively, where u = All m 2 i −ŝ−t. Incorporating PDFs for the quark parton q(x ), we have The inelasticity y is connected to the zenith angle in the parton COM frame with y =ŝ + m 2 q − m 2 l − λ 1/2 ( √ŝ , m l , m q ) cos θ 2ŝ , where 4ŝ| p l | 2 = λ( √ŝ , m l , m q ) = (ŝ − (m l − m q ) 2 )(ŝ − (m l + m q ) 2 ) with p l being the momentum of l in the COM frame. By converting all Mandelstam variables into the Bjorken x and inelasticity y, the final cross section takes the form d 2 σ νN dx dy = d 2 σ νN dx d cos θ · ∂(x , cos θ) ∂(x , y) .
One can check the correctness by considering the extreme case m l = 0 andŝ → m 2 q , such that q is produced nearly at rest in the COM frame. In this case, the hadron q should take away almost all the initial neutrino energy, i.e., y → 1.
The neutral-current cross section at the parton level can be similarly obtained dσ νq d cos θ = where g V = 1/2−4 sin 2 θ w /3 and g A = 1/2 for q = U , and g V = −1/2+2 sin 2 θ w /3 and g A = −1/2 for q = D, with sin 2 θ w ≈ 0.231 being the weak mixing angle. The rest of the derivation is similar to the charged-current case.

B Leptoquark Models
The relevant Yukawa interactions of S 1 LQ with matter is as follows: where V represents the CKM mixing matrix, y 1ij denote the elements of an arbitrary complex 3×3 Yukawa matrix, and the flavor and SU (2) indices are denoted by i, j = 1, 2, 3 and a, b = 1, 2, respectively. Here we set the Yukawa texture such a way that S 1 LQ can dominantly decay via the following modes: S 1 → τ c, ν τ s. In this way, the bounds on the LQ mass will become less stringent, i.e., M LQ 400 GeV, if the branching ratios to τ c and ν τ s are set to 80% and 20%, respectively.
Similarly, for S 3 LQ, the relevant Yukawa Lagrangian can be written as: where τ k , k = 1, 2, 3 denotes Pauli matrices. The branching ratios of S 3 LQ are the same for neutrino and charged lepton final states, which implies a lower limit M LQ 650 GeV from the searches of LQ pair production. For R 2 LQ, the relevant part of the Yukawa Lagrangian can be expressed as: There are two states of R 2 LQ, one with electric charge Q = 2/3 (R 2/3 2 ) and other one with Q = 5/3 (R 5/3 2 ). Since R 2/3 2 state interacts with neutrino and up-type quark as well as charged lepton and down-type quark, it provides a significant imprint at tau neutrino telescopes. For R 2 LQ, the Yukawa structure is chosen in such a way that the state R 2/3 2 from R 2 LQ can dominantly decay to the following modes: R 2/3 2 → τ s, ν τ c. ForR 2 LQ, the relevant Lagrangian is given by R 2 LQ also comprises of two states with electric charges Q = 2/3 and −1/3. However, we find that due to the Yukawa structure,R −1/3 2 coupled to neutrinos dominantly decays toR −1/3 2 → ν τ d, ν τ s, which requires a LQ mass 1 TeV at least to satisfy existing collider constraints.