Nuclear effects in high-energy neutrino interactions

Neutrino telescopes like IceCube, KM3NeT and Baikal-GVD offer physicists the opportunity to study neutrinos with energies far beyond the reach of terrestrial accelerators. These neutrinos are used to study high-energy neutrino interactions and to probe the Earth through absorption tomography. Current studies of TeV neutrinos use cross sections which are calculated for free nucleons with targets which are assumed to contain equal numbers of protons and neutrons. Here we consider modifications of high-energy neutrino interactions due to two nuclear effects: modifications of the parton densities in the nucleus, referred to here as shadowing, and the effect of non-isoscalar targets, with unequal numbers of neutrons and protons. Both these effects depend on the interaction medium. Because shadowing is larger for heavier nuclei, such as iron, found in the Earth's core, it introduces a zenith-angle dependent change in the absorption cross section. These modifications increase the cross sections by 1-2\% at energies below 100 TeV (antishadowing), and reduce it by 3-4\% at higher energies (shadowing). Nuclear effects also alter the inelasticity distribution of neutrino interactions in water/ice by increasing the number of low inelasticity interactions, with a larger effect for $\nu$ than $\bar\nu$. These effects are particularly large in the energy range below a few TeV. These effects could alter the cross sections inferred from events with tracks originating within the active detector volume as well as the ratio $\nu/\bar\nu$ inferred from inelasticity measurements. The uncertainties in these nuclear effects are larger than the uncertainties on the free-proton cross sections and will thus limit the systematic precision of future high-precision measurements at neutrino telescopes.


