The ultra-high-energy neutrino-nucleon cross section: measurement forecasts for an era of cosmic EeV-neutrino discovery

Neutrino interactions with protons and neutrons probe their deep structure and may reveal new physics. The higher the neutrino energy, the sharper the probe. So far, the neutrino-nucleon ($\nu N$) cross section is known across neutrino energies from a few hundred MeV to a few PeV. Soon, ultra-high-energy (UHE) cosmic neutrinos, with energies above 100 PeV, could take us farther. So far, they have evaded discovery, but upcoming UHE neutrino telescopes endeavor to find them. We present the first detailed measurement forecasts of the UHE $\nu N$ cross section, geared to IceCube-Gen2, one of the leading detectors under planning. We use state-of-the-art ingredients in every stage of our forecasts: in the UHE neutrino flux predictions, the neutrino propagation inside Earth, the emission of neutrino-induced radio signals in the detector, their propagation and detection, and the treatment of backgrounds. After 10 years, if at least a few tens of UHE neutrino-induced events are detected, IceCube-Gen2 could measure the $\nu N$ cross section at center-of-mass energies of $\sqrt{s} \approx 10-100$ TeV for the first time, with a precision comparable to that of its theory prediction.


I. INTRODUCTION
Neutrino interactions with matter are powerful probes of particle physics: they map the deep structure of nuclei and nucleons, and may unearth evidence of new physics. Broadly stated, the higher the energy of the interacting neutrino, the sharper its probing power. Today, highenergy neutrino-matter interactions-in the form of the neutrino-nucleon (νN ) cross section-are known experimentally up to PeV neutrino energies, the highest detected so far. Yet, a trove of further insight likely lies in the measurement of the νN cross section at higher energies. Presently, those energies are practically out of the reach of existing detectors, but this limitation will likely be overcome in the coming years. Figure 1 shows the current landscape of measurements of the νN cross section. At neutrino energies from about 100 MeV to 350 GeV, the cross section is measured precisely using artificial neutrino beams from particle accelerators [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Soon, the planned accelerator-neutrino experiment FASERν [19] will reach TeV-scale energies, but not more. Beyond TeV neutrino energies, there is no existing or planned artificial neutrino beam. Thus, TeV-PeV cross-section measurements [20][21][22] used instead the high-energy cosmic neutrinos discovered by the IceCube neutrino telescope [23][24][25][26][27][28][29]. These are the most energetic neutrinos known so far, though not the most energetic neutrinos predicted. At even higher energies, of 100 PeV and above, the existence of ultra-high-energy (UHE) cosmic neutrinos was firmly predicted more than fifty years * vvalera@nbi.ku.dk † mbustamante@nbi.ku.dk ‡ christian.glaser@physics.uu.se ago [30,31]. They represent the only feasible way to extend νN cross section measurements to higher energies. Yet, because their predicted flux is low, they have so far evaded discovery [32,33]. Fortunately, a host of new neutrino telescopes [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49], designed to discover UHE neutrinos even if their flux is low in the next 10-20 years, may provide a way forward. For astrophysics, discovering UHE neutrinos would bring critical progress in understanding the long-standing origin of ultra-high-energy cosmic rays [49,50]. For particle physics, discovering UHE neutrinos, in general, would allow access to tests of fundamental physics in a new energy regime and, in particular, would allow us to further cross-section measurements [48,49,51,52].
However, detailed and realistic predictions for the capability of upcoming neutrino telescopes to measure the νN cross section, considering their design elements, are still lacking; see, however, Refs. [53] for important first estimates. To address this, and in order to capitalize on this upcoming opportunity, we present state-ofthe-art forecasts for the measurement of the UHE νN cross section. We gear our forecasts to the promising case of radio-detection of UHE neutrinos in the planned IceCube-Gen2 neutrino telescope [39], one of the leading next-generation detectors under design.
To make our forecasts realistic, we use state-of-theart ingredients in every stage of the calculation (see Section II). We model the propagation of UHE neutrinos inside the Earth and their radio-detection in detail. For the latter, we estimate the radio-detection response of the detector, via dedicated simulations of in-medium shower development and radio emission in IceCube-Gen2, and its energy and directional resolution. To capture the large uncertainty that exists in the prediction of UHE neutrinos, our cross-section forecasts factor in the uncertainty Accel. ν (projected) High-energy cosmic ν Ultra-high-energy cosmic ν (this work) IceCube-Gen2 Radio (10 yr, projected) FIG. 1. Neutrino-nucleon (νN ) charged-current (CC) cross section, measurements and predictions. Sub-TeV measurements are from accelerator-neutrino experiments [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Few-TeV measurements will be covered by the upcoming FASERν acceleratorneutrino experiment [19]. TeV-PeV measurements use high-energy astrophysical neutrinos detected by the IceCube neutrino telescope [20][21][22]. We forecast cross-section measurements above 100 PeV-the ultra-high-energy (UHE) range-where no measurement exists presently. Our forecasts are for radio-detection of UHE neutrinos in the planned IceCube-Gen2 [39], and are based on sophisticated simulations of UHE neutrino propagation inside Earth (Section V) and detector response (Section VI).
In this plot, we showcase 10-year forecasts for three attractive flux scenarios that may afford the most precise cross-section measurements [54,55] (models 2, 4, and 7 below). Table II shows results for many more flux models [29,[54][55][56][57][58][59][60][61][62] (Section IV). The precision of cross-section measurements that we report accounts for the significant uncertainty on the flux normalization (Section VII), and detector resolution (Section VI). For comparison, we show the state-of-the-art BGR18 calculation of the νN deep-inelastic-scattering cross section [63] (Section III A) on isoscalar matter, averaged between neutrinos and anti-neutrinos.
The goal of our forecasts is double. On the one hand, they are intended to showcase the reach of IceCube-Gen2 to make particle-physics measurements in a new energy regime, in as realistic a way as it is presently possible. On the other hand, and more generally, our forecasts are intended to stimulate the development of the particlephysics research programs of upcoming high-energy neutrino telescopes. The calculation framework that we introduce as part of our analysis can be adapted to other neutrino telescopes, and other measurement goals. Because the design of telescopes that will run in 10-20 years is being decided upon presently, our forecasts are timely. This paper is organized as follows. Section II showcases the salient points and strengths of our analysis. Section III presents a brief introduction to neutrino-nucleon deep inelastic scattering and outlines the strategy that we use to measure the cross section. Section IV introduces the various benchmark models of the cosmic UHE neutrino flux that we adopt. Section V describes the effect of neutrino propagation through Earth and how we compute it. Section VI gives an overview of radiodetection of UHE neutrinos, introduces the response of IceCube-Gen2, and shows how we estimate event rates in it. Section VII introduces the statistical analysis that we use to forecast cross-section measurements. Section VIII shows our resulting forecasts for the different benchmark flux models, and the effect on them of changing detector parameters. Section IX points out potential future research directions. Section X summarizes and concludes. Figure 2 illustrates the flow of calculations involved in producing our forecasts of the measurement of the UHE νN cross section. For a given prediction of the diffuse UHE neutrino flux, we propagate it through the Earth, where neutrinos interact with matter, and model its detection in the radio component of IceCube-Gen2. To make our forecasts realistic, timely, and representative of current unknowns, we use state-of-the-art ingredients at every stage of the calculation. In brief, these are:

II. SYNOPSIS AND CONTEXT
UHE νN cross section: UHE neutrinos interact with matter while propagating through the Earth, and when detected. The sensitivity to the νN cross section stems from an interplay of both effects; see Section III C 2. In our forecasts, we use a stateof-the-art prediction of the UHE cross section [63], and deviations from it, to compute neutrino propagation and detection. See Section III for details.
Diffuse neutrino flux: We have not discovered a flux of UHE neutrinos yet, so we are forced to consider different possibilities. In our forecasts, we assume a wide variety of UHE neutrino flux predictions that span the breadth available in the literature, from pessimistic to optimistic [29,[54][55][56][57][58][59][60][61][62]. We model in detail the flavor content of neutrinos and antineutrinos in the flux, and keep track of it throughout. See Section IV for details.
Neutrino propagation through Earth: When propagating UHE neutrinos through the Earth, we compute how neutrino interactions with matter modify the neutrino flux that reaches the detector. The modifications are energy-, direction-, and flavor-dependent, and are slightly different for neutrinos and anti-neutrinos. In our forecasts, we account for the dominant neutrino interaction-νN deep inelastic scattering-and for other interactions that are collectively non-negligible, via NuPropEarth [78,79]. See Section V for details.
Neutrino detection: We gear our forecasts to the detection of UHE neutrinos in the radio component of IceCube-Gen2, optimized for neutrino detection above 10 7 GeV. In our forecasts, we model the detector geometry, simulate the development of particle showers in the ice, the emission and propagation of radio signals from them, the detector response, including the direction-dependent response of the different antenna types in the array, via NuRa-dioMC [80] and NuRadioReco [81], the calculation of the trigger condition, and the uncertainties in reconstructing the energy and direction of detected events. See Section VI for details.
Non-neutrino backgrounds: In the radio-detection of UHE neutrinos, the main backgrounds that may mimic neutrinos are due to showers induced by atmospheric muons in the ice [82,83] and to showers induced in the atmosphere, mainly by cosmic rays, that penetrate into the ice. In our forecasts, we model the detection of atmospheric muons and show its effect on our ability to measure the νN cross section. We comment on the tentative effect of showers induced by cosmic rays, whose importance in neutrino radio-detection is still unclear.
Below, we expand on each of the above elements. Previous works studied the sensitivity of upcoming neutrino telescopes to the UHE νN cross section [74], or their capability to measure it [53,84]. They incorporated some of the above elements, often partially or in less detail; none incorporated all of them, or in full detail. Our analysis is the first detailed, realistic forecast of the capability to measure the UHE νN cross section. It has the flexibility needed to explore how sensitive the measurements are to changes in detector features. It is geared at IceCube-Gen2, but serves as a template for other neutrino telescopes under planning [48,49,61,85,86].

III. NEUTRINO-NUCLEON DEEP INELASTIC SCATTERING
A. The neutrino-nucleon DIS cross section At neutrino energies above 1 TeV, the neutrino-nucleon (νN ) cross section is dominated by deep inelastic scattering (DIS) [87][88][89][90], where the interacting neutrino scatters off the partons-i.e., quarks and gluons-that make up the nucleon. In the process, the nucleon, N , is broken up, and the final-state parton promptly hadronizes into hadrons, X. In charged-current (CC) DIS, mediated by a W boson, the final state contains in addition a charged lepton of the same flavor as the incoming neutrino, i.e., ν α + N → α − + X (α = e, µ, τ ). In neutral-current (NC) DIS, mediated by a Z boson, the final state contains instead a neutrino, i.e., ν α + N → ν α + X. Anti-neutrinos undergo the same DIS interactions, charge-conjugated.
In a νN DIS interaction, a fraction y-the inelasticity-of the initial neutrino energy, E ν , is transferred to X; the remainder fraction (1 − y) is transferred to the final-state lepton. The distribution of inelasticity The νN DIS cross sections are computed on the basis of the parton distribution functions (PDFs) inside the nucleon, i.e., the probability densities of finding valence and sea quarks of different flavors, and gluons, inside the nucleon [87][88][89][90]. They depend on two kinematic variables: the four-momentum transferred to the interacting parton, Q 2 , and the Bjorken scaling parameter x, the fraction of the nucleon moment carried by the interacting parton. PDFs are measured predominantly in charged lepton-nucleon scattering [93,94]; because they describe the nucleon, and not the lepton that probes it, they apply also to neutrino-nucleon scattering.
To compute νN DIS cross sections at an energy E ν , we need to evaluate PDFs at roughly x m W /E ν , where m W is the mass of the W boson. Presently, PDFs are known down to x ∼ 10 −5 . At TeV-PeV energies, this is sufficient to compute the cross section with low uncertainty. At EeV energies, relevant to our work, PDFs must be extrapolated to x ∼ 10 −7 ; this extrapolation is the main source of uncertainty in the calculation of EeV cross sections. Broadly stated, different calculations of the high-energy νN DIS cross section [63,72,74,91,[95][96][97][98][99][100] differ in three aspects: the perturbative order to which the cross section is computed, the PDF set that they use, and the procedure they use to extrapolate PDFs to low values of x. Competing calculations are close at TeV-PeV, but diverge at EeV; see, e.g., Fig. 3 in Ref. [21].
In our analysis, we adopt the state-of-the-art BGR18 νN DIS cross section calculation from Ref. [63], computed to next-to-next-to-leading order, as our baseline. It uses the recent NNPDF3.1sx PDFs [101], informed by D-meson data from LHCb [102][103][104], including the effect of nuclear corrections and heavy quarks, and treats consistently the behavior of the PDFs at small values of Bjorken-x. The uncertainty in the BGR18 cross section calculation ranges from 3% up to 100 PeV, to 10% at 10 EeV; see Fig. 1  total cross section used in our analysis, during the propagation of neutrinos inside the Earth using NuPropEarth [78,79] and in the computation of event rates at IceCube-Gen2. See Section III A for details. This is the central value of the BGR18 calculation [63], extracted from the HEDIS [78] module of GENIE [92]. The shaded region is the approximate energy window where the radio array of IceCube-Gen2 will be sensitive. Top: Cross sections for νe and νµ; they are equal. Bottom: Ratio of the ντ to νe CC cross sections showing the low-energy suppression of the former.
certainties in the PDFs [78]. Later, in Section VIII, we find that for the most optimistic UHE neutrino flux predictions, we might be able to measure the cross section to within comparable uncertainty; see Fig. 1. (The theory uncertainty in the cross section can be larger in the presence of nuclear corrections [63], which we ignore for Fig. 1, but comment on later.) Figure 3, and also Fig. 1, shows the central value of the BGR18 CC νN DIS cross section (in the notation introduced later, in Section VII A, this is σ std ). The cross section grows with E ν , and the softer-than-linear dependence on E ν due to the W mass is apparent from a few TeV on. At EeV energies, the cross section grows roughly as ∝ E 0.3 ν [91]. Within the energy window where IceCube-Gen2 will be sensitive, the cross sections are flavor-universal, yet, at low energies, Fig. 3 shows that the ν τ CC cross section is kinematically suppressed due to the large mass of the tauon. Figure 4 shows the corresponding inelasticity distributions, for NC and CC interactions. Later, in Section VI D, we use them to compute the expected event rate at IceCube-Gen2. The distributions peak at y = 0, but they are broad, i.e., there is a large spread in how the incoming neutrino energy in a νN DIS is split into CC, E ν = 10 7 GeV NC, E ν = 10 7 GeV CC, E ν = 10 10 GeV NC, E ν = 10 10 GeV FIG. 4. Inelasticity distributions in neutrino-proton (νp) DIS for the BGR18 cross section [63], for two illustrative neutrino energies. We use the inelasticity distributions to compute the expected shower rate at IceCube-Gen2 in Section VI D. At the energies relevant to our analysis, the distributions are equal for να andνα of all flavors. the final-state particles. This, in turn, generates a large spread in the energies of the neutrino-initiated showers at the detector, and on their detectability.
To make our forecasts self-consistent, we use the same BGR18 cross sections to compute the propagation of neutrinos inside Earth and their detection. At the neutrino energies relevant to us, E ν 10 7 GeV, the CC cross section is roughly twice the NC cross section, and differences between the cross sections of ν α andν α and between different flavors are small. Below 100 TeV, differences are larger, but those energies do not come into play in our analysis. We always use the CC and NC cross section of each flavor of ν α andν α individually. Below, as part of our forecasts, we use the central value of the BGR18 cross section as a baseline (σ std ), and also versions of it shifted up (> σ std ) and down (< σ std ) by constant factors.

B. Other neutrino interactions
At ultra-high energies, νN DIS is the dominant neutrino interaction. However, other sub-dominant interactions also affect the propagation of neutrinos inside the Earth. During propagation, we account for the following interactions, as implemented in NuPropEarth [78,79]; we defer to Ref. [78] for details: • νN DIS on the partons of the nucleon: The dominant interaction channel at ultra-high ener-gies; see Section III A.
• Neutrino DIS on the photon field of the nucleon: The neutrino interacts with a lepton generated by the photon field of the nucleon. This is negligible except when the neutrino can produce an on-shell W boson that enhances the cross section resonantly, i.e., when √ 2m N E ν m W , or E ν 3 TeV. It can account for a correction of up to 3% of the total DIS cross section [105][106][107][108][109].
• Coherent neutrino-nucleus scattering: The neutrino interacts coherently with the photon field of the target nucleus [110][111][112]. The cross section is ∝ Z 2 , where Z is the atomic number of the nucleus; thus, it matters mostly for heavy nuclei. This process is important only in scatterings with small transferred momentum, Q 1 GeV, where it can contribute up to 10% of the total cross section.
• Elastic and diffractive neutrino scattering on nucleons: The neutrino interacts with the photon field of individual nucleons. This process is important in scattering with small Q ∼ GeV, where the elastic component of the form factors of the proton become relevant. When E ν 3 TeV, the process may create an on-shell W boson resonantly, contributing to the resonant cross section of the neutrino DIS on the photon field of the nucleon.
• Neutrino scattering on atomic electrons: Neutrino scattering on atomic electrons is negligible except for high-energyν e . When the center-ofmass energy is √ 2m e E ν ≈ m W , i.e., when E ν ≈ 6.3 PeV, theν e e scattering is resonance and produces an on-shell W ; this is known as the Glashow resonance [113,114]. Around the resonance energy, theν e e cross section dominates; it is roughly 200 larger than the νN DIS cross section. We adopt the Glashow-resonance cross section from Ref. [107], computed to next-to-leading order.
The sub-dominant interactions increase the attenuation of the UHE neutrino fluxes by up to 10%, when compared to νN DIS only, especially for neutrinos coming into the detector from around the horizon [78]. Below, when computing the propagation of neutrinos through the Earth and their resulting fluxes at the detector, we always account for all of the above interactions; see Section V. In all of the interactions, because of the high energies, finalstate neutrinos are nearly co-linear with initial-state neutrinos; any transverse momentum is negligible compared to the forward momentum, and we ignore it.
C. High-energy νN DIS using cosmic neutrinos

Motivation
Measuring the high-energy νN DIS cross section offers the possibility to probe fundamental physics on two fronts. First, the higher the energy of the neutrino, the smaller the value of x that it probes; see Section III A. This allows us to probe the structure of nucleons deeper, testing predictions of potentially non-linear behavior of the PDFs, such as from BFKL theory [64][65][66][67] and color-glass condensates [68], and to look for electroweak sphalerons [69]; see, e.g., Refs. [69][70][71][72]. Second, the higher the neutrino energy, the higher the energy scale probed where new physics could affect the cross section. Possibilities include, e.g., leptoquarks, extra dimensions, and new gauge bosons [53,[73][74][75][76][77].
Existing measurements of the νN DIS cross sections using accelerator neutrinos reach 350 GeV [12]; see Fig. 1. Upcoming accelerator-neutrino experiments, like FASERν [19], should measure the cross section up to a few TeV. To measure it at higher energies, we must use neutrinos from natural sources: atmospheric neutrinos, from tens of TeV to roughly 100 TeV; high-energy cosmic neutrinos from 10 TeV to 10 PeV; and ultra-highenergy neutrinos, from 100 PeV on. Measurements using the former two exist [20][21][22]; we forecast measurements using the latter, where no measurement exists yet.

Overview
We base our forecasts of cross-section measurements on estimates of the expected number of detected neutrinoinduced events in IceCube-Gen2. To understand where the sensitivity to the cross section comes from, we estimate the number of detected neutrinos of energy E ν coming into the detector from zenith angle θ z as Here, Φ ν is the diffuse, isotropic flux of UHE cosmic neutrinos that arrives at the surface of the Earth, σ is the νN cross section, L is distance from the surface of the Earth to the detector, L νN ≡ (σn N ) −1 is the mean free path, and n N is the average number density of nucleons encountered by the neutrino along its way inside Earth. The latter depends on θ z and on the internal matter density of Earth. Equation (1) is merely a simplified expression for the purpose of providing insight, and is not used to produce our results. The full treatment of neutrino propagation and detection that we use to produce results is in Sections V and VI, respectively. Equation (1) shows that the number of events depends on the νN cross section doubly. During propagation, the cross section acts via the exponential, attenuating the flux of neutrinos as they go through Earth (in the full calculation, there is also regeneration of lower-energy neutrinos, which we ignore momentarily). Higher neutrino energies-and, therefore, larger cross section-and longer distances traveled inside Earth lead to stronger attenuation. At detection, the cross section acts proportionally; the larger it is, the higher the chances of detecting a neutrino that arrives at the detector. The sensitivity of neutrino telescopes to the high-energy νN DIS cross section stems from the interplay of these two effects.
We extract the cross section by examining the angular distribution of detected events [20-22, 74, 115-118]. For events induced by neutrinos arriving from above the detector, i.e., downgoing events, where L is small, the attenuation is negligible. In this case, the right-hand side of Eq. (1) becomes ∝ Φ ν σ and, therefore, the sensitivity to the cross section is mild. This is due to the degeneracy between Φ ν and σ: for a fixed event rate N ν , a higher flux can be traded off for a lower cross section, and vice versa. For events induced by neutrinos arriving from well below the horizon, i.e., upgoing events, where L is large, the attenuation is strong. In this case, the event rate is low: the higher the neutrino energy-i.e., the higher the cross section-the stronger the attenuation.
For events induced by neutrinos arriving horizontally and nearly horizontally into the detector, the two effects above balance out. Neutrinos from this direction travel tens to hundreds of kilometers underground, enough for the flux to be attenuated, but not eliminated. The higher the neutrino energy, the narrower the solid angle around the horizon from where neutrinos arrive to the detector. At ultra-high energies, neutrinos can only arrive at the detector from a few degrees around the horizon, where the distance traveled underground is not too long (see, e.g., Fig. A2 in Ref. [21]); these are called Earthskimming neutrinos [119].
Combining events from all directions breaks the degeneracy between Φ ν and σ in downgoing events, and, via the attenuation of upgoing and near-horizontal events, grants us sensitivity to σ [74,[115][116][117][118]. Yet, to make use of this, the detector requires a good angular resolution around the horizon. This is necessary to infer the direction of the incoming neutrino precisely, and, therefore, the column density of matter that it traversed on its way to the detector. In IceCube, at TeV-PeV energies, the angular resolution varies from sub-degree for ν µ -induced track events, to tens of degrees for shower events induced by all neutrino flavors. In our forecasts, geared at radiodetection of EeV-scale neutrinos in IceCube-Gen2, we adopt a baseline resolution of 2 • to produce our main results, and also explore the effect of alternative choices; see Section VI and Fig. 21 for details. Figure 1 shows the three existing measurements of the TeV-PeV νN DIS cross section, all based on IceCube data [20][21][22]. They used different analysis strategies and data sets; they agree with Standard Model (SM) predictions, though the measurement uncertainties are large.

Existing measurements: TeV-PeV
In addition to the TeV-PeV measurements, there are complementary studies to measure the cross section in IceCube from 100 GeV to a few TeV [120].
Reference [20] used roughly 10800 through-going neutrino-induced muon tracks, predominantly of atmospheric origin, to measure the cross section in a single neutrino energy bin spanning 6.30-980 TeV. In a through-going track event, a ν µ interacts at an unknown position outside the detector and creates a high-energy muon that crosses part of it. This analysis let the normalization of the CC and NC cross sections float, fit it to the data, and found the cross section to be 1.30 +0. 30 −0.26 times the SM prediction from Ref. [121]. The error is approximately equal parts statistical and systematic; the latter is mainly due to the difficulty in reconstructing the neutrino energy from the measured muon energy. Improvements on this measurement strategy are ongoing [122].
References [21,22] used instead High Energy Starting Events (HESE), predominantly of cosmic origin, to measure the cross section from 20 TeV to a few PeV. Unlike a through-going track, in a HESE event the neutrino interacts inside the detector. The ensuing shower deposits a large fraction of its energy in the instrumented volume. As a result, the neutrino energy is reconstructed accurately, which facilitates measuring the cross section in multiple energy bins. However, because HESE events are relatively rare, these analyses are limited by low event rates. The CC interaction of ν e or ν τ , or the NC interaction of a neutrino of any flavor, induces a shower that is typically contained in the detector. The NC interaction of ν µ induces in addition a muon track that starts inside the detector, but ranges out of it. Reference [21] used 58 HESE showers collected in 6 years to measure the cross section in bins in the range 18 TeV-2 PeV, with an accuracy of roughly half an order of magnitude in each bin. Reference [22] used 60 HESE showers and tracks collected in 7.5 years to measure the cross section in bins in the range 60 TeV-10 PeV with increased accuracy, thanks to reduced detector systematics. A related analysis [123] used contained events to make the first measurement of the multi-TeV inelasticity distribution.

This work: forecasts at EeV
To produce forecasts of UHE cross-section measurements in IceCube-Gen2, we adopt the BGR18 cross section as the baseline, and explore measurement prospects for a variety of diffuse UHE flux models. Our procedure is reminiscent of analyses that use HESE events [20,21]. Section VII describes it in detail. In our forecasts, variations in the cross section affect the neutrino propagation inside Earth and their detection. To account for the degeneracy between the flux and the cross section, we always forecast measuring both simultaneously. We adopt a Bayesian approach in our statistical analysis, and account for statistical fluctuations in the number of neutrino-induced events and non-neutrino backgrounds.
IceCube HESE (7.5 yr) extrapolated IceCube-Gen2 Radio energy range A u g e r u .l .  [29,[54][55][56][57][58][59][60][61][62]. Highlighted models 2, 4, 6, and 7 receive special attention in the main text, but results for all models are shown in Tables I and II, and in Appendix A. The upper limits on the UHE neutrino flux are from IceCube [32] and the Pierre Auger Observatory [33]. The projected sensitivity of the radio component of IceCube-Gen2 is from Ref. [43]. See Section IV for an overview of the flux models.

IV. ULTRA-HIGH-ENERGY NEUTRINOS
Ultra-high-energy neutrinos, with energies of 100 PeV and above, are expected to come from interactions of ultra-high-energy cosmic rays (UHECRs), of EeV-scale energies, on radiation or matter. UHECR interactions may occur inside the cosmic accelerators that are their sources or outside them, en route to Earth. In the former case, the resulting UHE neutrinos are dubbed source neutrinos (or astrophysical neutrinos); in the latter, they are dubbed cosmogenic neutrinos. Because of unknowns in the properties of UHECRs and their sources, there is a large spread in the predicted neutrino flux normalization and the shape of the neutrino energy spectrum. In our analysis, we consider a wide breadth of benchmark flux models from the literature in order to represent this spread. Below, we introduce our benchmark flux models, and the choices that we make in building them.

A. Overview
Cosmic accelerators are expected to generate a population of non-thermal UHECRs with a power-law energy spectrum [124]. The interaction of UHECR protons on matter (pp) or radiation (pγ) often creates a short-lived ∆(1232) resonance that promptly decays into charged pions. Upon decaying, they create neutrinos, i.e., π + → µ + + ν µ , followed by µ + → e + + ν e +ν µ , and their charge-conjugated processes. Each final-state neutrino receives 5% of the energy of the parent proton, on average. En route to Earth, neutrino oscillations change the flavor composition of the flux, i.e., the proportion of ν e , ν µ , and ν τ in it; we account for this in Section IV B.
For UHE neutrinos produced in pp interactions, the resulting neutrino energy spectrum follows the power-law energy spectrum of the parent protons, and may extend unbroken down to low energies [58]. The spectral index of the neutrino spectrum is inherited from the spectral index of the parent proton spectrum.
For UHE neutrinos produced in photohadronic, i.e., pγ, interactions, the resulting neutrino energy spectrum is determined by the energy spectra of the interacting protons and photons. Because the proton spectrum is a power law and the photon spectrum is peaked or is a power law with a spectral break, the resulting neutrino spectrum is peaked. The neutrino spectrum peaks at an energy determined by the ∆ resonance energy; the width of the neutrino peak is determined by the widths of the photon and proton spectra.
Cosmogenic neutrinos, or GZK (Greisen-Zatsepin-Kuzmin) neutrinos, were first predicted in the late 1960s [30,31,125]. They are expected to be made during the extragalactic propagation of UHECRs, in photohadronic interactions on the cosmic microwave background (CMB), for neutrinos of energies typically in the EeV-scale, and on the extragalactic background light (EBL), for neutrino of energies typically in the tens of PeV. (Cosmogenic anti-neutrinos have an additional contribution from the beta-decay of neutrons produced in photohadronic interactions, typically around PeV energies, outside the region of interest for neutrino radiodetection in IceCube-Gen2.) Because the CMB photon spectrum is well-known, the uncertainty in the prediction of the cosmogenic neutrino flux comes mainly from uncertainties in the energy spectrum, mass composition, and maximum energies of UHECRs, as measured by the Pierre Auger Observatory [126][127][128] and the Telescope Array (TA) [129,130], and in the abundance of the UHECR sources at different redshifts. See, e.g., Refs. [49,131] for an overview. Generally, a harder UHECR energy spectrum, lighter mass composition, higher maximum energy, and a source number density that peaks at intermediate redshifts lead to a higher cosmogenic neutrino flux; see, e.g., Refs. [59,[132][133][134].
UHE source neutrinos are expected to be made in either pp interactions, pγ interactions, or both, inside UHECR sources. When photohadronic interactions are dominant, the spectrum of UHE source neutrinos has a similar shape to that of cosmogenic neutrinos, except that it contains a single pγ bump, since there is typically a single relevant target spectrum of photons inside the sources. (UHE source anti-neutrinos also have an additional contribution from the beta-decay of neutrons produced in photohadronic interactions, typically at energies too low to be relevant for our analysis.) In some models of UHE neutrino production in cosmic-ray reservoirs [58], the contribution of neutrinos from pp interactions extends to low energies, and the contribution of pγ interactions is dominant at high energies.
In realistic models of high-energy neutrino production, including some of the ones that we consider in our analysis, different neutrino production channels become accessible at different energies. In pγ interactions [135][136][137], neutrino production occurs dominantly via the ∆(1232) resonance at the lowest energies, with a sub-leading contribution from direct (t-channel) production, via heavier resonances at intermediate energies, and via multipion production at the highest energies. In pp interactions [138], the neutrino yield evolves with energy as a result of the evolving pion multiplicity. Moreover, the physical conditions in the production region affect the energy of charged particles-protons, muons, pions, and kaons-whose decay yields neutrinos. For instance, intense magnetic fields may induce important synchrotron energy losses [139][140][141][142] that cap the high-energy neutrino yield, while re-acceleration of charged particles might counteract these losses [143]. Further, the presence of nuclei heavier than protons, and the nuclear cascades initiated by their interactions with source environments, introduce additional nuance [144,145].

B. Flavor and ν vs.ν composition in our analysis
Because, at different energies, different neutrino production channels dominate (see Section IV A) and the physical conditions at the sources affect charged particles differently, the flavor composition of the UHE cosmic neutrinos, i.e., the proportion of ν e +ν e , ν µ +ν µ , and ν τ +ν τ in the total flux, and the ratio of ν α toν α , are expected to evolve with neutrino energy. This matters for the purpose of propagating neutrinos through the Earth, on their way to the detector, and of forecasting their detection rates. Section V shows that neutrinos of different flavor are affected differently by their passage through Earth. Differences between ν α vs.ν α are small, though we keep track of them. The exception where differences are significant is the case of ν e vs.ν e , since onlȳ ν e interact via the Glashow resonance [113]. Section VI shows that neutrinos of different flavors have different interaction rates and deposit energy differently at IceCube-Gen2; there, differences between ν α vs.ν α are also small, though we keep track of them.
In Section IV C, we introduce the benchmark UHE neutrino flux models that we later use to forecast νN cross-section measurements in IceCube-Gen2. In order to make our predictions as informed as possible, we model the flavor composition and ν α vs.ν α content of the benchmark fluxes as accurately as possible. In doing so, neutrino flavor transitions play a key role. Below, we explain how we compute them.
Formally, the probability P αβ of a flavor transition ν α → ν β oscillates as a function of the distance traveled by the neutrino. However, for high-energy cosmic neutrinos, the oscillation length, which is ∝ 1/E ν , is tiny compared to the typical traveled distance of Mpc-Gpc, so the probability oscillates rapidly. In addition, neutrino telescopes have limited energy resolution [147]. As a consequence, in practice, oscillations average out, and we are sensitive only to the average probability [148], Below, we compute the flux of each neutrino flavor at Earth for our benchmark flux models by evaluating P βα using values of the mixing parameters forecast for the 2030s [149], anchored on present measurements [150,151]. For this purpose, our benchmark flux models fall into three categories, depending on what information is available to us to build the model with. For each, we compute the flux of each neutrino species at Earth differently: (a) Flux models for which we have available the preoscillation flux of each neutrino species separately (models 3-7 below). In this case, we compute the flux of ν α at Earth, Φ να , from the pre-oscillation fluxes that we have available, Φ ν β ,S , as and similarly for the flux Φν α ofν α , but changing Φ ν β ,S → Φν β ,S . Because in all of our benchmark flux models neutrinos are produced by pion, kaon, and neutron decays, only ν e ,ν e , ν µ , andν µ exist preoscillation; however, after oscillations, all six species are populated in the flux at Earth.
(b) Flux models for which we only have available the sum of the oscillated fluxes of all neutrino species at Earth (models 1, 8-11 below). In this case, we consider the flavor composition to be energy-independent and split the flux of each flavor evenly between ν α toν α at all energies. The latter assumption holds approximately, but can have large deviations at high energy, depending on model-dependent details of the neutrino production; see, e.g., Ref. [136] and Fig. 6 below. To estimate the flavor composition, we assume that all neutrinos are made in pion decays, i.e., π + → µ + + ν µ , followed by µ + → e + + ν e +ν µ , and their charge-conjugated processes. Hence, pre-oscillation, the flavor composition is f π e,S , f π µ,S , f π τ,S ≡ 1 3 , 2 3 , 0 , where f π β,S is the ratio of ν β +ν β to the total. After oscillations, at Earth, the flavor ratios become Thus, starting, from the all-species oscillated flux at Earth that we have available, Φ 6ν , we estimate the oscillated fluxes of ν α andν α at Earth as where the factor of 1/2 splits the flux of ν α +ν α evenly between them.
(c) Flux models for which we only have available the ν µ +ν µ oscillated flux at Earth (model 2 below). Like with fluxes in category (b), we consider the flavor composition to be energy-independent and split the flux of each flavor evenly between ν α toν α at all energies. Starting from the flux of ν µ +ν µ at Earth that we have available, we estimate the oscillated fluxes of ν α andν α at Earth as where the factor of 1/2 splits the flux of ν α +ν α evenly between them.
For benchmark flux model 12, the flux of each neutrino species at Earth is directly available [152]; we adopt them without modification. In our analysis, we forecast measurements in IceCube-Gen2 in the 2030s. By then, the values of the mixing parameters are expected to be known more precisely than today [150,151], thanks to the upcoming oscillation experiments DUNE [153], Hyper-Kamiokande [154], and JUNO [155]. Assuming that the true values of the mixing parameters are equal to their present-day best-fit values from the NuFit 5.0 global fit to oscillation data [150,151], and that the neutrino mass ordering is normal, by 2030 we expect that [149] ignoring correlations between the mixing parameters.
Our selection of benchmark flux models is representative of the breadth of theoretical predictions available in the literature at the time of writing. The highest of our benchmark fluxes-models 4, 6, and 12-saturate the present upper limits on the UHE neutrino flux. The lowest-models 1, 3, and 5-lie below the 10-year differential sensitivity of IceCube-Gen2. The remaining flux models lie in-between these two extremes. Later, in Section VIII, we find that measuring the UHE νN cross section should be possible for all but the lowest flux models.
Below we present an overview of the benchmark flux models. We defer to their original publications for details (see also Ref. [156]).
2. IceCube ν µ (9.5 yr) extrapolated [55]: Using 6.5 × 10 5 through-going muon tracks collected over 9.5 yr, IceCube fit the ν µ +ν µ astrophysical neutrino flux using a power law, To build our benchmark flux model 3, we extend this flux to 10 10 GeV, unbroken, using the best-fit values of Φ νµ+νµ,0 and γ νµ+νµ . We assume that the flux of each neutrino species shares this common value of γ νµ+νµ , and estimate it using Eq. (6).    sources, distributed in redshift, and their flux is normalized by fitting the predicted UHECR energy spectrum and mass composition at Earth to recent data from the Pierre Auger Observatory [157,158]. References [133,134] predict similar fluxes, also from fits to Auger data. UHECR interactions on the CMB and EBL, including photodisintegration and photohadronic processes, are computed using PriNCe [159]. Because the best-fit value of the maximum UHECR rigidity is low, 2.5 × 10 9 GV, there are relatively few UHECRs at the highest energies, and so the neutrino yield is low. To build our benchmark model 3, we compute the pre-oscillation fluxes of ν e ,ν e , ν µ , andν µ as functions of energy using PriNCe, and use Eq. (3) to transform them into the oscillated fluxes of all species at Earth.

4.
Bergman & van Vliet, fit to TA UHECRs (cosmogenic) [61]: Cosmogenic neutrinos are generated in the same way as for the benchmark flux model 3, but instead fitting the UHECR energy spectrum and mass composition at Earth to recent data from TA [160,161]. Reference [162] predicts a similar flux. Because the TA data is compatible with a lighter UHECR mass composition and a higher maximum rigidity, the neutrino flux inferred using TA data is larger than with Auger data (model 3). To build our benchmark model 4, the preoscillation fluxes of ν e ,ν e , ν µ , andν µ as functions of energy are computed using CRPropa3 [163,164], and use Eq. (3) to transform them into the oscillated fluxes of all species at Earth.

5.
Rodrigues et al., all AGN (cosmogenic) [54]: Neutrinos are produced by the entire population of active galactic nuclei (AGN), which serve as UHECR accelerators. AGN are divided into three sub-populations: low-luminosity BL Lacs, high-luminosity BL Lacs, and flat-spectrum radio quasars (FSRQs). The number density of each sub-population evolves differently with redshift and luminosity [165,166]. Before escaping, UHECRs interact in the AGN jets, via photodisintegration and photohadronic processes [144,167]. Cosmogenic neutrinos are produced by the UHE-CRs that escape, and their flux is computed using PriNCe [159]. Source neutrinos are produced inside the jets, and their flux is computed using Neu-CosmA [136,137,168]. The predicted UHECR energy spectrum and mass composition at Earth agree with Auger data [157], while the UHE neutrino flux satisfies the IceCube upper limit [32]. Low-luminosity BL Lacs explain the UHECR flux, while FSRQs dominate neutrino production. To build our benchmark model 5, we adopt the maximum allowed cosmogenic neutrino flux from the entire population of AGN ( Fig. 2 in Ref. [54]). We take the pre-oscillation fluxes of ν e ,ν e , ν µ , and ν µ as functions of energy [169], and use Eq. To build our benchmark model 6, we adopt the maximum allowed source neutrino flux from the entire population of AGN ( Fig. 2 in Ref. [54]). We take the preoscillation fluxes of ν e ,ν e , ν µ , andν µ as functions of energy [169], and use Eq. (3) to transform them into the oscillated fluxes of all species at Earth.

Rodrigues et al., HL BL Lacs (cosmogenic) [54]:
UHECRs and neutrinos are produced only by highluminosity (HL) BL Lacs. The predicted UHECRs agree with the Auger energy spectrum above the ankle, but are lighter than the Auger mass composition above a few EeV. We adopt the cosmogenic neutrino spectrum from HL BL Lacs ( Fig. 5 in Ref. [54]) as benchmark because it peaks at energies higher than the benchmark models 5 and 6, and has a normalization in-between theirs. We take the pre-oscillation fluxes of ν e ,ν e , ν µ , andν µ as functions of energy [169], and use Eq. (3) to transform them into oscillated fluxes of all species at Earth.

8.
Fang & Murase, cosmic-ray reservoirs (cosmogenic + source) [58]: Neutrinos are produced in a grandunified multi-messenger model of high-energy emission where UHECRs are accelerated in the jets of supermassive black holes of radio-loud AGN embedded in galaxy clusters that act as cosmic-ray reservoirs. There, UHECRs remain confined for 1-10 Gyr and produce UHE neutrinos via pγ and pp interactions. From 100 TeV to 100 PeV, neutrinos are primarily made inside the clusters, in UHECR interactions on the intra-cluster medium; above 10 9 GeV, neutrinos are primarily cosmogenic. The neutrino flux normalization results from fitting the predicted UHECR energy spectrum and mass composition to Auger data [170], and the predicted TeV-PeV neutrino flux to Ice-Cube data [27,171]. Reference [58] provided the all-species neutrino flux, Φ 6ν . To build our benchmark flux model 8, we use it to estimate the flux of each neutrino species using Eq. (5).
9. Fang et al., newborn pulsars (source) [56]: Fastspinning newborn pulsars that harbor intense surface magnetic fields, of up to 10 13 G, may efficiently accelerate charged particles in the pulsar wind during pulsar spin-down. Accelerated particles propagate through the expanding supernova ejecta that surrounds the pulsar; as they do, pp interactions on the ejecta produce neutrinos. Reference [56] computed the diffuse flux of neutrinos produced by the cosmological population of newborn pulsars, integrated over their neutrino-producing lifetimes, with a spread in magnetic field intensity and spin period, and distributed in redshift following the star formation rate (SFR). (We consider only the neutrino contribution from the sources, not the contribution of cosmogenic neutrinos produced by UHE-CRs emitted by the pulsars, which is of the same order [56].) Reference [56] provided the all-species neutrino flux, Φ 6ν . To build our benchmark flux model 9, we use it to estimate the flux of each neutrino species using Eq. (5).
10. Padovani et al., BL Lacs (source) [57]: The neutrino flux is obtained within the framework of the simplified view of blazars [172]. Neutrinos are produced in photohadronic interactions inside the jets of BL Lacs, whose population is simulated using a spread of redshifts and source features like synchrotron peak energies and X-ray flux. The flux prediction was originally constructed to explain the TeV-PeV neutrino range; we adopt it because it spills into the UHE regime. A key parameter of the model is Y νγ , the ratio of the neutrino intensity to the gamma-ray intensity. Following Ref. [173], to satisfy the IceCube upper limit on the UHE neutrino flux, we set Y νγ = 0.13. Reference [57] provided the all-species neutrino flux, Φ 6ν . To build our benchmark flux model 10, we use it to estimate the flux of each neutrino species using Eq. (5).
11. Muzio et al., maximum extra p component (cosmogenic + source) [60]: Cosmogenic and source neutrinos are produced via photohadronic interactions within the UFA15 multi-messenger framework [174], where sources emit UHECRs whose energy spectrum and mass composition at Earth are fit to Auger data. Reference [60] added a subdominant UHECR pure-proton component that escapes the sources with energies above 10 9 GeV, motivated in part by the observation, in Auger, of a slowdown in the increase of average nuclear mass with energy [126], and that enhances the neutrino flux. We adopt the maximum allowed neutrino flux that results from the joint single-mass UFA15 plus pure-proton components, using Sybill2.3c [175] for the hadronic interaction of UHECRs in the atmosphere ( Fig. 9 in Ref. [60]). Reference [60] provided the all-species neutrino flux, Φ 6ν . To build our benchmark flux model 11, we use it to estimate the flux of each neutrino species using Eq. (5).
12. Muzio et al., fit to Auger & IceCube (cosmogenic + source) [62]: Cosmogenic and source neutrinos are produced via photohadronic interactions within the UFA15 multi-messenger framework (see above); in addition, source neutrinos are produced via pp interactions of UHECRs in the source environment. Neutrinos from hadronic interactions dominate at low energies, below the ∆ resonance energy. UHECR predictions are fit to Auger data. The total neutrino flux from Ref. [62] includes contributions from UHECR sources and non-UHECR sources; the total flux is fit to the IceCube TeV-PeV neutrino flux measurement [114,176]. We adopt the best-fit total neutrino flux ("UHECR ν" plus "Non-UHECR ν from Fig. 1 in Ref. [62]). To build our benchmark model 12, we use the oscillated fluxes of each neutrino species at Earth as a function of energy, which are available directly from the calculation [152]. Figure 6 shows the breakdown into the flux of each neutrino species for the benchmark models; see Section IV B. For models 1, 2, 8-11, the ratio of each species to the total flux is constant in energy. For models 3-7 and 12, for which non-trivial energy evolution of the flux of each species separately is available, common trends are evident. At low energies, typically below the energy range of the IceCube-Gen2 radio component,ν α dominate due to the presence and oscillation ofν e produced in the beta-decay of neutrons and neutron-rich isotopes created in UHECR interactions, mainly during their extragalactic propagation. At higher energies, neutrinos are produced by pion decay; throughout the IceCube-Gen2 energy range, flavor equipartition holds approximately. Roughly within the range 10 7 -10 9 GeV, there is a slight excess of neutrinos over anti-neutrinos, because more π + than π − are produced. At the highest energies, multipion production dominates, and the excess flips. Later, in Fig. 11, we show how the flavor composition is affected by neutrino propagation through the Earth.

V. NEUTRINO PROPAGATION INSIDE EARTH
To compute the expected rate of neutrino-initiated showers at a neutrino telescope, first we compute the flux of neutrinos that reaches it, after propagating through the Earth across different directions.

A. Computing neutrino propagation
Above TeV energies, the dominant interaction that neutrinos undergo while propagating inside the Earth is νN DIS, NC and CC; see Section III. NC scatterings pile up originally high-energy neutrinos at low energies, while CC scatterings remove neutrinos from the flux altogether. The exception is the CC scattering of ν τ , where the finalstate tauon may propagate some distance before decaying into a ν τ , albeit with an energy lower than that of the original ν τ ; this is known as "ν τ regeneration". Because of it, the flux of ν τ is less attenuated than that of other flavors. This becomes especially important above 10 PeV, where the final-state ν τ are still high-energy.
We propagate UHE neutrinos towards the detector along different directions, parametrized as cos θ z , where θ z is the zenith angle measured from the location of the detector. For us, this is the South Pole (cos θ z = 1), where IceCube-Gen2 will be located; see Fig. 2. We compute the energy spectra of each neutrino species that reach the detector from −1 ≤ cos θ z ≤ +1. Figure 7 shows examples of neutrinos hitting the surface of the simulated IceCube-Gen2 detector volume from different directions. We model the detector volume as a cylinder of radius 12.6 km and height 1.50 km, buried underground at a distance of 1.51 km from the surface of the Earth to the bottom of the cylinder. The top surface area of the cylinder is 500 km 2 [39]. Once a neutrino reaches the surface of the cylinder, we stop its propagation. Inside the cylinder, the propagation of the neutrino is computed separately, in the detection step of our calculation, where any further interaction that occurs within the detector volume initiates a particle shower that emits a radio signal that might be detected by underground antennas; see Section VI. Modeling the detector volume as a cylinder vs. as a point impacts by up to 10% the attenuation of neutrinos that reach it from directions around the horizon, i.e., Earth-skimming neutrinos, for which the detector is of comparable size to the distance traveled inside the Earth. This is especially relevant because these neutrinos offer the greatest sensitivity to the νN cross section; see Section III C.
Depending on the direction of the neutrino, it will encounter a different matter column density. To account for this, for the internal matter density of Earth, we adopt the Preliminary Reference Earth Model [177], built from seismographic data, which models the density radially out from the center of a spherical Earth, as a series of concentric layers of increasing density towards the center. For our calculations, since IceCube-Gen2 will be embedded in the Antarctic ice, we add a layer of ice of thickness 3 km at the surface of the Earth. In addition, the composition of matter inside the Earth changes with radial distance: deeper layers contain heavier elements-iron, nickel-than shallower layers. Further, matter is, in general, not isoscalar, though this affects mainly neutrino energies below 1 TeV [78]. When propagating neutrinos inside the Earth, we account for the changes in density and composition as a function of position inside Earth.
As an illustration only, and not accounting for ν τ regeneration, the exponential dampening in Eq. (1) describes the attenuation of the neutrino flux inside Earth. (We do not use those simplified expressions to produce our results, but more sophisticated methods; see Section VI D.) The attenuation due to CC interactions is stronger the higher the neutrino energy and the longer the length of the path traveled by the neutrinos inside the Earth. The low-energy pile-up due to NC interactions has a similar dependence on energy and direction As a result of the interactions inside the Earth, while the neutrino flux is isotropic at the surface of the Earth, it has become anisotropic by the time it reaches the detector. At EeV energies, no detectable flux reaches the detector from below; instead, the flux comes from above, where it is only lightly attenuated by the detector overburden, and from around the horizon, where the attenuation is significant to modify the shape of the spectrum, but not enough to eliminate it.
In NuPropEarth, for the νN DIS cross section, we Number of ν reaching IceCube-Gen2,  [181]. We have modified NuPro-pEarth and HEDIS to be able to use versions of the BGR18 cross section that are scaled up and down by a constant factor, as part of our method of measuring the cross section; see Section VII. For the cross sections of the sub-dominant neutrino interactions (see Section III B), we use their implementations in HEDIS, unmodified. Figure 8 shows energy histograms at the detector resulting from propagating a mono-energetic Earthskimming neutrino beam inside the Earth. For ν µ , the cascading down to lower energies due to multiple interactions inside the Earth is evident. For ν τ , the effect of regeneration is evident as a pile-up at lower energies. Figure 8 shows that changes in the νN DIS cross section affect the propagation inside Earth significantly. For ν µ , and also for ν e , not shown, the main effect of a higher cross section is to attenuate the flux further. For ν τ , a higher cross section shifts the peak of the pile-up to lower energies, due to a larger number of neutrino interactions.

B. Computational speed-ups
As part of our statistical analysis below (see Section VII D), we need to propagate many different UHE spectra of ν e ,ν e , ν µ ,ν µ , ν τ , andν τ separately through the Earth, for different values of the νN DIS cross section, along different directions, and with high accuracy. This is a computationally taxing task.
We circumvent this limitation as follows. First, we inject a large number N surf να = 10 6 of ν α at the surface of the Earth, each with initial energy E surf ν , and propagate it along different directions towards the detector, using NuPropEarth. Upon arriving at the detector, the neutrinos are no longer mono-energetic, but their final energies, E det ν , are spread out as a result of interactions inside Earth, i.e., N det . Second, we compute the transmission coefficients as and save T να as look-up tables. We do this separately for ν α andν α of all flavors. Third, given any UHE neutrino spectrum at the surface of the Earth, Φ να , we use the precomputed T να to estimate the average flux at the detector within an interval of final energy [E ν , E ν + ∆E ν ] as and similarly forν α . Finally, we approximate the true spectrum by its binned average, i.e., Φ det να (cos θ z , E ν ) ≈ Φ det να (cos θ z , E ν , ∆E ν ). We pre-compute the look-up coefficients T να in fine grids of E surf ν and E det ν . We repeat the above procedure to generate look-up tables for different values of the νN cross section, since we need them for our statistical procedure; see Section VII D. Figure 9 shows sample transmission coefficients T νµ and T ντ , for directions at and around the horizon. For ν µ , and also for ν e , not shown, because the νN cross section grows with energy, the cascading down to lower energies becomes more important the higher the injected energy E surf ν . It is most significant when neutrinos arrive from below the horizon, due to the larger matter column density that they traverse. For ν τ , in addition, the presence of regeneration is evident for E surf ν 10 PeV, for neutrinos coming from the horizon and below it. At the detector, regenerated ν τ are concentrated in a band centered around E det ν ≈ 10 PeV. When producing our results, we use Eq. (12) to compute neutrino spectra. By doing this, we circumvent the computationally intensive need to propagate every time each benchmark neutrino flux from scratch, for each value of the νN cross section.
C. UHE neutrino flux at the detector Figure 10 illustrates the effect of the propagation through Earth on the benchmark flux model 4, for an example arrival direction. We choose a direction from below the horizon because the column density traversed along it is large enough that changes in the νN cross section imprint sizeable changes in the flux that reaches the detector. For ν µ , and also for ν e , not shown, even the central value of the BGR18 νN cross section (σ std ) is large enough to suppress the flux at the detector by more that one order of magnitude relative to the flux at the surface of the Earth, at the highest energies. Larger cross sections vanish the flux altogether. For ν τ , the suppression is mitigated, though not counterbalanced, by the pile-up of low-energy, regenerated ν τ . For downgoing neutrinos (θ z 80 • ), not shown, the fluxes are unaffected even by large cross sections, due to the small column densities. Conversely, for upgoing neutrinos (θ z 100 • ), not shown, the fluxes vanish even if the cross section is small, due to the large column densities. Figure 11 shows how the flavor ratios, i.e., the proportion of ν α +ν α to the all-flavor flux, are affected by the propagation through Earth, also for benchmark flux model 4. Above 10 8 GeV, neutrinos of all flavors are attenuated equally, and so their flavor ratios at the detector are equal; at these energies, they also match the flavor ratios at the surface of the Earth; see Fig. 6. The effects of the propagation become apparent at lower energies. There, the neutrinos of all flavors regenerated from NC interactions plus the ν τ regenerated from CC interactions pile up. As a result of the latter, the tau-flavor ratio dominates over the electron-and muon-flavor ratios. The dip in the electron-flavor ratio around 6.3 PeV is due to the Glashow resonance experienced only byν e . There are corresponding bumps at the same energy in the muon- and tau-flavor ratios. The above features become more prominent the higher the νN DIS cross section. While these features are an implicit part of our analysis, we do not make use of flavor identification, since the capabilities of IceCube-Gen2 in this direction are still under study [82,83,182]. Section IX comments on this possibility for future work.

VI. RADIO-DETECTING UHE NEUTRINOS
We gear our forecasts to the case of UHE neutrino radio-detection in the planned radio array of IceCube-Gen2 [39]. Once neutrinos have propagated through the Earth and reached the detector, they might interact inside the detector volume via NC and CC νN DIS. See Section III for details on the νN DIS process. In most of these interactions, a substantial fraction of the neutrino energy is transferred into a high-energy particle shower. As the shower develops, it produces coherent radio emission that might be detected by antennas of the array.

A. Neutrino-induced radio emission in ice
In the NC or CC DIS interaction of a neutrino or anti-neutrino of energy E ν of any flavor, the final-state hadrons initiate a hadronic shower, rich in pions and muons [183]. The energy of the hadronic shower is E sh = (1 − y)E ν , where y is the inelasticity; see Section III. In the CC DIS interaction of a ν e orν e , the final-state electron or positron initiates an additional electromagnetic shower, rich in photons, electrons, and positrons, co-located with the hadronic shower. The energy of the electromagnetic shower is E sh = yE ν .
At neutrino energies below roughly 10 9 GeV, the hadronic and electromagnetic showers develop in phase and appear as a single shower with energy E sh = E ν ; see, e.g., Ref. [80,182]).
At higher energies, the electromagnetic shower is subject to the Landau-Pomeranchuk-Migdal (LPM) effect [184,185], which reduces the cross section of the high-energy electrons and positrons. The precise role of the LPM effect in the radio-detection of neutrinos is under study [182], but seems to be significant: if present, it delays the first interactions of electrons and positrons in the shower or leads to multiple spatially displaced sub-showers. As a result, the hadronic and electromagnetic showers may develop differently, and the shower energy might not match the neutrino energy anymore. The simulations of ν e -induced CC showers that we perform to describe the detector response, described in Section VI B, include the LPM effect. Nevertheless, when computing ν e -induced CC event rates, as described in Section VI D, we maintain the relation E sh = E ν across all energies, since further work is needed to find an equivalent form of it that accounts for the changing dominance of the LPM effect with energy.
Separately, in the CC DIS interactions of ν µ and ν τ , we also ignore the contribution of secondary interactions of final-state muons and tauons to the event rates, because they are challenging to simulate. Since their inclusion would increase the event rates by up to 25% [82,83], depending on energy and spectral shape, ignoring them makes our forecasts conservative. Future revised estimates might isolate the contribution of the LPM effect and include secondary leptons, via changes to the relation between shower and neutrino energies, Eq. (14), and to the simulated detector effective volume; see below.
IceCube-Gen2 will be built in the Antarctic ice; see below for details. Because ice is dielectric, as a neutrinoinduced shower develops inside it, it builds up a timevarying negative charge excess in the shower front which produces coherent radio emission. This Askaryan radiation [186] is strongest when the shower is observed along a cone with half-angle of arccos(1/n) ≈ 56 • , centered on the shower axis, where n = 1.78 is the index of refraction of deep ice. Along this "Cherenkov angle", the radiation emitted by the shower interferes constructively. The radio signal is a broadband, bipolar pulse a few nanoseconds long, predominantly in the frequency range 100 MHz-1 GHz. While the shower track itself is only a few tens of meters long, the radio emission can propagate over kilometers. However, the emission strength decreases quickly if the shower is observed from angles smaller or larger than the Cherenkov angle, even from only a couple of degrees away from it. Reference [187] first pointed out that Askaryan radiation can be used to detect neutrinos. Reference [188] contains an in-depth description of radio emission from high-energy particles.
Accelerator measurements have demonstrated the existence of in-ice Askaryan emission, and found agreement with theoretical predictions [189][190][191]. Additional evidence comes from the observation of radio emission from extensive air showers-i.e., particle showers in the atmosphere-to which Askaryan radiation contributes in a sub-dominant capacity [188,[192][193][194].

B. The radio component of IceCube-Gen2
The main advantage of radio-detection of UHE neutrinos is the large attenuation length of radio signals in polar ice, of roughly 1 km [195], versus the attenuation and scattering lengths of optical Cherenkov light, of roughly 100 m. This allows for a cost-efficient instrumentation of huge detector volumes with a sparse array of compact radio-detection stations. Each station contains several antennas buried in the ice at depths of up to 200 m and acts as an autonomous unit. To maximize the overall sensitivity of the array, the detector stations are separated by more than 1 km, so that coincidences between stations are rare. Each detector stations measures enough information to determine the shower energy and neutrino arrival direction; see, e.g., Refs. [182,[196][197][198][199][200].
The envisioned design of IceCube-Gen2 [39] capitalizes on this opportunity by including an array of radio antennas covering a total surface area of 500 km 2 [43]. The radio array will be of critical importance in the quest for the discovery of UHE neutrinos. Its design combines the advantages found in the pathfinder experiments ARA [201], ARIANNA [202], and RNO-G [35] into a hybrid array. The IceCube-Gen2 array consists of two types of radio detector stations that measure and reconstruct neutrino properties with complementary accuracy, intended to maximize the discovery potential by mitigating risks and adding multiple handles for the rejection of rare backgrounds [43]. The design blends the hybrid stations explored by RNO-G [35]-narrow bicone and quad-slot antennas on three strings up to a depth of 200 m plus high-gain log-periodic-dipole-array (LPDA) antennas close to the surface-with shallow-only stations explored by ARIANNA-LPDA antennas close to the surface with one additional dipole antenna at a depth of 15 m to aid event reconstruction [203].
A key ingredient in our estimates of the neutrinoinduced event rates below is the detector response of the IceCube-Gen2 radio component. We compute it via numerical simulations using the same open-source tools that are used by the IceCube-Gen2 Collaboration, Nu-RadioMC [80] and NuRadioReco [81]. Below we describe how we simulate the response of one radio station. After, we scale up its response, expressed in terms of the effective volume, to the size of the full radio array.
NuRadioMC simulates the neutrino interaction in the ice, the generation of the radio emission and its propagation to the antennas, and performs a full detector and trigger simulation. We simulate showers of varying energy, E sh , that enter the detector from different directions, cos θ z . To compute the Askaryan emission, we adopt the ARZ prescription [204], combined with a representative shower library of charge-excess profiles [80], which provides realistic modelling of the LPM effect for ν e CC interactions. The deep antennas of the station are triggered by a four-channel phased array [35,205] and the shallow antennas by a high/low-threshold trigger with a time-coincidence trigger requiring coincident detection of two out of the four LPDA antennas, using an optimized trigger bandwidth [206]. We use the exact same simulation settings as in Ref. [43].
We report the detector response via its effective volume, i.e., the simulation volume multiplied by the ratio of the number of detected showers to the number of simulated showers; see, e.g., Ref. [207]. We do this separately for the shallow and deep detector components of the simulated station, and separately for the NC interactions of all neutrino flavors and CC interactions of ν µ and ν τ , and for the CC interactions of ν e . Results for ν α andν α are the same, since their UHE νN DIS cross sections and inelasticity distributions are nearly equal; see Section III and Refs. [63,74,121]. The resulting effective volumes are functions of E sh and cos θ z .
There are two notable differences in the effective volumes that we have generated compared to previously reported results. First, we do not fold in the effects of in-Earth propagation into the effective volumes. Those effects are instead imprinted on the neutrino flux that arrives at the detector; see Section V. As a result, the directional dependence of the effective volumes is due solely to the geometry of the detector components and the propagation of radio signals in the Antarctic ice. This facilitates exploring the effect of using different values of the cross section during neutrino in-Earth propagation and detection. Second, the effective volumes that we use are functions of shower energy, not of neutrino energy. In our calculation of shower rates below, this choice allows us to account for all possible combinations of E ν and y that produce a shower of a given energy E sh (see Section III A), and to assess the detectability of that shower. Figure 12 shows the resulting effective volumes for single stations and their dependence on shower energy and direction. The effective volumes are higher for higher shower energy, since the radio signal intensifies with en-ergy. The directional dependence of the effective volumes is more nuanced. Broadly stated, shallow antennas have a higher effective volume around the horizon (cos θ z ≈ 0) and deep antennas have a higher effective volume for downgoing directions (cos θ z 0). This is because deeper antennas are less affected by shadowing due to the changing index of refraction in the 200 m of overhead ice that leads to a downward bending of signal trajectories; see, e.g., Ref. [80]. Since the capability to measure the UHE νN cross section is contingent on the observation of near-horizontal events (see Section III C), shallow antennas are key. Later, in Fig. 19, we quantify their importance. (In Fig. 12, the effective volume for the deep detector component dips at cos θ z ≈ −0.3 because the beams of the phased-array trigger system have not been optimized for upgoing neutrinos, as they will anyway be attenuated by propagating through the Earth; see Section V and Fig. 14.) Figure 13 shows the effective volume of the full IceCube-Gen2 radio array. We adopt the baseline array design from Ref. [43]: 144 hybrid stations, containing shallow and deep detector components, plus 169 shallowonly stations. Because few showers are expected to be detected simultaneously by multiple radio stations [43], it is straightforward to scale the single-station effective volumes up to the size of the full array: it is equal to the single-station volume of shallow antennas times the number of stations that contain shallow antennas (313) plus the single-station volume of deep antennas times the number of stations that contain deep antennas (144). Below, when forecasting event rates, we use the full-array effective volumes from Fig. 13, computed separately for the NC interaction of all neutrino flavors or the CC interaction of ν µ and ν τ , V NC eff,να , and for the CC interaction of ν e , V CC eff,νe . At the time of writing, the number of shallow and deep stations in the IceCube-Gen2 array is under consideration; our results might provide guidance to optimize it for cross-section measurements.

C. Reconstruction of energy and direction
The reconstruction of the shower direction and energy requires the measurement of the signal arrival direction, polarization, viewing angle, i.e., the angle under which the shower is observed, and the distance to the neutrino interaction vertex, as well as accurate modelling of the ice properties to correct for the bending of the signal trajectories due to the changing index of refraction [196].
The signal arrival direction can be determined to subdegree precision via the signal arrival times to the different antennas of a detector station [199,208,209]. The index-of-refraction profile is known well enough to achieve sub-degree precision on the signal direction, as tested in measurements at the South Pole [208]. The viewing angle can typically be reconstructed to within 1 • [197,200,210]. The dominant uncertainty on the direction comes from the measurement of the signal polar-  Table I and Fig. 15, and in our main results, we use the baseline choices for the energy and angular resolution, σ = 0.1 and σ θz = 2 • , respectively. See Section VI for details.
ization [197,208,210]. Polarization is generally measured better in shallow detector components because their orthogonal LPDA pairs provide equal sensitivity to both polarization states. For a shallow station, the capability to measure the polarization has been tested in-situ [208] and via measurements of cosmic rays [211]. The ability to measure the viewing angle and to combine all individual measurements to estimate the arrival direction was quantified in simulation studies using the forward folding technique [81,197,199,210] and deep neural networks [182,212]. The best resolution of the direction was obtained with a shallow detector station with a 68% quantile of 3 • , which translates into an uncertainty of the zenith angle of σ θz = 2 • . We adopt this angular resolution as the baseline in our forecasts below. Later, in Fig. 21, we explore how our results depend on the angular resolution.
Similar studies have estimated the uncertainty of the reconstructed shower energy [200,203,210]. They yield a resolution of approximately 0.1 in log 10 (E rec sh /E sh ), where E rec sh is the reconstructed shower energy and E sh is the real shower energy, or roughly 30% on a linear scale in the energy range of interest. We adopt this energy resolution as the baseline in our forecasts below.

D. Neutrino-induced event rates in IceCube-Gen2
In a realistic experimental setting, only the shower energy is measured, not the energy of the neutrino that initiated it. Thus, in predicting the detected shower rate, we must account for the fact that a shower measured with a certain energy E sh could have been initiated by any combination of neutrino energy, E ν , and inelasticity, y, that satisfies E sh = E ν y for NC interactions of all neutrino flavors and CC interactions of ν µ and ν τ , or E sh = E ν for CC interactions of ν e . See Section VI A for details. The relative contribution of each possible combination is weighed by the neutrino flux-which determines how many neutrinos of energy E ν reach the detectorand by the differential cross section-which determines the chances that the interaction of a neutrino of energy E ν has inelasticity y. (While the neutrino energy may be reconstructed from the measured shower energy, doing so requires making additional assumptions [196,200], and is unnecessary for our goals, so we do not attempt that.) Thus, the differential rate of showers induced in the radio array of IceCube-Gen2 by ν α arriving at the detector with flux Φ det να after propagating inside Earth is where T is the exposure time of the detector, n t ≡ N Av ρ ice /M ice is the number density of water molecules in ice, N Av is Avogadro's number, ρ ice = 0.9168 g cm −3 is the density of ice, and M ice = 18.01528 g mol −1 is the molar mass of water. Equation (13) accounts for the contributions of NC and CC interactions. On the right-hand side, the term E NC να (E sh , y)/E sh re-scales the units of the flux from neutrino energy to shower energy. The effective volume, V NC eff,να , is the full-array volume described in Section VI B. The cross section, σ NC ναw , is for neutrino DIS interaction on a water molecule (H 2 O), i.e., σ NC ναw = 10σ NC ναp + 8σ NC ναn , where σ NC ναp and σ NC ναn are the BGR18 ν α p and ν α n cross sections, respectively, and similarly for the CC cross sections [63]. In Eq. (13), the differential shower rate on the left-hand side is evaluated at a shower energy E sh . This determines, on the righthand side, the neutrino energy at which the integrand is evaluated, i.e., for i = NC, CC. Equation (14) applies also to antineutrinos. (In our numerical solution of Eq. (13), we set the lower limit of the y integral to 10 −8 instead of 0 to prevent the integral from diverging.) The differential shower rate in Eq. (13) is the "true" shower rate. Because the detector measures energy and direction imperfectly, the true rate is unobservable. Instead, to produce our results, we use exclusively the "reconstructed" shower rate, i.e., the rate expressed in terms of the reconstructed shower energy, E rec sh , and the reconstructed direction, θ rec z . These are the shower quantities that would be observed in a realistic experiment, after accounting for the measurement uncertainties. To do this, we fold the true shower rate with resolution functions in shower energy and direction, R E sh and R θz , respectively, that capture the measurement precision of the detector.
For the resolution function in energy, we adopt a Gaussian function centered at the real shower energy, E sh , i.e., where the normalization constant, ensures that the integral of R E sh from E sh = 0 to ∞ equals 1. The width of the resolution function is σ E sh = 10 σ E sh , where σ is the spread in the ratio ≡ log 10 (E rec sh /E sh ). Reference [200] reported approximately σ = 0.1, based on simulations performed for RNO-G; we take this value to be representative also of IceCube-Gen2, and use it to produce our results; see Section VI C. References [203,210] reported similar results for a shallow detector component in ARIANNA [203,210]. The choice of σ has little effect on our results, as long as it is under control, since the cross section is extracted mainly from the angular distribution of the showers, rather than from their energy distribution; see Section III C.
For the resolution function in direction, we adopt a Gaussian function centered at the real direction, θ z , i.e., where the normalization constant, ensures that the integral of R θz from θ z = 0 to π equals 1. We use σ θz = 2 • to produce our main results; see Section III C. In Section VIII, we study the effect of varying σ θz on measuring the νN cross section. Folding the resolution functions, Eqs. (15)- (18), with the true shower rate, Eq. (13), yields the reconstructed shower rate, i.e., Finally, the total shower rate is the contribution from ν α andν α of all flavors, i.e., We use Eq. (20) to forecast shower rates throughout our analysis, integrated in bins of E rec sh and cos θ rec z . In Section VII C, we comment on the effect of the choice of binning on the measurement of the cross section. Figure 14 shows the mean expected number of neutrino-induced showers for benchmark UHE neutrino flux model 4 [61] as illustration, broken down into the contribution of each interaction type and neutrino species. Because flux model 4 has comparable fluxes of each neutrino species at the surface of the Earth (see Fig. 6), the differences between the rates of different channels in Fig. 14 are due mainly to differences in the in-Earth propagation of the different species, the connection between neutrino and shower energy, Eq. (14), and the effective volumes for ν e CC interactions and for all other interaction channels (see Fig. 13). Figure 14 displays features that are common to all interaction channels and all benchmark flux models. Barring differences in the flux of the different species that reach the detector, the following observations hold generally. First, the dominant contribution comes from the CC interactions of ν e andν e , the only two cases for which E sh = E ν ; see Eq. (14). (In Fig. 14, there are more showers in the range 10 9 -10 10 GeV because those showers receive the full neutrino energy, and the flux model 4 peaks within that range; see Fig. 5.) Second, the NC interactions of all species, and the CC interactions of ν µ ,ν µ , ν τ , andν τ , contribute at about the same level, since E sh = yE ν for all of them; see Eq. (14). Third, in all channels, because of in-Earth attenuation (see Section V), events are mainly downgoing, i.e., 0.2 ≤ cos θ rec z 1, and, to a lesser extent, near-horizontal, i.e., −0.2 ≤ cos θ rec z ≤ 0.2. Fourth, all channels have a deficit around cos θ rec z = 1, because the neutrino interaction takes place almost always below the deepest antenna of the station, and all radio signals are emitted downward for vertical neutrino directions; see Fig. 13. Fifth, and finally, all channels have a deficit around E rec sh = 10 7 GeV, because the detector loses sensitivity quickly at low shower energies, where radio emission is weak; see Fig. 13. Table I shows that the IceCube-Gen2 all-sky shower rate grows with the νN cross section, regardless of the choice of benchmark flux model. On the one hand, for all values of the cross section, the shower rate is dominated by downgoing showers, with θ rec z ≤ 80 • , as expected from I. Expected rates of neutrino-induced showers in the radio component of IceCube-Gen2, after T = 10 years of exposure time, for the benchmark UHE diffuse neutrino fluxes used in this analysis (see Section IV), and varying the νN DIS cross section σ ≡ σν αN , where α = e, µ, τ . The cross section σ std is the central BGR18 prediction [63] (see Section III), against which we measure variations. The variation is the same for all flavors of να andνα. Possible flux types are: extrapolation to ultra-high energies (•), cosmogenic ( ), source ( ), and cosmogenic + source ( ). Results are obtained using the baseline choices for the detector resolution: shower energy resolution of σ = 0.1, with ≡ log 10 (E rec sh /E sh ), and angular resolution of σ θz = 2 • ; see Section VI D for details. In this table, the shower rates shown are binned in a single bin of reconstructed shower energy, 10 7 ≤ E rec sh /GeV ≤ 10 10 . Rates are all-sky, i.e., summed over all values of reconstructed direction, −1 ≤ cos θ rec z ≤ 1; we show separately the fraction of showers that are upgoing or near-horizontal, i.e., with θ rec z > 80 • . (To produce our main results, instead we bin finely in reconstructed energy and direction; see Section VII C.) # Type UHE ν flux model Shower rate in IceCube-Gen2 radio array (10 yr) Section III C. On the other hand, the fraction of nearhorizontal and upgoing showers (θ rec z > 80 • ) decreases with the cross section. These two observations reveal that, in the all-sky shower rate, the main dependence on the cross section comes from its role in the detector, rather than in the in-Earth attenuation of the neutrino flux; see Eq. (1) for a simplified picture. Table I confirms what Eq. (1) professed: as the attenuation becomes more important, with rising cross section, the rate of near-horizontal and upgoing showers decreases. This is where the sensitivity to the cross section comes from. Table I shows that the changes in the shower rate with the cross section are roughly comparable for all benchmark flux models; differences between them are due to their different energy spectra. Figure 15 illustrates the behavior seen in Table I for one benchmark flux model; other models behave similarly. The shape of the shower distribution in E rec sh is largely insensitive to the value of the cross section: changes in the cross section merely scale the distribution down-if the cross section is smaller-or up-if the cross section is larger. This is because the all-sky shower rate is dominated by downgoing showers, initiated by neutrinos whose flux is unattenuated due to traveling a short distance inside Earth. Equation (1) shows that for downgoing showers, because there is little attenuation, the rate scales roughly linearly with the cross section. In contrast, the shape of the shower distribution in cos θ rec z is strongly affected by the value of the cross section. Equation (1) shows that, due to the exponential attenuation, a smaller cross section mitigates in-Earth attenuation and flattens the distribution in cos θ rec z , while a larger cross section enhances the attenuation and tilts the distribution in cos θ rec z . This exponential angular dependence of the shower rate on the cross section-rather than the linear dependence of the all-sky shower rate-is where most of the sensitivity to the cross section comes from. The tilting of the angular distribution becomes stronger the larger the value of θ rec z , i.e., the longer the distance neutrinos travel inside the Earth. However, for upgoing neutrinos, the flux is nearly fully attenuated, even for the most optimistic diffuse flux models. Therefore, in practice, the sensitivity to the cross section comes from the angular distribution of showers coming from around the horizon. We pay special attention to them later, in Section VII C, by binning them in fine angular bins.

VII. MEASURING THE CROSS SECTION
We forecast the capability of IceCube-Gen2 to measure the UHE νN DIS cross section by means of its effect   Tables I and II). Each column shows a different choice of the νN DIS cross section that affects the propagation of neutrinos through the Earth and their detection. The cross section σ std is the central value of the BGR18 calculation [63], which we take as our baseline. See Section VI D for details about the calculation of shower rates. We include the background of showers due to atmospheric muons, which factors into our statistical analysis, computed using the hadronic interaction model Sybill 2.3c [175]. (In the lower row of panels, the background is multiplied by a factor of 10 to make it visible.) In our analysis, we use exclusively the background after the application of the surface veto ("veto on"); in this plot, we show the case without the veto ("veto off") only for comparison.
on the predicted shower rates. Because of the degeneracy that exists between the cross section and the neutrino flux normalization when computing shower rates (see Section III C 2), our forecasts are always for their simultaneous measurement. We use the methods introduced in Section VI D to generate samples of showers for different choices of their values. We perform the analysis below separately for each of the UHE neutrino flux models 1-12 introduced in Section IV.

A. Overview
For a choice of flux model, we start by generating the "real" shower sample, i.e., a sample of showers distributed and binned in E rec sh and cos θ rec z , computed using the "real" values of the cross section and flux normalization. For the cross section, the real value is given by the central value of the BGR18 DIS calculation, σ std . For the flux normalization, the real value, Φ 0,std , is given by the model flux evaluated at a reference energy of E ν,0 = 10 8 GeV. The resulting real shower sample is the mean expected sample predicted for IceCube-Gen2. Later, in Section VII D, we use a statistical procedure to assess how well we can recover the real values of the cross section and flux normalization in a realistic experimental setting. As part of that procedure, we generate test shower samples using values of the cross section, σ ναN , and flux normalization, Φ 0 , different from the real values, and compare those shower samples to the real one, including the presence of a non-neutrino background.
To generate test shower samples, we parametrize deviations from the real value of the νN DIS cross section, within the range E ν = 10 7 -10 10 GeV, via We consider only energy-independent modifications of the cross section, as in Ref. [20], so f σ is constant in energy. Changing f σ shifts the value of the cross section up or down from the central BGR18 value. The same value of f σ applies to all flavors of ν α andν α , and to interactions on protons (N = p) and neutrons (N = n). The real shower sample is computed using f σ = 1. Later, in Section VII D,we allow its value to float generously, from 0.01 to 100, and generate a large number of test shower distributions for its different values. Changing f σ affects the neutrino flux that reaches the detector-by changing the propagation through the Earth (see Section V)-and the cross section at the detection (see Section VI D). For each choice of f σ , we compute both effects. In our analysis, changing f σ affects only the νN DIS cross section off the nucleon partons; the remaining, sub-leading neutrino interaction channels that act during neutrino propagation through the Earth (see Section V) remain unchanged. Section IX comments on possible generalizations.
Similarly, for each choice of flux model out of models 1-12, we parametrize deviations from the real value of the flux normalization, via where the value of Φ 0,std is unique to each flux model. Changing f Φ shifts the flux up or down from its central value, but the shape of the neutrino spectrum, flavor composition, and ν vs.ν content are unchanged. The real shower sample is computed using f Φ = 1. The same value of f Φ applies to all flavors of ν α andν α . Later, in Section VII D, we allow its value to float generously, and generate a large number of test shower distributions for its different values. We vary f Φ indirectly: for a given choice of flux model out of 1-12, we allow the value of E 2 ν,0 Φ 0 to float from 10 −12 to 10 −7 GeV cm −2 s −1 sr −1 , and compute f Φ for each value of E 2 ν,0 Φ 0 . Section IX comments on possible generalizations of this procedure.

B. Non-neutrino backgrounds
In the radio component of IceCube-Gen2, in addition to neutrino-initiated events, we expect a background of irreducible non-neutrino-induced events that could mimic a neutrino signal. There are two main such backgrounds: in-ice showers induced by high-energy atmospheric muons created by cosmic-ray interactions in the atmosphere [82], and air-shower cores, i.e., cores of showers that start in the atmosphere, induced mainly by cosmic rays, and that continue developing downwards in the ice [213]. The latter lead to in-ice Askaryan radiation that can be reflected back upwards by reflection layers that are known to be present in deeper ice. In our forecasts, we account only for the background of showers due to atmospheric muons, by generating their expected distribution in reconstructed energy and direction, as described below. We do not account for the background of cosmic-ray air-shower cores, since their estimates are presently uncertain, but comment on them in Section IX. Figure 16 shows the mean distribution of muoninduced shower events in the radio component of IceCube-Gen2 [43,82,83]; see also Fig. 15. We estimate it using the hadronic interaction model Sybill 2.3c [175] and apply a veto coming from the detection of atmospheric muons by a surface array of cosmic-ray detector tanks, that helps to mitigates the rate. We expect roughly 0.2 showers per year in the range 10 4 -10 11 GeV, and < 0.1 showers per year in the range 10 7 -10 10 GeV relevant for our forecasts. These rates are subject to significant uncertainties, due mainly to the high-energy muon flux, but the ability to measure the cross section in IceCube-Gen2 will not be affected by the uncertainties; we comment on this below. The background of showers from atmospheric muons is direction-and energydependent. Figure 15 shows that most of the showers are in the range 10 7 -10 9 GeV; thus, the highest-energy end of the shower spectrum is nearly free from this background. Since it is unlikely that a muon produced in the atmosphere reaches the detector from below the horizon, the background for θ z > 90 • is zero (this changes slightly after applying angular smearing below).
Similarly to how we use the reconstructed neutrinoinduced shower rate (see Section VI D) in our analysis, we also use the reconstructed rate of muon-induced showers; this is what Fig. 16 shows. We compute it in the same way: by convolving the true rate of muon-induced showers with the resolution functions in Eqs. (15)- (18), and expressing the result in terms of E rec sh and cos θ rec z . At all times, we use the same energy and angular resolution for muon-induced showers as for neutrino-induced showers. As a result, a tiny fraction of the muon-induced showers leak into the highest energies and into directions just below the horizon, i.e., for θ rec z 90 • ; see Fig. 15.

C. Choice of shower binning
To produce our forecasts, we bin the real and test shower samples (see above) in fine bins of E rec sh and cos θ rec z . Accordingly, in our statistical analysis we use a binned likelihood function; see Section VII D. The bins of reconstructed shower energy and direction that we use are equal in size to or larger than the width of the detector energy and angular resolution, i.e., σ θz and σ E sh , respectively, from Eqs. (15) and (17). Section IX comments on the possible future use of an unbinned likelihood.
When binning in cos θ rec z , we pay special attention to Earth-skimming neutrinos, which have the strongest sensitivity to the cross section; see Section III C. For downgoing showers, we expect a large rate; but, since they are largely unattenuated by their passage through Earth, they do not carry a significant angular dependence induced by the cross section. Thus, for them we use a single bin covering θ rec z ∈ [0 • , 80 • ]. For horizontal and near-horizontal showers, our region of interest, we use ten bins, equally spaced between θ rec z = 80 • and 100 • , each with a bin size equal to the angular resolution of the detector, σ θz . To produce our main results, we use the baseline angular resolution of σ θz = 2 • ; later, we explore alternative choices (see Fig. 21). For upgoing showers, of which there are few due to the attenuation of the flux inside Earth, we use two large angular bins: one for θ rec z ∈ [100 • , 110 • ], which is populated only if the flux is high and the cross section is low (see Fig. 15), and one for θ rec z ∈ [110 • , 180 • ], which is typically nearly empty. For upgoing showers, using finer binning would lead to a large number of empty, uninformative bins.
When binning in E rec sh , we use twelve bins equally spaced in logarithmic scale from log 10 (E rec sh /GeV) = 7 to 10. Even though in our analysis we allow only for energy-independent re-scaling of the cross section from its predicted central value (see Section VII A), the distribution of showers in energy still matters when extracting the cross section, for three reasons. First, during neutrino propagation inside the Earth an energy-independent rescaling of the cross section affects the flux attenuation differently at different energies, due to the non-linear dependence on energy of the regeneration of neutrinos in NC interactions and ν τ regeneration; see Fig. 8. Second, the distribution of showers in energy reflects the distribution of neutrino energies in the diffuse flux. Third, showers from the background of atmospheric muons are concentrated at low energies: above a few times 10 7 GeV, there are nearly none, so binning in energy enhances the statistical power of the neutrino-induced showers in that range. Naturally, the choice of energy binning is relevant only when the predicted number of showers is relatively large. If shower rate is low, like for benchmark flux models 1, 3, and 5, we have checked that the choice of E rec sh binning does not affect cross-section measurements.

D. Statistical analysis
To produce our forecasts, we use a Bayesian statistical analysis. Section VII A outlined the procedure; below, we present it in detail. We carry the procedure separately for each UHE neutrino flux model, 1-12 (see Section IV).
First, for a given choice of flux model we generate the real shower sample (r), i.e., the mean expected number of showers assuming the real values f σ = 1 and f Φ = 1 (see Section VII A), binned in E rec sh and cos θ rec z as described in Section VII C. In the i-th bin in E rec sh and the j-th bin in cos θ rec z of the real shower sample, the mean number of neutrino-induced showers isN (r) ν,ij . To this bin, we associate a Poisson probability distribution, with meanN (r) ν,ij , of observing N (r) ν,ij showers in a particular realization of the shower distribution, i.e., Similarly, in each bin, the mean number of background muon-induced showers isN (r) µ,ij , and we associate to it a Poisson probability distribution, given by Eq. (23) but withN µ,ij (θ). Given a specific random realization of the real shower distribution, generated following the above procedure, we quantify its compatibility with the mean test shower distribution expected for a particular choice of θ, in each bin, via the likelihood function The full likelihood function is the product of all the partial likelihood functions, i.e., where N E rec sh and N θ rec z are, respectively, the number of bins in E rec sh and θ rec z ; see Section VII C for details. The posterior probability distribution is where π(θ) ≡ π(f σ )π(f Φ ), π(f σ ) is the prior on f σ , and π(f Φ ) is the prior on f Φ . We use wide priors in log 10 f σ and log 10 f Φ to facilitate exploring large possible deviations from the real values of f σ = 1 and f Φ = 1. The priors are flat to avoid introducing unnecessary bias. For log 10 f σ , the prior spans the interval from -2 and 2; this range covers all reasonable modifications of the cross section. For log 10 f Φ , the prior spans a different interval for each flux model. For a given flux model, it is defined by varying the flux at E ν,0 = 10 8 GeV from E 2 ν,0 Φ(E ν,0 ) = 10 −12 to 10 −7 GeV cm −2 s −1 sr −1 and using Eq. (22) to compute f Φ . The low end of this interval corresponds, in practice, to an unobservable flux. The high end exceeds the upper limits on the flux from IceCube and Auger (see Fig. 5). However, this does not put our analysis in tension: the IceCube and Auger limits were obtained under the assumption that the νN cross section has its nominal value (f σ = 1) while, in our analysis, when we explore large reductions of the cross section (f σ 1), they can be traded off for large enhancements of the flux while keeping the number of showers unchanged; see Section III C.
For a particular random realization of the real shower distribution, we compute the posterior, Eq. (26), as a function of f σ and f Φ . We generate a large number of different random realizations of the real shower distribution using the above prescription; for each, we compute the associated posterior. Since our goal is to produce forecasts for the average sensitivity of IceCube-Gen2, we take the arithmetic average, P , of all the posteriors thus generated. This is equivalent to marginalizing over the Poisson distribution in each bin of the real shower distribution. In this way, our forecasts reflect the average sensitivity of the detector, taking into account statistical fluctuations. In an actual detection in IceCube-Gen2, and assuming that f σ = 1 and f Φ = 1 are indeed the real values, the posterior computed with real data would be one of the random realizations that we averaged over.
Below, to produce our forecasts, we maximize P with respect to f σ and f Φ , and find their best-fit values and credible intervals. We report results either as two-dimensional posteriors (see Figs. 17 and A1), or as one-dimensional posteriors obtained by marginalizing the full posterior over one of the two parameters (see Table II and Fig. A2), i.e., P(f σ ) fσ = df Φ P(f σ , f Φ ) for f σ , and equivalently for f Φ .

VIII. RESULTS
Below, to produce our main results, we use the baseline design of the IceCube-Gen2 radio array, i.e., 144 hybrid (deep + shallow) stations plus 169 shallow-only stations, the baseline detector energy resolution, σ = 0.1, and angular resolution, σ θz = 2 • , and T = 10 years of exposure time. We state explicitly when we vary these parameters. Figure 17 shows our results for the simultaneous measurement of the cross section and flux normalization for a selection of benchmark UHE neutrino flux models that is representative of the full variety explored in our analysis; see Fig. 5. Appendix A contains similar results for all benchmark flux models. In Fig. 17, the elongated shape of the two-dimensional posteriors follows from the partial degeneracy between the cross section and the flux normalization: the number of showers depends on their product; see Eq. (1). The posteriors are narrower for flux models with a higher associated mean expected shower rate, i.e., in Fig. 17, for flux models 4 ("Bergman & van Vliet, fit to TA UHECRs") and 6 ("Rodrigues et al., all AGN"). This is because, in those cases, the number of horizontal showers is large enough to partially break the degeneracy between cross section and flux; in Eq. (1), this happens via the exponential attenuation term. Figure 17 reveals an asymmetry between the region of high flux and small cross section in the upper-left corner-which is disfavored-and the region of low flux and large cross section in the lower-right corner-which is allowed. The asymmetry becomes more evident for lower-flux models. This preference for low flux and large cross section stems from the effect of the cross section on the shower distribution, as illustrated in Fig. 15 for flux model 4. A small cross section (f σ 1) implies that the in-Earth attenuation is small and, therefore, the flux at the detector, and the angular distribution of the associated showers, are closer to isotropic. In contrast, a high cross section (f σ 1) implies a tilted angular distribution, with a sharp drop of the rate around the horizon. For the real values of the cross section and flux normalization, f σ = 1 and f Φ = 1, the predicted near-horizontal All-flavor ν flux (10 8 GeV), log 10 (E 2  shower rate is low: for most fluxes, it is only about one shower per bin in 10 years (exceptionally, for flux model 4, the highest among the benchmarks, there is a handful of near-horizontal showers). Therefore, high values of f σ 1, which also result in a sharp drop, are favored over f σ 1; the asymmetry in Fig. 17 reflects this. Figure 18 shows that the sensitivity to f σ and f Φ comes predominantly from near-horizontal and horizontal showers, since they are the ones initiated by neutrinos that undergo significant, but not complete, attenuation inside the Earth. Using only downgoing showers reveals the underlying degeneracy between cross section and flux in the calculation of shower rates (see Section III C), but is unable to break it. In short, downgoing showers, more numerous, help in the discovery of the diffuse UHE neutrino flux, while near-horizontal showers, more rare, help to jointly perform studies beyond flux discovery, like crosssection measurements. Figure 19 shows that, using the baseline array designmade up of 144 hybrid (shallow + deep) plus 169 shallow-only stations-the precision in the measured values f σ and f Φ is dominated by shallow stations. The caveat to this is that in our analysis we assume that the angular resolution of all detected showers is the same, whereas in reality there will be differences. For instance, detection in deep antennas is expected to yield poorer resolution (see Section VI C), which would further boost the importance of shallow stations. Section IX comments on future improvements of the analysis in this direction. Table II shows the one-dimensional marginalized ranges of f σ and f Φ for all the benchmark flux models. The results show that the radio component of IceCube-Gen2 could indeed jointly measure f σ and f Φ , to varying precision, as long as at least a few tens of showers are detected in 10 years. Appendix A shows the corresponding one-dimensional marginalized posteriors.
We classify the results in Table II in three groups, according to the size of the flux and its associated shower rate. First, for low-flux models 1, 3, and 5, we expect fewer than one shower in 10 years and so no measurement is possible. Second, for medium-flux models 2, 7, 8, 10, 11, and 12, we expect tens of events; with them, we can measure f σ to within 50-180% and f Φ to within 70-110%, depending on the flux model. Third, for highflux models 4, 6, and 9, we expect 100-300 events; with them, we can measure f σ and f Φ to within 40%. Table II reveals that these groups are not dominated by a single flux type; rather, different flux types-extrapolated, cosmogenic, source, and cosmogenic + source (see Section IV)-are represented in all groups.
Table II reveals important nuance associated to the shape of the UHE neutrino spectrum. For instance, even though flux model 12 has a lower shower rate than flux model 7, it leads to more accurate measurements of f σ and f Φ . Figure 5 shows that model 12 peaks at lower energies than model 7. Since the νN cross section grows with energy, this implies that in model 12 neutrinos are less attenuated by their passage through Earth, and so more of them reach the detector from near-horizontal directions. Conversely, in model 7 neutrinos are higherenergy, so they are more strongly attenuated and a relatively larger number of them are downgoing. This, again, highlights the importance of Earth-skimming showers to perform studies beyond flux discovery.
We have tested that the capability of IceCube-Gen2 to measure the cross section is not limited by the uncertainty in the background shower rate of atmospheric muons. We did this by repeating our analysis using a rate of muon-induced showers five times higher than the one in Fig. 16. The results are practically indistinguishable from those in Table II. There are three reasons for this. First, for the benchmark flux models for which the cross section can be measured, the rate of muon-induced showers is tiny compared to the rate of neutrino-induced showers. Second, muon-induced showers are concentrated at lower energies, leaving the higher energies backgroundfree. Third, muon-induced showers affect almost exclusively downgoing bins, whereas the sensitivity to the cross All-flavor ν flux (10 8 GeV), log 10 (E 2 section comes from bins around the horizon. In our analysis, we set the real values of the cross section and flux normalization to be f σ = 1 and f Φ = 1; see above. However, Table II shows that, for all flux models, their average best-fit values have a consistent bias toward f σ 1 and f Φ 1 (except for low-flux models 1, 3, and 5, for which f Φ is compatible with zero due to the low shower rate). The bias toward an average bestfit f σ 1, though slight, stems from the interplay in shower rates between the single bin that contains downgoing showers-whose mean expected shower rate is relatively high-and the multiple bins that contain horizontal and near-horizontal showers-whose mean expected shower rates are relatively low. As a result, the Poisson probability distribution of the rate, Eq. (23), of the single downgoing bin is rather symmetric, while Poisson probability distributions of the horizontal bins are rather asymmetric. Therefore, in any one of the random realizations of the shower distribution that we generate as part of our statistical analysis (see Section VII D), in the downgoing bin under-fluctuations and over-fluctuations relative to the mean shower rate are equally likely, while in the horizontal bins under-fluctuations are more likely than over-fluctuations. As a result, on average, underfluctuations in the all-sky shower rate are more likely than over-fluctuations. Because lower shower rates correspond to smaller values of the cross section, via Eq. (1), the best-fit value of f σ 1 on average. Nevertheless, within 68% C.L., the inferred values of f σ and f Φ are compatible with their real ones (again, except for the lowflux models). Further, their best-fit values approach their All-flavor ν flux (10 8 GeV), log 10 (E 2 real values closer the higher the mean expected shower rate, and the longer the exposure time; see Fig. 20. Figure 1 places our results in the context of existing cross-section measurements, including in IceCube, by showcasing a selection of three optimistic, representative benchmark flux models, 2, 4, and 5. In the best-possible scenario, i.e., for flux model 4 ("Bergman & van Vliet, fit to TA UHECRs"), after 10 years we could be able to measure the cross section to within the theory uncertainty of the BGR18 prediction [63]. For lower, but still detectable, flux models, Fig. 1 shows that the measurement worsens, but remains informative. In particular, for flux model 2 ("IceCube ν µ flux (9.5 yr), extrapolated to UHE")-a conservative benchmark based on the known TeV-PeV diffuse neutrino flux-the cross section could be measured to roughly the same accuracy at ultra-high energies as it is currently measured in the TeV-PeV range. Overall, medium-flux and highflux models (see above) offer enough accuracy to test a diversity of proposed deviations of the cross section, due to new physics, e.g., large extra dimensions [73], electroweak sphalerons [69], or to non-perturbative effects in the nucleon structure, e.g., color glass condensate [70]; see Fig. 2 in Ref. [51] for an example. Figure 20 shows the expected improvement in the measurements over time, for selected medium-flux and highflux benchmark models. In the most optimistic scenario, for flux model 4 ("Bergman & van Vliet, fit to TA UHE-CRs"), the 2.5-year uncertainty in the inferred value of the cross section, f σ , of 75 %, is already comparable to the baseline 10-year uncertainty, of 35 %. The 20-year II. Expected rates of neutrino-induced showers and measurement capability of the cross section and flux normalization in the radio component of IceCube-Gen2, after T = 10 years of exposure time, for the benchmark UHE diffuse neutrino fluxes used in this analysis; see Section IV. Possible flux types are: extrapolation to ultra-high energies (•), cosmogenic ( ), source ( ), and cosmogenic + source ( ). Measurement capabilities are expressed in terms of fσ ≡ σ/σ std for the cross section, and fΦ ≡ Φ0/Φ 0,std for the flux normalization at a neutrino energy of 10 8 GeV. The real values are fσ = 1 and fΦ = 1. Results are obtained using the baseline choices for the detector resolution: shower energy resolution of σ = 0.1, with ≡ log 10 (E rec sh /E sh ), and angular resolution of σ θz = 2 • ; see Section VI D for details. The shower rates shown are all-sky, i.e., summed over all values of reconstructed direction, −1 ≤ cos θ rec z ≤ 1, and binned in a single bin of reconstructed shower energy, 10 7 ≤ E rec sh /GeV ≤ 10 10 . However, all-sky rates are only referential, and in the statistical analysis used to infer the cross section and flux normalization, we use instead binned shower rates; see Section VII D for details. For the inferred cross section and flux normalization, we show the best-fit values and 68% credible intervals of their one-dimensional marginalized posteriors. The statistical procedure returns results for log 10 fσ and log 10 fΦ; we show the corresponding results in terms of fσ and fΦ to facilitate their interpretation. uncertainty, of 22 %, is smaller than the BGR18 theory uncertainty shown in Fig. 1. Results for model 2 ("Ice-Cube ν µ flux (9.5 yr), extrapolated to UHE") are comparable. At short exposure times, increased exposure leads to drastic improvements in the measurement uncertainty because the measurement is dominated by statistical uncertainties. At long exposure times, improvements slow down, as the measurement becomes increasingly dominated by systematic uncertainties, i.e., the energy and angular detector resolution and, to a lesser degree, the background of atmospheric muons. Figure 21 shows that angular resolution, σ θz , is key to making precise measurements of f σ and f Φ , especially in medium-flux models where the degeneracy between them is more resilient on account of a lower shower rate. A good angular resolution allows the detector to distinguish between similar arrival directions around the horizon, where the joint sensitivity to f σ and f Φ comes from; see Sections VI D and VII. In Fig. 21, we repeat our analysis, but for different choices of values of σ θz . For each choice, we set the size of the angular bins used for horizontal and near-horizontal showers, with θ rec z ∈ [80 • , 100 • ], equal to σ θz . Figure 21 shows, e.g., for flux model 4 ("Bergman & van Vliet, fit to TA UHECRs"), that using a resolution of σ θz = 5 • yields f σ = 0.94 +0. 73 −0.31 and f Φ = 1.06 +0.49 −0.44 , and using σ θz = 10 • yields f σ = 0.88 +1.54 −0.38 and f Φ = 1.14 +1.28 −0.64 , an increase of 50% and 178% in the uncertainty of f σ and of 41% and 191% in the uncertainty of f Φ , respectively, compared to the baseline results with σ θz = 2 • from Table II. Results for other flux models change similarly.

IX. FUTURE DIRECTIONS
When generating our forecasts above, we have used state-of-the-art predictions and dedicated simulations to describe the incoming UHE neutrino flux, its propagation through Earth, and its detection in the radio component of IceCube-Gen2, and to jointly measure the cross section and flux normalization. Yet, there are potential improvements to be implemented in future, revised versions of our calculation. None of them represent a fundamental limitation of our present analysis. Below, we present possible future improvements, listed roughly in ascending order of implementation challenge, with the hope of sparking progress in various directions: Raising the maximum neutrino energy in in-Earth propagation.-To compute the propagation of UHE neutrinos through the Earth, we have used NuPro-pEarth [78,79], based on the state-of-the-art BGR18 IceCube-Gen2 Radio νN DIS scattering cross section [63]; see Section III A However, the version of NuPropEarth that we used was unable to propagate neutrinos beyond energies of 10 10 GeV. Fortunately, for the majority of the benchmark flux models 1-12, the flux beyond 10 10 GeV contributes marginally to the shower rate; see Fig. 5. Only for models 4 and 11 we expect that contribution to be slightly more sizeable. We recommend raising the maximum allowed neutrino energy in future versions of NuPropEarth.
Including the LPM effect in the relation between neutrino and shower energies.-The LPM effect is included in the NuRadioMC simulations of ν e -induced CC showers and radio propagation that we use to generate the detector effective volume; see Section VI D. However, the energy-dependent relevance of the LPM is not accounted for in the relation between neutrino energy, E ν , and shower energy, E sh , Eq. (14), because a parametrization of it was not available at the time of writing. As a result, for ν e -initiated CC showers, we have assumed that E sh = E ν at all energies, whereas in reality the LPM should make E sh < E ν at the highest energies. This may be improved in future work by building a suitable parametrization from dedicated simulations.
Including the contribution of secondary leptons in the effective volumes for ν µ -and ν τ -initiated CC showers.-When generating the effective volumes for ν µ -and ν τinitiated CC showers, we have ignored the contribution of the radio emission coming from showers induced by the secondary charged leptons, i.e., the final-state muons and tauons in the CC interactions, due to the additional computational expense involved in simulating them and in parametrizing the relation between neutrino to shower energy in the presence of secondary interactions; see Section VI D. Reference [82] showed that they may enhance the shower rate by up to 25%, which makes our results above conservative. Their contribution may be included in future work via dedicated simulations.
Improving the modeling of the angular resolution.-To obtain our results, we assumed the same angular resolution for all detected showers; see Section VIII. However, in reality, the angular resolution of a detected shower depends on its event quality, and on which detector component measured it, i.e., only the shallow one, only the deep one, or multiple components in coincidence. With a better modelling of the expected experimental uncertainties, our analysis can be redone using different resolution functions depending on the event class. Further, we have only considered the resolution in zenith angle, not in azimuth angle; the latter might be relevant when analyzing real shower events, with definite celestial coordinates.
Using an unbinned likelihood analysis.-In our analysis above, we binned the expected IceCube-Gen2 showers in reconstructed shower energy and direction (see Section VII C), and we used a binned likelihood function, Eq. (25), to infer the values of the νN cross section and the neutrino flux normalization. We chose the bin sizes to be no smaller than the detector energy and angular resolution. Yet, when inferring the cross section and flux normalization from actual future data detected by IceCube-Gen2, with each shower carrying its own energy and direction resolution, an unbinned likelihood would be preferable, in the style of IceCube TeV-PeV cross-section measurements in Refs. [21,22]. Such unbinned prescription may also be used in forecasts of the cross-section measurement, with modifications of our above procedure.
Including the background of cosmic ray-induced showers.-Section VII B stated that the two main nonneutrino backgrounds in the radio array of IceCube-Gen2 are showers induced by atmospheric muons, which we account for in our analysis, and air-shower cores propagating into the ice, which we do not account for, since their rate is still largely unknown. Early estimates of this background put its rate anywhere between a handful and tens of showers per year, before the application of any cuts or of a surface veto that could mitigate them. If these estimates are accurate, and if vetoes do not mitigate the rate significantly, then they might render cross-section measurements in even medium-flux models unfeasible. Ongoing work should clarify the relevance of this background.
Including nuclear effects in the cross section.-To obtain our results, we have used the BGR18 cross section built using free-nucleon PDFs. However, the nucleons with which neutrinos interact are bound in heavy nuclei; this modifies the PDFs and, therefore, the νN DIS cross section. Reference [78] showed that using nuclear- corrected PDFs in in-Earth propagation could modify the fraction of neutrinos that reach the detector by more than 50%, especially at high energies. This is comparable to or greater than the measurement precision that we have forecast; see, e.g., Fig. 1. This should motivate future forecasts of UHE cross-section measurement to devote more attention to the effect of nuclear corrections.
Looking for energy-dependent new physics in the cross section.-In our analysis above, we have looked for energy-independent modifications of the Standard Model νN DIS cross section. However, if there are beyond-the-Standard-Model contributions to the cross section, they may instead be energy-dependent and possibly even localized only within a narrow energy window. In particular, UHE neutrinos could allow us to probe new νN resonant interactions via new heavy mediators with masses up to, roughly, (m q /GeV) 0.5 × 100 TeV, where m q is the mass of the quark involved in the interaction. Reference [53] explored the effect of new energy-dependent contributions to the νN cross section, and of resonant contributions to the neutrino-electron cross section, in the context of GRAND, POEMMA and Trinity. For IceCube-Gen2, attempting to explore energy-dependent modifications of the νN cross section within the same sophisticated end-to-end framework that we have developed above quickly escalates the complexity of the tasks involved, and is best left for future work.
Jointly measuring the νN cross section, flux normalization, and spectral shape.-In our analysis above, we forecast the joint measurement of the νN cross section and flux normalization, but kept the shape of the neutrino energy spectrum fixed to that of one of the benchmark flux models 1-12 taken from the literature; see Section IV. While this procedure provides us already with useful information about the capability to measure the cross section, ultimately we would also like the shape of the neutrino energy spectrum to be determined from data, rather than assumed. To do this, future studies could exploit the fact that the sensitivity to the cross section comes almost exclusively from showers around the horizon; see Fig. 18. This suggests that we may use the large number of downgoing, unattenuated showers to measure the shape of the neutrino energy spectrum.

X. SUMMARY AND OUTLOOK
Measuring the neutrino-nucleon (νN ) cross section at ultra-high neutrino energies, above 100 PeV, offers novel insight into the deep structure of protons and neutrons, and sensitivity to potentially transformative new physics. Yet, it hinges on using ultra-high-energy (UHE) neutrinos of cosmic origin, long sought, but so far undiscovered. Fortunately, upcoming neutrino telescopes, currently under planning and testing, will have a realistic chance of discovering them in the next 10-20 years, even if their flux is low, as is expected. Accordingly, we have put forward the first detailed forecasts of the capabilities of IceCube-Gen2 [39], one of the leading upcoming neutrino telescopes, to measure the UHE νN cross section. We have endeavored to make our forecasts complete, robust, and realistic. Our results are encouraging.
The sensitivity to the UHE νN cross section stems from the strong directional dependence of the UHE neutrino flux that reaches the detector after propagating through Earth. Due to significant νN interactions underground, UHE neutrinos only reach the detector from above-where they undergo nearly no interactions before reaching the detector-and from directions around the horizon-where significant interactions imprint the νN cross section on the angular dependence of the flux. It is from these latter Earth-skimming neutrinos that the sensitivity to the UHE νN cross section comes from.
In our forecasts, we have used state-of-the-art ingredients at every stage. For the νN cross section, we adopted as a baseline the recent BGR18 deep-inelastic-scattering (DIS) calculation [63]. For the fluxes of UHE neutrinos, we explored a wide variety of predictions [29,[54][55][56][57][58][59][60][61][62], including extrapolations of the known IceCube TeV-PeV fluxes to ultra-high energies, cosmogenic neutrinos, neutrinos produced inside cosmic-ray sources, and selfconsistent models of joint production of the latter two. We treat each neutrino species separately, ν e ,ν e , ν µ , ν µ , ν τ , andν τ . For the neutrino propagation through Earth, we compute νN DIS plus sub-dominant interactions, via NuPropEarth [78]. For detection, we focus on the radio-detection of neutrino-initiated showers. We use the same simulation tools as the IceCube-Gen2 Collaboration, NuRadioMC [80] and NuRadioReco [81], to simulate neutrino interactions in the detector, ensuing particle showers, emission of coherent radio signals-Askaryan radiation-its propagation in the Antarctic ice, and its detection in the radio stations. We account for the detector resolution in shower energy and direction, and for the background of showers from atmospheric muons.
To properly represent the significant uncertainty in the predicted neutrino flux, in our forecasts we have jointly extracted the cross section and the neutrino flux normalization. To do this, we used a Bayesian statistical approach that accounts for random fluctuations in the number of expected detected showers. We report results in terms of mean sensitivity, averaged over many random realizations of the predicted shower distributions.
We find that it may be possible to measure the UHE νN cross section to within 50% of the BGR18 prediction within 10 years of exposure of IceCube-Gen2, as long as at least a few tens of neutrino-induced showers are detected in this period; see Fig. 1. In the most optimistic case, we expect that comparable precision could be achieved within 5 years; see Fig. 20. By far, the largest systematic uncertainty in our forecasts comes from the UHE neutrino flux. Regardless, the level of precision that may be achieved is enough to test the standard prediction of the νN DIS cross section [63], to look for non-linear effects in the nucleon structure [64][65][66][67], color-glass condensates [68], and sphalerons [69], and to identify deviations introduced by a host of new-physics models [53,[73][74][75][76][77].
Our forecasts are comprehensive, but we have identified areas where further work is needed. Most pressing are estimating the effect of the background of showers induced by cosmic rays, improved modelling of reconstruction uncertainties based on individual event quality, and extending our analysis to let both the neutrino flux normalization and the shape of the neutrino energy spectrum be inferred from data. Thanks to the flexibility of our calculation framework, these improvements can be easily incorporated. Work in these directions is ongoing.
We provide our forecasts in the hope of informing the design choices of IceCube-Gen2 and other upcoming neutrino telescopes. Ultimately, we hope that detailed forecasts like ours provide guidance to tap into their full potential for fundamental-physics research.
well as IJCLab, CEA, IPhT, APPEC, the IN2P3 master project UCMN, and EuCAPT. This work used resources provided by the High Performance Computing Center at the University of Copenhagen. The computations and data handling were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX partially funded by the Swedish Research Council through grant agreement no. 2018-05973.
Appendix A: Posterior distributions for all benchmark UHE neutrino flux models In the main text we used twelve benchmark UHE neutrino flux models to make our forecasts of cross-section measurements; see Section IV. For all of them, we showed expected shower rates in IceCube-Gen2 in Table I and one-dimensional marginalized posteriors of f σ and f Φ in Table II, but we showed two-dimensional posteriors of f σ and f Φ only for four of them, models 2, 4, 6, and 7, in Fig. 17. Figure A1 shows the two-dimensional posteriors of f σ and f Φ for all benchmark flux models, for the same baseline detector parameters as in Fig. 17. From top to bottom and left to right, results are arranged in order of increasing sensitivity, from the most pessimistic to the most optimistic scenarios. Figure A2 shows the corresponding one-dimensional marginalized posteriors of f σ and f Φ . The 68% C.L. ranges in this figure match those shown in Table II.  Flux normalization, FIG. A2. One-dimensional marginalized posteriors of the UHE νN cross section, fσ, and the UHE neutrino flux normalization, fΦ, for all the benchmark UHE neutrino flux models adopted in our analysis. See Fig. A1 for the corresponding two-dimensional posteriors and Table II for numerical values.