I. INTRODUCTION
Neutrino telescopes have observed neutrinos with energies well above 10 6 GeV [1]. Future experiments will use radio-Cherenkov techniques to search for neutrinos with energies up to 10 20 eV and beyond [2] due to interactions of ultra-high energy cosmic-rays with the cosmic microwave background radiation. These neutrinos are important probes of the cosmos; Astrophysical neutrinos should point back to the locations of high-energy cosmic-ray accelerators in the universe [3] while measurements of the diffuse, 4π, flux are sensitive to the properties of the accelerators [4]. Atmospheric neutrinos are sensitive to the composition of cosmic rays and also to some aspects of hadron physics [5,6]. High-energy neutrinos are also used to study neutrino interactions and to search for beyond the standard model (BSM) physics [7]. The 1 km 3 IceCube neutrino telescope [8] has collected large samples of neutrino events, sometimes comprising more than 500,000 neutrino events [9]. Optical sensors observe the Cherenkov radiation from relativistic charged particles produced in the interactions. These events may be through-going muons from neutrino interactions outside the detector, or starting events, neutrino interactions within the detector. More complex topologies are also possible. Through-going muon analyses usually bin the data by muon energy (a proxy for neutrino energy) and zenith angle then fitting these data to models making different assumptions about the astrophysical neutrino flux and angular distribution (isotropicity), possibly including neutrino propagation and interaction, including BSM interactions. The muon energy in the detector is inferred by measurement of its specific energy loss, dE/dx [10].
At energies above a few TeV, absorption in the Earth alters the zenith angle distribution. The absorption correction should be modeled with a precision similar to the other uncertainties in the zenith angle distribution, reaching a few percent. It is also important to accurately model the neutrino inelasticity distribution since it is important for inferring the neutrino energy spectrum from the muon energy spectrum. In charged-current ν µ and ν µ starting event interactions, both the hadronic cascade and the outgoing lepton are observed.
Thus inelasticity becomes a third variable in astrophysical neutrino flux fits. The energy of the hadronic cascades is inferred through calorimetry [11].
Natural neutrinos are also used to study neutrino interactions at energies above those that are accessible at accelerators, where the highest available neutrino energy is 500 GeV.
IceCube has measured the neutrino-nucleon cross sections for ν µ [12] and for a primarily ν e mixture [13,14]. They have also measured the neutrino inelasticity distribution [6] and provided interaction data used to produce a tomographic image of the mass distribution within the Earth [15,16]. These data have been used to search for BSM processes [17,18].
The cross section and Earth tomography analyses rely on measurements of neutrino absorption within the Earth while the inelasticity measurements and most BSM studies use direct observations of neutrino interactions in Antarctic ice. The tomography effort is noteworthy because seismographic measurements of the Earth's density are already quite precise [19]. Therefore any neutrino-based measurement must have a comparably small statistical and systematic errors to be competitive. For example, in Ref. [12] IceCube assigned a 1% systematic uncertainty to account for uncertainties in the Earth's mass distribution.
These studies assume that high-energy neutrinos interact with isospin-invariant, unshadowed targets via Deep inelastic Scattering (DIS). The purpose of this paper is to examine how nuclear corrections (compared to an isospin-invariant non-shadowed reference) to DIS cross sections and inelasticity distributions affect the assumptions used in neutrino-telescope analyses. We focus on energies above 1 TeV.
The reference cross sections for these studies are next-to-leading order (NLO) perturbative QCD calculations [20,21]. Because neutrino interactions involve W ± and Z 0 exchange, higher-order corrections to the cross sections should be small and the dominant uncertainties should arise from the uncertainties in the parton distribution functions (PDFs). The uncertainties on these cross sections due to the proton PDFs was found to be 5% for E ν from 50 GeV to 10 9 GeV, rising significantly at higher energies [20].
These calculations have some limitations. While they assume free-nucleon targets, the Earth contains significant contributions from elements as light as hydrogen to as heavy as iron. Shadowing may be significant for the heavier nuclei. The calculations assume that the targets are isoscalar, with equal numbers of protons and neutrons. However, H 2 O is 62.5% protons. Some other elemental components of the Earth are neutron rich. The GENIE model [22], which is commonly used to model neutrino interactions at energies up to a few hundred GeV, includes non-isoscalar targets and some other nuclear effects [22]. At these lower energies, neutrino interaction phenomenology is quite rich [23]. This paper will consider how shadowing and non-isoscalar targets affect neutrino absorp-tion in the Earth and the observed inelasticity distribution of interactions in ice. We explore the change to neutrino absorption for different paths through the Earth, accounting for the different elemental components of the core, mantle and crust; determine the nuclear modifications for the elements in different chords through the Earth (i.e. different zenith angles); and calculate the overall change in cross sections. We use a leading-order pQCD calculation, forming a ratio to a parallel calculation assuming an unshadowed, isoscalar, deuterium target. We do not consider initial-state phenomena such as the colored glass condensate (CGC) [24], or other factors that might lead to larger changes in the cross sections. The presence of a CGC was predicted to greatly reduce the neutrino-nucleon interaction cross sections at energies above ∼ 10 10 GeV [25] because of the large change in the gluon density. The calculation in Ref. [25] considered free nucleons. However, heavier nuclear targets increase the minimum Bjorken-x required for the onset of the CGC regime as A 1/3 [26] where A is the atomic number. In the case of iron, A 1/3 ≈ 3.8. Thus, in addition to reducing the overall neutrino cross section, a CGC initial state would also alter in the zenith-angle dependence of neutrino absorption in the Earth, in particular, a larger reduction in cross section in the Earth's core. Another calculation considered the effects of nonlinear QCD evolution, at sufficiently high parton densities for recombination to be important, and found a smaller, factor-of-2, reduction in the cross section [27] at energies as low as 10 4 TeV. This result seems to exhibit at least mild tension with the IceCube ν µ analysis [12]. A third calculation found a similar reduction at higher energies, above 10 8 TeV [28]. They noted that, in this regime, nonlinearities lead to significant uncertainties in the cross section.
We consider the effects on the cross section and inelasticity distribution for neutrinos interacting in ice or water, particularly the effects on hydrogen with a nucleus of a single proton. We discuss how these nuclear effects, and their uncertainties, affect measurements that can be performed using ice-based neutrino telescopes.

II. CALCULATIONS OF NUCLEAR EFFECTS
Both charged lepton and neutrino interactions in nuclei in deep-inelastic scattering (DIS) measurements show modifications of the quark parton densities in the medium, even for nuclei with masses as low as A = 4. Here we give the neutrino-nucleus cross sections for charged and neutral currents for non-isoscalar targets including nuclear modifications of the parton densities (nPDFs), relative to a deuteron. Although our calculation is leading order (LO), we use NLO nPDFs and PDFs to obtain a first estimate of the effect on absorption of neutrinos by the Earth. There are few strictly LO proton PDFs and those that exist do not include error sets. In addition, the latest nPDF set, EPPS16 [29] only has a NLO set available. We have checked the proton to deuterium ratio at LO calculated with both CTEQ61L (a fully LO calculation) and with CTEQ6M [30] (LO cross sections and NLO proton PDFs) and found agreement between the two calculations on the sub-percent level.
We chose these two sets because they were the CTEQ sets used in the global analyses of the EPS09 nPDF sets at LO and NLO respectively. (We note that the nature of the gluon densities in these sets may lead to some differences for gluon-dominated processes. This would not be the case for other nPDF sets based on proton parton densities with similar small x gluon behavior [31].) This justifies using a LO calculation to test the modifications in question using calculated ratios.
Neglecting the longitudinal structure function, F L , the charged current cross section for neutrino-or antineutrino-proton interactions is where is the square of the momentum transferred from the neutrino to the proton. The fraction of the proton momentum carried by a quark is x = Q 2 /(2m(E ν − E ν )) with proton mass m, incoming and outgoing neutrino energies E ν and E ν respectively, and inelasticity y = The structure functions for the proton, F νCC 2p and xF νCC 3p , for an exchanged W + are Here we suppress the dependence of the structure functions and the individual quark parton densities on x and Q 2 to be more concise. Note that u, d, s, c and b refer to the up, down, strange, charm and bottom quark distributions. The up and down distributions include contributions from the valence quarks that define the proton identity, u V and d V , while the perturbatively-generated sea quark and antiquark distributions can be referred to as q s or q interchangeably. We assume q s = q s = q since these distributions are produced pairwise through gluon splitting and do not contribute to the baryon number or charge of the hadron, unlike the valence quark distributions. These sea distributions dominate the quark PDFs at low momentum. We also define Likewise, for an exchanged W − , the structure functions for charged current interactions of antineutrinos with a proton are xF νCC To go from a proton target to a nuclear target, we have to define the structure functions for charged current neutrino interactions on a neutron. We take u p = d n , d p = u n , u p = d n and d p = u p . We make the distinction for the light quark sea because the u and d distributions were found to be different in studies of Drell-Yan dilepton production in p + p and p + d interactions [32][33][34]. This difference has been incorporated into global analyses of the proton PDFs since the MRS(A) [35] and CTEQ4 [36] sets. Writing the neutron structure functions in terms of the proton parton densities, we have We write the structure functions for neutrino-nucleus interactions with contributions from Z protons and N neutrons in a nucleus of mass number A = Z + N , leaving out, for the moment, the modifications due to the nuclear medium, as For an isoscalar target, Z = N = A/2, we have e.g. F νCC 2A = 0.5(F νCC 2p + F νCC 2n ) = x(u + u + d + d + 2s + 2c + 2b), the usual definition for these calculations. The results for the charged current structure functions with a nuclear target are To include nuclear modifications of the parton densities in the nucleus, we introduce the shadowing ratios, R i (x, Q 2 , A), discussed below. We assume that they are distinct for each quark flavor, as well as differentiating between valence and sea contributions for up and down quarks. The nPDFs are determined from global analyses of experiments with nuclear targets, including Drell-Yan and nuclear DIS with charged leptons and neutrinos (in some cases). The various contributing processes allow separation of effects on the valence and sea distributions. These analyses will be discussed further shortly.
We use the EPS09 [37] and EPPS16 [29] sets at next-to-leading order for the modifications. Thus, in Eq. (1), we replace σ CC (ν(ν)p) by σ CC (ν(ν)A) and the proton structure functions in Eqs. (2) and (3) with the non-isoscalar structure functions including these nuclear modification ratios: We assume above that the sea quark modifications are identical for quarks and antiquarks, We now turn to neutral currents. The cross sections are defined similarly where We have now directly written the cross sections and structure functions in terms of nuclear mass number A. In this case, there are additional u-like and d-like couplings on the structure functions [38]. With F 2 , the u-like couplings a 2 u + v 2 u multiply the up and charm parton densities while the d-like couplings a 2 d + v 2 d multiply the down, strange and bottom parton densities. In the case of xF 3 , the u-like couplings are 2a u v u while the d-like couplings are The subtraction of the identical q and q at the vertices for the F 3 structure function with isoscalar target and no shadowing leads to very simple structure in this case, depending only on the valence quarks with F νNC there is no cancellation when N = Z and the expression is therefore more complex. Including the nuclear modifications, the structure functions for the neutral current We now turn to nuclear shadowing parameterizations. We use the CT10 proton parton densities [39] for the free nucleons in our further calculations. We use the central set only and do not include uncertainties in the proton PDFs here to focus on the effect of the nuclear modifications. We will, however, discuss the sensitivity of the inelasticity distributions to the proton PDF choice in Sec. IV.
Both of the parameterizations of the nuclear modifications used in this work assume collinear factorization with DGLAP evolution. Both sets are optimized assuming that there is no shadowing effect present in deuterium (A = 2, Z = 1). While shadowing may depend on where the probe impacts the nucleus, for example, closer to the 'edge' of the nucleus where there might be only one nucleon in its path or more in the center where it may encounter multiple nucleons [40,41], this is not taken into account. Thus the parameterizations themselves are blind to the nuclear shape and density so a more loosely bound nucleus like 6 Li, which might be described as an alpha particle ( 4 He) with two neutrons is treated the same way as a tightly bound nucleus such as 56 Fe with two closed shells. Or, nuclei with neutron skins might experience enhanced shadowing for protons with correspondingly weaker shadowing for neutrons.
EPS09 [37] defines three different nuclear corrections at the initial scale Q 2 0 = 1.69 GeV 2 : R A V for both up and down valence quarks; R A S for all sea quarks; and R A G for gluons. Fifteen fit parameters were employed, resulting in 30 error sets determined by varying each parameter by one standard deviation in each direction from its optimized value, in addition to the best fit, central set. The nuclear dependence of each of the parameters is assumed to follow Outside of these ranges, the value of the required ratio at the minimum x or maximum Q is returned.
This set, developed before the LHC turned on, relied primarily on fixed-target DIS of electrons and muons from nuclear targets of He, Li, Be, C, Al, Ca, Fe, and Cu measured relative to scattering off deuterium [42][43][44][45][46]. Drell-Yan studies from the Fermilab E772 [47] and E866 [48] experiments produced ratios of C/D, Ca/D, Fe/D (E772) and Fe/Be (E866) that could be used to separate valence from sea contributions. None of these data were significantly above Q 2 = 100 GeV 2 . The Drell-Yan data were primarily in the range 16 < Q 2 < 81 GeV 2 (corresponding to the dimuon mass range of 4 < M < 9 GeV, between the J/ψ and the Υ spectral peaks) and probed 0.01 < x < 0.2 with the precise x range shifting for each mass bin. Only DIS data above the minimum Q 2 used by EPS09, 1.69 GeV 2 , were employed in the analysis, leaving an x range of 0.005 < x < 0.7 available for fitting the modifications. Therefore, one cannot expect great sensitivity to the individual valence and sea quark distributions. Any small differences between the up and down valence and sea distributions in nuclear targets are due to the Q 2 evolution. It is notable that only the valence quark distributions showed broad antishadowing. In contrast, even though they followed the same shape, the up and down sea quark distributions did not produce ratios significantly above unity and then only in the lower part of the antishadowing x range. While the sea quarks all evolve from gluons, the strange and charm sea quarks are less sensitive to the data used in the fit and thus follow the shape of the gluon ratio more directly, resulting in larger antishadowing. Counterintuitively, the charm ratio shows more antishadowing than does the strange quark ratio.
EPPS16 [29], the successor to EPS09, was innovative in several ways. The fit used LHC data from the 2012 p+Pb run at a center of mass energy of 5.02 TeV. The dijet [49] and gauge boson [50][51][52] data sets, while only a few points each, probed much higher Q 2 scales than previously possible. Although the same fixed-target electron and muon DIS data were used for EPPS16, the group also used the extensive CHORUS data [53] with ν and ν beams on a lead target. (Since it also employed a fixed target, CHORUS covered a similar region in the x and Q 2 plane as the charged lepton DIS data: 4 < Q 2 < 100 GeV 2 and 0.05 < x < 0.7). Along with the CMS W ± data [50], these data were sensitive to differences between the valence quark distributions as well as differences between the various sea quark distributions. Although all these newly-incorporated data sets employed a lead target, a much heavier nucleus than relevant here, these data informed the lower A results through a power-law A scaling of the parameters, similar to that used in EPS09.
Because there was more information available to distinguish between the quark distributions, each one was treated separately. The nuclear dependence of the parameters was also handled somewhat differently in EPPS16, adjusting the A dependence to ensure that nuclear effects are larger for heavier A. Due to the greater number of available constraints, the number of parameters increased. EPPS16 has 20 parameters, giving 41 total sets with one central set and 40 error sets. The error sets and the total uncertainty are produced the same way as in EPS09. The sets cover the range 1.3 < Q < 10000 GeV and 10 −7 < x < 1, reaching both lower x and higher Q than EPS09. The EPS09 valence and light sea ratios in Fig. 1 are very similar since they assume that, at the starting scale, . Some differences arise with Q 2 evolution according to the DGLAP equations giving, for example, a more pronounced EMC effect at large x for d ratios than the u. While the strange quark ratio R s (x, Q 2 0 , A) might have started out equal to that of the light quark sea, it evolves to a larger antishadowing than the u and d ratios. The charm quark sea ratios, more closely connected to the gluon ratio, show similar significant antishadowing, as does R g (x, Q 2 , A).
This may seem somewhat counterintuitive, however, because one might expect the heavier quarks, which enter the Q 2 evolution only above the quark mass scale, to be modified less in the presence of a medium. The x range is extended to 10 −7 to show that, for the very low x range, x < 10 −6 , the ratios are fixed to their value at the lowest x considered by EPS09, The effects of including the CHORUS neutrino data as well as the LHC gauge boson data in the EPPS16 sets are clearly illustrated by comparing Fig. 2 to Fig. 1. There are clear differences between the valence ratios and between the sea quark ratios. The nuclear modification is weaker for d V than for u V while there is no antishadowing (quantified as a ratio larger than unity) at all for d in this x range. On the other hand, the u ratio shows antishadowing on a level similar to that of the valence quarks, albeit over a narrower range in x. Because the level of antishadowing in the strange quark ratio was allowed to float, it shows stronger antishadowing than the charm ratio, closer to what might be intuitively expected.
As E ν increases, the cross sections probe successively lower values of x, with E ν > 10 5 GeV corresponding typically to x < 0.01, marking the transition from the antishadowing to the shadowing region where a suppression relative to deuterium should be observed. The  subtracted, as for antineutrinos. In the case of CC interactions, the differences between ν and ν should be larger, especially for non-isoscalar targets, because of the difference in valence quark content between protons and neutrons. which are better constrained in the global analyses of the nPDFs.
All four cases show a decrease in R with increasing E ν . The fairly sharp drop between 100 and 500 TeV is due to the transition from the antishadowing region to the shadowing region. At small E ν , there is considerable variation in shadowing for E ν < 500 TeV for CC interactions, while for NC interactions, the nuclear effects are much smaller, as we now discuss.
The charged current results display a clear separation between isoscalar and non-isoscalar targets. The ratios for the isoscalar targets (He, Li, C, O, Ne, Si, and Ca) in charged current interactions are all somewhat higher than unity for E ν ≤ 10 5 GeV, for both ν and ν.
However, non-isoscalar targets (Be, V, Fe and Cu) are different, with a larger enhancement for ν and a slight suppression for ν because N > Z. From Eqs. (19)- (22), we see that the valence distributions are weighted differently in CC interactions. In neutrino-initiated CC interactions, Z multiplies d V while N multiplies u V . This is opposite for antineutrinoinduced CC interactions. When Z = N , there is no difference but one arises when N > Z.
In the case of NC interactions, Eqs. (27) and (27) At low E ν the effects of isoscalar targets dominate over the effects of shadowing. At higher E ν , the enhancement is washed out by larger contributions from sea quarks. When E ν ≤ 10 5 GeV, there should be a stronger dependence on isospin for more neutron-rich nuclei such as lead. These results also demonstrate a non-negligible effect from the inclusion of nPDFs, beyond the difference between isoscalar and non-isoscalar targets.

III. NUCLEAR EFFECTS ON CROSS SECTIONS
Neutrino telescopes measure the neutrino interaction cross section by observing neutrino absorption in the Earth [54] as a function of neutrino energy E ν and zenith angle θ; spherical symmetry is assumed. The absorption measurements can be unfolded to obtain the distri- This corresponds to a zenith angle θ = 180 • , while horizontal incidence is θ = 90 • and θ = 0 • corresponds to vertically-downward incidence. Figure 5 shows how the zenith angle corresponding to a chord length equal to 1Λ decreases with increasing neutrino energy.
The logarithm of the angular distance from the horizon (θ horizon = 90 • ), log(θ 1Λ − θ horizon ), is shown. Large zenith angles are required to measure absorption at TeV-scale energies.
As the energy rises, more horizontal angles become more important. The Earth is almost opaque to neutrinos of energies ∼ 10 9 GeV and experiments are most sensitive to absorption near the horizon, primarily in the mantle and crust. At these energies, the Glashow resonance dominates the ν e cross section, complicating DIS measurements.
To quantify the nuclear effects on neutrino absorption along different paths, we first determine the cross sections without nuclear effects. We then calculate the same probability with nuclear effects, using EPPS16 shadowing. We quantify the effect by R, the ratio of the probabilities with and without shadowing.
We use a simplified model of the Earth, based on Table 2.2 of Ref. [56], which gives the elemental composition of the Earth in four regions: crust, mantle, outer core and inner core.
We use a single core region, based on an average of the inner and outer core components.  An average R, R region , is calculated in each region based on the percentage abundance of each element and using an interpolation for elements not included in EPPS16.
A neutrino traversing the Earth at zenith angle θ has a path length L through the three layers (crust, mantle and core), as shown in Fig. 6. Nuclear shadowing effects for each zenith angle are based on the average over the whole path length, where R region in each region is weighted by its percentage of the total path length.
We consider CC and NC interactions separately. CC interactions are straightforward: the neutrino interacts and disappears. NC interactions are more complicated because the neutrino loses a fraction of its energy without disappearing. The average fractional energy loss is 20% at high neutrino energies. The effect of this energy loss can be handled by treating absorption as a two dimensional matrix, relating the energy distribution of neutrinos entering the Earth to the energy distribution of detected neutrinos. This is expressed as a relationship between two fluxes, Φ in a function of E ν when the neutrino enters the Earth, and Φ out , a At angles above 147 • the path length transverses the Earth's core, resulting in an observable effect on neutrino interactions, as highlighted in Figure 5.
function of the detected neutrino energy [12]. The result can be expressed in terms of an apparent transparency, T , at a single energy E ν , T = Φ out (E ν , θ)/Φ in (E ν , θ), but, because of energy loss by more energetic neutrinos, this apparent absorption at a given E ν has some dependence on the assumed neutrino spectrum. We avoid this spectral dependence here by considering only monoenergetic neutrinos. Figure 7 shows R as a function of neutrino energy and zenith angle, for ν µ and ν µ . The shift from shadowing to antishadowing is visible at E ν ∼ 500 TeV. The effect of the Earth's core, with its high iron and likely nickel content, enhances nuclear effects for large zenith angles. These results are shown more quantitatively in Fig. 8 where R is presented as a function of θ for several discrete neutrino energies. R drops from about 1.02 to 0.96 for increasing neutrino energies. With trajectories traversing the Earth's core, the spread is slightly larger due to the greater nuclear effects. The differences between ν and ν are shown to be relatively small.
The probability for a neutrino to survive passage through the Earth is where L is the path length through the Earth for a given zenith angle θ while Λ(E ν ) ∝ 1/σ(E ν ) is the absorption length. Nuclear effects modify the cross section so that σ(E ν ) ∼ Rσ 0 (E ν ) where σ 0 (E ν ) is the cross section for an isoscalar target without shadowing. At an energy and zenith angle where L ≈ Λ(E ν ), the change survival probability due to the inclusion of nuclear effects, ∆P (E ν ) is roughly comparable to −(R − 1). As L rises, dP (E ν )/dR ≈ −L/Λ so that for long chords through the Earth with a high absorption probability, L/Λ 1, small changes in R lead to larger changes in the absorption probability.
The larger modification of P (E ν ) with R due to shadowing could affect the interpretation of anomalous events, such as the possible ν τ events observed by ANITA [57]. For example, ANITA-III event 15717147 traverses a 7,000 km chord through the Earth, corresponding to ∼ 15Λ without nuclear effects, requiring contributions from BSM physics. Reducing the cross section by 4% due to shadowing roughly doubles the survival probability. While this reduction is too small to alter the overall conclusion, a larger reduction in cross section, such as from a colored glass condensate [25] might allow this event to be interpreted without requiring BSM physics.
It is important to note that, although we have focused on the best-fit values of EPS09 and EPPS16, the uncertainties are significant. The uncertainties on these R values are considerable, more than ±5 % and thus larger than |R − 1|. This uncertainty will limit many quantitative analyses.

IV. NUCLEAR EFFECTS ON INELASTICITY
The inelasticity of high-energy ν µ and ν µ charged-current interactions on ice targets has been measured by the IceCube Collaboration [6]. They separately measured the energy of the hadronic cascade from the struck nucleus, E cas , and the energy of the outgoing muon, E µ , and calculated the inelasticity as y = E cas /(E µ + E cas ). There are two major reconstruction challenges to measuring neutrino inelasticity in a neutrino telescope. First, a muon with energy above 1 TeV travels many kilometers before losing its energy, making the energy determination difficult. Second, there are corrections to account for missing (carried off by the neutrino) and mismeasured (hadronic instead of electromagnetic) energy.
Measurements using starting tracks, defined as a track originating in the active volume of the detector, have many applications. In addition to inelasticity, starting tracks are important for flavor ratio measurements; determining ν/ν from the difference in inelasticity distributions; and searching for charm production in neutrino interactions [58]. Inelasticity is also relevant for ν τ decays since it helps determine the division of energy between the ν τ interaction and the τ , i.e. the energy division between the first and second showers in double-bang (ν τ events with E ν > 10 6 GeV, characterized by a hadronic shower originated when the ν τ interacts, a ∼ 100 m minimum ionizing track, and a subsequent cascade when the τ decays) and double-pulse (similar to double-bang but at lower energy, resulting a ∼ 10 m minimum ionizing track) events. Changes in the inelasticity distribution may also affect the relationship between E ν and muon energy in through-going muon events.
Ice consists of H 2 O. As we discuss, hydrogen and oxygen exhibit different nuclear effects.
While oxygen is subject to nuclear shadowing, the nuclear effects are dominated by hydrogen because it strongly violates isospin invariance. As Figs. 9 and 10 show, hydrogen has a large |R − 1|, with a rather distinctive energy and inelasticity dependence. The inelasticity can be related to x and Q 2 by where m is the nucleon mass. The inelasticity at a given y and fixed E ν is integrated over x and Q 2 . The integral over Q 2 starts from a minimum Q 2 of 9 GeV 2 . The cross sections increase with Q 2 so that at higher E ν the Q 2 integral is dominated by Q 2 ≈ M 2 W . Figures 9 and 10 show R relative to an isospin invariant deuterium target for both protons and H 2 O, using the CT10 and CT14 [59] proton PDFs respectively with the EPPS16 nPDFs.
Results are given for values of E ν separated by an order of magnitude for E ν from 10 2 to 10 7 GeV.
In all cases, as y decreases, x increases. The shape of R for fixed E ν is balance of contributions from x and Q 2 . The Q 2 range at a given y is constrained by the requirement that x < 1. For a given E ν , as y → 1, Q 2 can be large, near the maximum of the available range, and x will remain less than unity. On the other hand, when y → 0, the denominator of Eq. (29) is small so that Q 2 nust be near its minimum value for x to be less than unity. At large y, x is small and decreases with increasing E ν . Thus, in this range, the sea quark distributions dominate R and R approaches unity, both as y → 1 and E ν → ∞. In models with a colored glass condensate, or strong nuclear shadowing, the cross section in this region would be reduced.
As y → 0, on the other hand, the valence distributions dominate R and, here, a difference between the proton PDFs appears at y < 0.1. The change in the slope of R in this region is due to the behavior of the d valence distribution of the CT10 PDFs. The origin of this rather abrupt change in slope in CT10 is not clear but it is absent in the CT14 sets. We have checked the older proton PDF sets, CTEQ6M and GRV98 [60], as well and the calculated R values for hydrogen with these sets agree with the CT14 results.
The bottom panels of the figures show R for H 2 O. In this case R is dominated by the effect on oxygen. As a consequence of the interplay of x and Q 2 as a function of y discussed earlier, the inelasticity curves in these figures trace the inverse of the EPPS16 shadowing ratios in Fig. 2. The rise at low y visible for E ν = 100 GeV is the result of the high x behavior of the shadowing function; the decrease to the minimum of R corresponds to the EMC region and the subsequent rise at large y is the effect of antishadowing. As E ν increases, the EMC and antishadowing regions become more compressed at low y and a decease is seen instead of a rise at larger y because x is in the shadowing region. Although the effect is relatively independent of whether the interaction is initiated by ν or ν, there is a slightly stronger rise for ν-initiated interactions because of the increased antishadowing in the CC interactions. In addition, as seen in Fig. 4, the u = u V + u distribution dominating antineutrino interactions includes antishadowing on both the valence and sea distributions in this region while the d = d V + d dominant for neutrino interactions shows effects of antishadowing on the valence distribution but shadowing in the same x region for the d distribution. Finally, we note that the rather large uncertainties visible in both cases for E ν = 100 GeV arise here because, even though higher Q 2 values are probed at large y, large uncertainties remain due to the 40 error sets for EPPS16.
The uncertainties on R due to the nPDF uncertainties are much larger than for the proton PDFs, on the order of 20%. Because the uncertainties on the nPDFs decrease with Q 2 , as does the overall nuclear modification, the largest uncertainties are on the smallest E ν values shown for fixed E ν . For fixed y, the largest nPDF uncertainties are at small y because the smallest Q 2 range is probed here.
Since the effect on oxygen dominates R for H 2 O, the difference between the CT10 and CT14 PDF sets does not play a significant role here. The relative independence of the nuclear effects on the proton PDFs seen here shows that the choice of proton PDFs did not play a role in the cross section effects discussed in the previous section. Larger nuclear effects, due to color glass condensates or nonlinear parton dynamics, could also alter the inelasticity distribution [27]. Since x and y are inversely related, reductions in the cross section are likely to be most prominent at large y, corresponding to small x.
These nuclear effects have implications for neutrinos and antineutrinos with energies from 100 GeV to a few TeV. At larger energies, the central values of R − 1 are fairly small but with significant uncertainties due to shadowing.
There are two key points arising from the result that R − 1 = 0: the nuclear effects are large at small inelasticities and they behave differently for ν and ν, particularly for hydrogen.
The differences in ν and ν may cancel for some effects if there is an equal admixture of ν and ν. At E ν ∼ 1 TeV, the flux of atmospheric neutrinos (the dominant source) is expected to be in the proportion ν/ν ∼ 1.55, rising to 1.75 at 100 TeV [12], resulting in an incomplete cancellation. We will consider the implications of these nuclear effects on a number of different analyses.
IceCube found rather good agreement between their data and inelasticity expectations [12]. In their lowest energy bin, 1 < E ν < 3 TeV, the rise in R led to an increase in low inelasticity events. IceCube has limited low inelasticity acceptance in this energy region because of the requirement that the number of detected photons, summed over the IceCube volume, N pe , is large enough to observe the event. The absence of a significant hadronic shower means that most of the neutrino energy is transferred to the muon and escapes the detector. There is room in these data for nuclear effects.
Many other neutrino telescope measurements rely on an implicit knowledge of the inelasticity distribution because many starting-track analyses have an inelasticity-dependent inefficiency. These analyses also require a minimum N pe , or other similar criteria, making them insensitive to events with y ≈ 0.
Inelasticity measurements have been used to measure ν/ν. At energies below about 10 TeV, ν and ν have different inelasticity distributions because of the contribution from the valence quarks. At these energies, the flux is dominated by atmospheric ν which then dominate the ν/ν measurement. However, with a low-threshold surface air shower detector to veto neutrinos accompanied by an air shower, it might be possible to eliminate most downgoing ν from the atmosphere. Because ν and ν experience different nuclear modifications, measurements of ν/ν may require a correction for nuclear effects.
In concert with cascades, starting tracks are also used to determine the neutrino flavor ratio. High inelasticity starting tracks will be mistaken for cascades while low inelasticity starting tracks will be missed because of the N pe cut.
Inelasticity also affects the energy that incident ν τ transfer to nuclear targets, affecting the energy seen in the first bang of double-bang events.
The excess of low-inelasticity events could be an unexpected background in searches for electromagnetic neutrino interactions which do not involve a nuclear recoil [61,62] and may be relevant in searches for new, BSM interaction topologies.

V. CONCLUSIONS
We have examined the role of two nuclear effects: shadowing and violation of isospin invariance in ultra-high energy neutrino interactions in both the Earth and in polar ice packs.
Antishadowing decreases the neutrino cross section by about 2% at energies below about 500 TeV while shadowing increases it by about 4% at higher energies. These corrections should be included in new high-accuracy measurements of neutrino absorption in the Earth, including cross section measurements and Earth tomography. These estimated corrections are based on standard QCD evoluion of the nPDFs; new phenomena like the colored glass condensate could perhaps lead to larger modifications.
Nuclear effects in both hydrogen and oxygen affect the inelasticity distribution of neutrino interactions in ice, particularly at low inelasticity, in an energy range that probes quarks at large Bjorken-x. The Fermi motion of nucleons in oxygen increases the number of highmomentum quarks and thus the cross section at very low y. On the other hand, the difference in sign of the nuclear effects for ν and ν on hydrogen will affect measurements of the ratio ν/ν. These effects are largest for E ν from 100 GeV up to a few TeV. At higher energies, the modifications are shifted to very low y and do not contribute greatly.
Although these central shifts are relatively small, the uncertainties are quite significant.
Until the uncertainties are reduced, they will limit the precision of many measurements. The uncertainties rise with neutrino energy and are particularly large at the energies targeted by next-generation radio-detection systems [63].
With the advent of higher-precision neutrino measurements [64], it will be necessary to account for nuclear effects in cross section and tomography measurements. This need will only increase in larger, next-generation neutrino observatories [65][66][67]. The introduction of large (order 100 km 3 ) radio-pulse-based neutrino telescopes, with energy thresholds above 10 7 GeV will also open the door to cross section measurements at higher energies [63,68], provided that good angular resolution is achieved for near-horizontal events. These experiments may also be able to measure inelasticity in ν e interactions at very high energies by taking advantage of the LPM effect, which elongates electromagnetic showers, allowing them to be separated from the hadronic shower on the basis of Cherenkov cone widths [69] or via multiple showers [6].