Measurements using the inelasticity distribution of multi-TeV neutrino interactions in IceCube

Inelasticity, the fraction of a neutrino ’ s energy transferred to hadrons, is a quantity of interest in the study of astrophysical and atmospheric neutrino interactions at multi-TeVenergies with IceCube. In this work, a sample of contained neutrino interactions in IceCube is obtained from five years of data and classified as 2650 tracks and 965 cascades. Tracks arise predominantly from charged-current ν μ interactions, and we demonstrate that we can reconstruct their energy and inelasticity. The inelasticity distribution is found to be consistent with the calculation of Cooper-Sarkar et al. across the energy range from ∼ 1 to ∼ 100 TeV. Along with cascades from neutrinos of all flavors, we also perform a fit over the energy, zenith angle, and inelasticity distribution to characterize the flux of astrophysical and atmospheric neutrinos. The energy spectrum of diffuse astrophysical neutrinos is described well by a power law in both track and cascade samples, and a best-fit index γ ¼ 2 . 62 (cid:2) 0 . 07 is found in the energy range from 3.5 TeV to 2.6 PeV. Limits are set on the astrophysical flavor composition and are compatible with a ratio of ð 13 ∶ 13 ∶ 13 Þ ⊕ . Exploiting the distinct inelasticity distribution of ν μ and ¯ ν μ interactions, the atmospheric ν μ to ¯ ν μ flux ratio in the energy range from 770 GeV to 21 TeV is found to be 0 . 77 þ 0 . 44 − 0 . 25 times the calculation by Honda et al. Lastly, the inelasticity distribution is also sensitive to neutrino charged-current charm production. The data are consistent with a leading-order calculation, with zero charm production excluded at 91% confidence level. Future analyses of inelasticity distributions may probe new physics that affects neutrino interactions both in and beyond the Standard Model.

the flux of astrophysical and atmospheric neutrinos.The energy spectrum of diffuse astrophysical neutrinos is well-described by a power-law in both track and cascade samples, and a best-fit index γ = 2.62 ± 0.07 is found in the energy range from 3.5 TeV to 2.6 PeV.Limits are set on the astrophysical flavor composition that are compatible with a ratio of 1  3 : 1 3 : 1 3 ⊕ .Exploiting the distinct inelasticity distribution of νµ and νµ interactions, the atmospheric νµ to νµ flux ratio in the energy range from 770 GeV to 21 TeV is found to be 0.77 +0.44  −0.25 times the calculation by Honda et al.Lastly, the inelasticity distribution is also sensitive to neutrino charged-current charm production.The data are consistent with a leading-order calculation, with zero charm production excluded at 91% confidence level.Future analyses of inelasticity distributions may probe new physics that affects neutrino interactions both in and beyond the Standard Model.

I. INTRODUCTION
The observation of astrophysical neutrinos [1,2] was a landmark in high-energy astrophysics.It introduced a new probe that is directionally sensitive to high-energy hadronic particle accelerators in the universe.Neutrinos provide both good directional information, unaffected by magnetic fields, and extremely long range, allowing us to probe accelerators at cosmologically significant distances.Measurements of the flux and energy spectrum [3,4] and flavor composition [5,6] are so far fully compatible with conventional acceleration models, but more exotic production mechanisms cannot be ruled out.
At the same time, the observation of high-energy astrophysical and atmospheric neutrinos by detectors like IceCube has opened up the study of neutrino interactions at energies orders of magnitude above those accessible at terrestrial accelerators.Already, the 1 km 3 IceCube neutrino observatory has used atmospheric and astrophysical neutrinos to measure neutrino absorption in the Earth, and from that determined the neutrino-nucleon cross section at energies from 6.3 TeV to 980 TeV to be in agreeement with the Standard Model prediction [7].
In this paper, we report on a new study of high-energy charged-current (CC) ν µ interactions contained within IceCube's instrumented volume.These interactions produce a cascade of hadrons and a muon, an event topology known as a starting track.By estimating the hadronic cascade and muon energies separately, we can estimate the inelasticity of each interaction -the ratio of hadronic cascade energy to total neutrino energy [8].The central 90% of neutrinos have estimated energies in the range from 1.1 TeV to 38 TeV, energies far beyond the reach of terrestrial accelerators.For example, the NuTeV data were used to measure inelasticity distributions at energies up to 250 GeV [9], while earlier experiments were limited to lower energies [10].
The starting track data, together with a similarly obtained set of cascades due to all neutrino flavors are binned by reconstructed energy and zenith angle and (for the tracks) inelasticity.This data is fitted to a neutrino flux model containing both atmospheric and as-trophysical neutrinos.From this, we present measurements of neutrino inelasticity and estimate the fraction of neutrino interactions that produce charmed particles.We also compare the track and cascade samples, study whether they have the same astrophysical neutrino flux spectral indices, and constrain the flavor composition of astrophysical neutrinos.Finally, we measure the ratio of neutrinos to antineutrinos in the atmospheric neutrino flux.

II. ν INTERACTIONS AND INELASTICITY
The inelasticity distribution of neutrinos is expected to be well described by the Standard Model for weak interactions.At TeV energies, the interactions are dominated by Deep Inelastic Scattering (DIS).Neutrinos interact with quarks in nuclear targets (hydrogen and oxygen) via charged-current (CC) and neutral-current (NC) reactions.In NC interactions, the neutrino interacts via Z 0 exchange, leaving the quark flavor unchanged.In CC ν interactions via W ± exchange, the quark charge changes by one, turning a charge +2/3 quark into a charge −1/3 one, while for ν, the reverse reactions occur, but with a different inelasticity distribution.W ± also interact with sea quarks and antiquarks in the nucleus in similar charge-changing reactions, so the differences between neutrinos and antineutrinos largely disappear at very high energies.The relative probability for producing a given final state quark depends on the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, so reactions involving u → d and s → c predominate.
The neutrino-nucleon cross section involves four main kinematic variables: s is the neutrino-nucleon center of mass energy squared, Q 2 is the negative of the 4momentum transfer squared, Bjorken−x, the fraction of the struck nucleon's momentum carried by the struck quark, and the inelasticity, y.These quantities are related through x = Q 2 /2M E ν y, where M is the nucleon mass and E ν is the incoming neutrino energy.Most of the interactions discussed in this paper involve struck partons with 10 −2 < x < 10 −1 , and the typical Q 2 is ∼ M 2 W,Z ≈ 6 × 10 3 GeV 2 .The double-differential cross section for a neutrino with energy E ν is where the + sign is for neutrinos, and the − is for antineutrinos.M V is the vector boson mass (M W for charged current interactions and M Z for neutral current interactions) and G F is the Fermi coupling constant [11].The cross sections depend on three nucleon structure functions: F 2 , F 3 , and F L .
Assuming standard vector minus axial vector (V − A) coupling and an isoscalar target, at leading order these structure functions are related to the quark parton distribution functions (PDFs) of the target nucleon, q i (x, Q 2 ), according to [11] The longitudinal structure function F L becomes non-zero in higher order calculations.It should be noted that measurements of the nuclear structure function inferred from neutrino-DIS on an iron target may be slightly different from those observed in charged lepton DIS at low Q 2 [12].This paper uses as a baseline a next-to-leading order calculation ('CSMS') which uses the DGLAP formalism for parton evolution [13].The calculation uses the HER-APDF1.5[14] parton distribution functions, with the MSTW2008 and CT10 [15] distributions used as a standard for comparison.In the relevant energy range, the calculation has an expected uncertainty of a few percent and is consistent with an independent calculation that used the MSTW2008 PDFs [16].
Equation 1 does not account for threshold behaviors which are important for heavy quark production.For this analysis, charm production is most important.In the relevant energy range, charm production is about 10% of the total production.Because of the quark mass, heavy quark production near its production threshold occurs at a larger average inelasticity than for light quarks at the same energy [17].In the case of charm quarks, the inelasticity also tends to be higher even at energies far above the production threshold since they originate primarily from scattering off sea s-quarks.In addition, heavy quarks may decay semileptonically, transferring some of their energy to a muon.This muon will not be separately visible in IceCube.Instead, it will increase the apparent brightness of the primary muon, leading to a higher measured muon energy and lower measured inelasticity.However, due to the small ∼ 10% branching ratio of muonic charm decays, IceCube is primarily sensitive to the signature of larger inelasticity rather than the di-muon signal.
The CSMS calculation is for neutrino DIS on a single nucleon target; there are a few other contributions to consider.Water contains hydrogen, so is not a perfectly isoscalar target.The MINERvA collaboration has seen that nuclear shadowing reduces the cross section for neutrino-heavy ion interactions with x < 0.1 [19].The suppression is expected to increase with decreasing x values [20].However, oxygen is a small nucleus, and we expect the cross section reduction to be small for the moderate x values probed here.
An additional contribution to the cross section is due to neutrino electromagnetic (diffractive) interactions with the Coulomb field of the nuclei, but this is small for low-Z nuclei like hydrogen and oxygen [21,22].These effects are not expected to be significant for this analysis.

III. NEUTRINO SOURCES AND SIMULATION
This analysis uses atmospheric and astrophysical neutrinos as sources.The signals observed in IceCube depend on the neutrino fluxes incident on the Earth, their absorption in the Earth, their interactions in and (for backgrounds) around IceCube, their propagation through the detector and the detector response.We will briefly discuss these factors, with a special emphasis on the interactions, which determine the inelasticity.
Since most of the results in this paper are based on comparisons of data with various simulations, this section will focus on the physics models used in the simulations.These models are implicit in the results presented, and the systematic errors depend on the assumptions used in the simulations; these uncertainties will be discussed when the individual results are presented.

A. Neutrino Sources
Conventional atmospheric neutrinos come from pions and kaons that are produced in cosmic-ray air showers.At the energies relevant for this analysis, the flux roughly follows a power-law spectrum dN/dE ν ∝ E −(γ+1) , where γ is the spectral index for the cosmic-ray energy spectrum.Below the cosmic-ray knee, γ ≈ 2.7, while at higher energies γ ≈ 3.0.The flux is highest for nearhorizontal incidence.This analysis uses the HKKMS [23] flux calculations extrapolated upward in energy, with a modification to account for the knee of the cosmic-ray spectrum [24].This calculation is in good agreement with previous IceCube measurements [25][26][27].
Prompt atmospheric neutrinos come from the decay of charmed mesons produced in cosmic-ray air showers.
They have yet to be observed but are expected to have a hard spectrum that follows that of the primary cosmicrays.This analysis uses as a baseline the BERSS [28] perturbative QCD calculation of the prompt flux, which is tied to recent data from RHIC and the LHC, and consistent with similar independent calculations [29,30].
The number of observed prompt and atmospheric neutrinos requires an important adjustment to account for the IceCube 'self-veto' -a downward-going atmospheric neutrino will be accompanied by an air shower and muon bundle, which may overshadow the neutrino and cause the event to fail the event containment cuts, and so not register as a starting event.This probability for a selfveto depends on whether the neutrino is prompt or conventional and on the neutrino energy and zenith angle.This analysis uses the probabilities calculated in Ref. [31].The muon threshold is taken to be 100 GeV; this is the minimum muon energy which is likely to trigger the selfveto.The appropriateness of the 100 GeV threshold was verified using detector simulations using the CORSIKA program [32].
Previous IceCube measurements have found that, in this energy range, the astrophysical flux is consistent with a single power law.Our fit is based on this single power law.
The neutrinos observed in IceCube may pass through the Earth before reaching the detector; the absorption in the Earth is simulated following the Standard Model cross sections [13].The Earth's density profile is assumed to follow the Preliminary Reference Earth Model (PREM) [33].

B. Neutrino Interactions and Cherenkov light emission
Neutrino interactions in and around the detector are modeled following the CSMS calculation, as described in the previous section.For NC interactions, the inelasticity is the fraction of the neutrino energy that is transferred to the struck nucleus; the remaining energy escapes from the detector.For CC interactions, the remainder of the energy is transferred to a charged lepton.IceCube observes the Cherenkov light emitted by the lepton and its secondary relativistic charged particles.The hadronic cascade from the struck nucleus also produces Cherenkov light.Each type of lepton produces light very differently in the detector.
Electrons produce electromagnetic cascades; the light output is proportional to the electron energy, with the relationship determined from detailed simulations [34].At low energies, they are treated as point sources, while at energies above 1 TeV, the longitudinal profile of a cascade is approximated using a sequence of uniformly spaced point sources.At higher energies above 1 PeV, the longitudinal profile includes the LPM effect [35].
As they traverse the detector, muons radiate energy via ionization, bremsstrahlung, direct pair production and photonuclear interactions; these are modeled following Refs.[36,37].
Tau leptons are propagated through the detector in a manner similar to muons, with adjustments for their higher mass.IceCube simulations allow the taus to decay into ν τ eν e , ν τ µν µ or ν τ plus hadrons.For the leptonic decays, the leptons are propagated through the detector starting at the τ decay point, while the neutrinos escape.For the hadronic decays, the hadronic energy is summed, producing another hadronic cascade, at the point where the τ decays.Because of the energy carried off by the escaping neutrinos, ν τ will deposit less energy in the detector and appear similar to lower energy ν e and ν µ .In muonic decays, the outgoing muon has a lower average energy than in corresponding ν µ interactions.This can affect the measured inelasticity distribution [38].
This analysis is sensitive to charm production in neutrino interactions.DIS interactions produce charm quarks, which then hadronize, forming charmed hadrons and baryons, with lifetimes of order 10 −12 s.We use the hadronization fractions into D + , D 0 , D + s , and Λ + c from the calculation in Ref. [39] at 100 GeV, which is also used in the GENIE event generator [40].If the charmed hadrons have energies above about 10 TeV, they may interact in the ice and lose energy before they decay.The interaction probability depends on the individual hadron-ice cross sections and hadron lifetimes.This analysis used parametrizations of the charm cross sections in CORSIKA [32].Since the parametrizations are only valid above 1 PeV, we extrapolate downward in energy using the kaon-nucleon and nucleon-nucleon cross sections scaled by 88% and 122%, respectively, to match CORSIKA parametrization at 1 PeV for charm mesons and baryons.This leads to critical energies for the D + , D 0 , D + s and Λ + c of 22 TeV, 53 TeV, 47 TeV and 104 TeV, respectively.When the charmed hadrons interact, they lose energy.We use the approach of Ref. [18] to parametrize their energy loss distribution.The observable effect of multiple charm interactions is the production of a low energy muon after a semileptonic decay, mimicking a high inelasticity track for ν e or ν τ interactions.However, due to the 10% muonic branching ratio, this is not a large effect, and the approximate treatment of charm interactions here does not significantly influence later results on charm production.
The conversion between hadronic cascade energy and light follows Ref. [41].Hadronic cascades produce less light than electromagnetic cascades of the same energy, with larger cascade-to-cascade variation.At 100 TeV, a hadronic cascade produces an average of 89% of the light of an equivalent energy electromagnetic cascade.The difference drops with increasing energy.Since electromagnetic and hadronic cascades cannot be readily distinguished in IceCube, we will refer to the visible cascade energy as the energy of an electromagnetic cascade producing an equivalent amount of light as the hadronic cascade.

C. Optical transmission through Antarctic ice
The emitted Cherenkov light travels through the ice, where it may scatter or be absorbed before reaching an optical sensor.Because the sensor array is sparse, only a tiny fraction of the produced Cherenkov photons are observed, and they are likely to scatter multiple times before reaching a sensor.Because of this, the signals seen by IceCube are sensitive to the optical properties of the ice.The scattering and absorption lengths in the ice depend on the position in the ice (largely, but not entirely on the depth below the surface) and on the direction of photon propagation [42,43].These optical properties have been measured using a variety of means, including laser and light emitting diode (LED) signals and cosmicray muons.The optical properties of the ice vary strongly with depth, with certain depth ranges containing "dust layers" with very large absorption and scattering, and other depths providing good transmission.This positional dependence is accounted for with an ice model.
This analysis used as a baseline the "SPICE Mie" ice model [42].It is based on measurements of the ice properties using LEDs mounted in the detector housings, supplemented with parameterizations to give the wavelength dependence of the scattering and absorption lengths.It divides the ice into 10 m thick layers and determines the scattering and absorption lengths separately for each depth range.It does not account for tilts of the dust layers (i.e.variation of optical properties depending on horizontal position), or for anisotropy in the ice.
Near the optical sensors, there is additional scattering in the 'hole ice' -the refrozen column of melted water created by the drill during deployment.This ice contains a central column of bubbles [44].In our baseline simulations, it is treated as having a scattering length of 50 cm.

IV. DETECTOR AND DATA
The Cherenkov light is detected with an array of 5,160 digital optical modules (DOMs), spread over 1 km 3 [45,46].The DOMs are deployed in 86 vertical strings, each containing 60 DOMs.Seventy-eight of the strings are laid out on a 125 m triangular grid.On those strings, the DOMs are deployed with 17 m vertical spacing between 1450 and 2450 m below the surface.The remaining strings are laid out in the middle of the array, with smaller string-to-string spacing [47].On those strings, the DOMs are deployed closer together at depths between about 2,000 and 2450 m.
Each DOM collects data independently, sending digitized data to the surface [48,49].The optical sensor is a 10 inch photomultiplier tube (PMT) [50].The PMT is read out with a data-acquisition system comprising two waveform digitizer systems which are triggered by a discriminator with a threshold of about 1/4 of a typical photoelectron pulse height.One records data with 14-bit dynamic range at 300 megasamples/s for 400 nsec, and the other collects data with 10-bit dynamic range at 40 megasamples/s for 6.4 µs.
All of the digitized data is sent to the surface where a trigger system monitors the incoming data and creates an event when certain conditions are satisfied.For this analysis, the main trigger required 8 hits within a sliding 5 µs window.When a trigger occurs, all of the data within 4 µs before or 6 µs after the trigger time is saved and sent to an on-line computer farm for further processing.
The farm applies a number of different selection algorithms to each event.Each algorithm tests for different classes of interesting events, albeit with significant overlap.This analysis used as input all events that passed either the cascade or muon track filters.These filters have very loose cuts.Simulation studies show that they capture more than 99.5% of the events that pass the other cuts that are applied here.
This analysis uses data collected between May 2011 and May 2016, a total live time of 1734 days during which the detector was in its complete 86-string configuration.During the design of the analysis, 10% of the data was used for testing and the remaining 90% was kept blind.

V. EVENT SELECTION
The analysis aims to select a high-purity sample of neutrino interactions that occur within the detector.Starting tracks are of the greatest interest here, but cascades are more numerous for an astrophysical flux with equal flavor ratio and provide important constraints for many of the fits discussed here.For both channels, the first step of the analysis is to select a clean sample of starting events.The initial cut selects events with more than 100 observed photoelectrons (PE), which captures most contained neutrino interactions above 1 TeV.
The next selection applies an outer layer veto to reject events that come from charged particles that enter from outside the detector.The veto uses DOMs on the top, bottom and sides of the detector.It also includes DOMs in a horizontal layer that passes through the detector, below the "dust layer" that allows charged particles to sneak in undetected.The veto is similar to that in Ref. [5] but with some changes to reduce the energy threshold, as discussed in Ref. [51,52].
First, the number of photoelectrons, Q start , required in a rolling 3 µs time window to define the start of the event, is chosen to depend on the total number of photoelectrons (NPE) of the event, Q tot , The horizontal veto layer that is placed just below the dust layer is made thicker at 120 m.The veto layer at the bottom is expanded to include the lowest DOM on each string.This leaves one remaining significant background, from two nearly-coincident air showers.Sometimes, a low-energy muon can sneak into the detector, producing only a small signal in the veto.Then, a higher energy muon can enter, producing a high enough number of photoelectrons to satisfy the Q tot requirement.To avoid this, an algorithm is used to count the number of causally disconnected clusters of light in the detector and reject events where more than one cluster is found.These cuts reduce the event rate to 0.36 Hz in the 10% testing sample, compared to an expected starting neutrino rate of 0.20 mHz.

A. BDT background rejection
The remaining background is harder to reject, particularly at lower muon energies where less light is produced.To further clean up the signal, a boosted decision tree (BDT) is used to classify the passing events as atmospheric muons, starting tracks or cascades.The BDT uses 15 variables as input, including log 10 (Q tot ) and the number of photoelectrons in the outer layer veto.The other variables are from the output of the track and cascade reconstructions.For the wrong hypothesis (i.e.tracks reconstructed as cascades), the output may vary greatly, providing considerable separating power.
The track-reconstruction variables are: cosine zenith of the reconstructed direction, the distance of the track's first visible energy loss from the edge of the detector, the distance between the first and last visible energy loss, the number of direct hits (photoelectrons with DOM arrival time consistent with zero scattering), the track length between the first and last DOMs that registered a direct hit, the estimated angular uncertainty of the track reconstruction, and the angle between the track direction and a simplified reconstruction where first arrival times are fit to a propagating plane wave.
For the cascade reconstruction, the variables were the depth of the cascade within the detector, the horizontal distance between the cascade and the edge of the detector, the reduced log likelihood of the fit, and the log likelihood ratio of the track and cascade fits.There are also two variables that use multiple track fits to the same event.The first considers 104 different in-coming down-going track directions, ending at the reconstructed cascade, and counts the number of photoelectrons in the time window from -15 ns to +1000 ns around the geometric first arrival time from the tracks.The largest number of photoelectrons among all the tracks is selected for use in the BDT, and this helps to discriminate against downgoing muons that do not have a well-reconstructed track.The second considers 192 outgoing track hypotheses from all directions and counts the number of photoelectrons in the time window from -30 to +500 ns.This is important to identify starting tracks with a low-energy outgoing muon [51,52].
The BDT was trained with a sample of 4.5 million simulated atmospheric muon events from CORSIKA [32] that pass the outer-layer veto using SIBYLL2.1 [53] to model hadronic interactions.The spectrum was weighted to the H3a cosmic-ray flux model [54].The neutrino input was a total of 734,000 simulated neutrinos, weighted to the sum of the HKKMS conventional atmospheric flux [23], the BERSS prompt neutrino flux [28] and an astrophysical flux with a spectral index γ = 2.5, using the flux from Ref. [3].Several quality criteria were applied to the training sample: neutrinos labeled as starting tracks were required to have a vertex that was contained within the detector and produce a muon having a path length more than 300 m within the detector and energy above 100 GeV.Further, the reconstructed muon direction was required to be within 5 o of the simulated direction.
Figure 1 shows the number of events as a function of BDT atmospheric muon score and number of photoelectrons, for atmospheric muons, cascade-like and track-like neutrinos, and data.There are an estimated 1,800 times as many cosmic-ray muons as neutrinos in the sample, concentrated at low NPE and fairly high BDT muon scores.To optimally select neutrino events, an elliptical cut was used, where s µ is the BDT atmospheric muon score and (a, b) are parameters describing the semi-major and semiminor axes of the ellipse.Values of a = 0.75 and b = 2.5 were chosen to eliminate muon background, and the resulting ellipse is shown in Fig. 1.The total event rate for data, atmospheric muons, and neutrinos as a function of a in Eq. 6 while keeping b fixed is shown in Fig. 2 and shows good agreement between data and simulation.
The chosen value of a = 0.75 reduces the data event rate to 0.024 mHz or 3615 events in the final data sample.The CORSIKA Monte Carlo statistics are inadequate to estimate the number of surviving cosmic-ray muons in the final sample; an extrapolation based on the distance from the cut line in Fig. 2 indicates that the number of events should be negligible.To check this, an additional simulation was made using a parametrization of the flux of single muons in the deep ice calculated from the H3a cosmic-ray model [52].The simulation contains 26.8 million single muons passing the outer layer veto, but since it does not contain multi-muon bundles, it underpredicts the total muon background rate at this level by a factor of 3.However, most events near the BDT cut threshold are single muons, as can be verified by looking at COR-SIKA events, so the simulation can still be used to calculate a background estimate.The predicted background of single muons from the simulation is 2.7 ± 1.0 events in the final sample.Even conservatively assuming that an unexpectedly large contribution from bundles increases this by a factor of 3 as observed for the outer layer veto, the resulting muon background of 8.1 events is negligible compared to the full sample size of 3615 events and will no longer be considered.

B. Cascade/track classification
The same BDT was used to split the neutrino sample into cascades and tracks.The BDT track score, s track , and cascade score, s casc.are combined into a variable, ŝtrack = s track /(s track + s casc.), which runs from 0 to 1 and is an estimate of an event's "trackness," the likelihood that it contains a track.The final selection between tracks and cascades depends on the threshold chosen for ŝtrack .Figure 3 shows the purity, efficiency and signal to noise ratio (SNR) as a function of ŝtrack , as determined using simulated track and cascade samples.The SNR is defined as the ratio of the true track rate to the size of Poisson fluctuations in the total rate of true tracks and misidentified cascades.Optimizing SNR, the criterion ŝtrack ≥ 0.52 is used to identify tracks.Events not satisfying this criterion are identified as cascades.
With this classification, there were 965 cascades and 2650 track events in the final sample.The larger number of tracks is largely due to the dominance of ν µ over ν e in the atmospheric neutrino flux.

VI. DIRECTION AND ENERGY RECONSTRUCTIONS
IceCube has previously developed algorithms for reconstructing the direction and energy of cascades and track events.Events classified as cascades are reconstructed using a maximum likelihood fit as in Ref. [5], using the full photoelectron timing recorded by each DOM.Simulations predict that the median angle between the simulated and reconstructed direction is 16 • .However, starting tracks have received much less attention and call for new approaches.

A. Starting track reconstruction and inelasticity
Most IceCube analyses use track reconstructions that fit the track to a straight line by maximizing the likelihood for the reconstruction, based on functions that give the light amplitude and arrival time distribution func-tion for a given DOM, given a track hypothesis [55].For starting tracks, these reconstructions have two significant limitations.First, they assume that the track is infinite, originating outside the detector and traversing entirely through it.Second, they assume that the muon energy loss is continuous, rather than stochastic.The latter assumption does not hold for through-going tracks with energies above 1 TeV.However, it is much more problematic for starting tracks, which are accompanied by a large hadronic cascade from the recoiling nuclear target.Nevertheless, it is the best reconstruction that we have for finding the directions of starting tracks.
The directional resolution depends on both the neutrino energy and the inelasticity; the higher the inelasticity, the more the cascade dominates the event.That said, the median angular error is less than 2 • degrees for events with a visible inelasticity up to 0.9, rising to 5 • at a visible inelasticity of 0.99.Overall, simulations predict a 1.5 • median angular error for the entire starting track sample.These resolutions do not significantly impact the current analysis.
After the direction is found, the energy loss profile is unfolded as a sequence of electromagnetic cascades along the reconstructed track [41], as is shown in the middle panel of Fig. 4 for the most energetic starting track found in the sample.Generally, the largest cascade is at the interaction vertex.
This unfolding places some of the cascades outside the detector; these cascades have significantly larger uncertainties than those within the detector and are not used for the reconstruction.The energy loss profile is then integrated to give the cumulative fraction of the energy loss as a function of position along the track, as shown in the bottom panel of Fig. 4. The quantiles where the first 0%, 1%, 2%, etc. of total energy loss are determined, where 0% corresponds to the first observed loss.
Another machine learning technique, a random forest, is used to estimate the energy of the starting tracks.It takes the positions of the 101 energy loss quantiles as input, along with the total deposited energy, the total track length contained in the detector and the normalized track BDT score -a total of 104 inputs.The random forest is trained using simulated events and validated using an independent sample.For each event, the forest produces two outputs: the estimated visible cascade energy at the vertex, E casc., and the estimated track energy, E track .These are combined to produce two new variables: the total visible energy, and the visible inelasticity, Since the visible cascade energy is less than the hadronic energy, E vis. and y vis.tend to be lower than the actual neutrino energy and inelasticity, for CC ν µ interactions.Figure 5 shows the resolutions for these variables as quantified through the root mean square (RMS) error.The resolution for log 10 E vis. is better than either of its components, because there is some anti-correlation between the cascade and track energies.With this algorithm, for starting tracks, the RMS error on log 10 E vis. is 0.18, better than the typical resolution of 0.22 for through-going muons [56].The resolution for the visible inelasticity is 0.19.These performance metrics are mildly dependent on the assumed neutrino energy spectrum, which here is assumed to follow from the HKKMS conventional atmospheric flux, BERSS prompt atmospheric flux, and the best-fit power-law astrophysical flux from Ref. [3].

B. Cascade angular resolution check
Starting track events offer an opportunity to study cascade directional reconstruction, using the track as an indicator of true direction.This is possible because the muon and cascade are boosted to nearly the same direction, and the track angular resolution is much better than the cascade resolution.Track events with visible inelasticity y vis.> 0.75 are chosen to minimize the effect of the outgoing muon on the cascade reconstruction.This comparison is sensitive to both systematic offsets, as may be caused by improper modeling of optical scattering in the bulk ice, hole ice, DOMs, or other detector phenomena.Figure 6 shows the distribution of the difference in zenith angle between cascade and track reconstructions for a sample of high-inelasticity tracks.The data are reconstructed using the SPICE Mie ice model and shows a significant mismatch between the cascade and track zenith angle distributions, with the cascades be- ing reconstructed as more downward-going.This zenith angle distribution is sensitive to the amount of optical scattering in the bulk ice and in the hole ice, and we find increased scattering can produce a downward bias in the distribution.Our fits, discussed below, also confirm this observation and find somewhat higher levels of optical scattering than the IceCube baseline.

VII. INELASTICITY FIT
With the reconstruction results from Section VI, we are now ready to determine the inelasticity distributions of the starting track sample.Ideally, we would use these distributions to unfold dσ/dy, but there are several complications.The detection efficiency drops at low energies for very large y because the track is no longer visible.It also drops for small y at low energies because there is not enough light visible in the detector.Because of these strongly varying efficiencies, the limited statistics and the limited resolution, we do not attempt to unfold the data to present dσ/dy distributions.Instead, we parameterize the inelasticity distribution and fit to find the  [57].The result of fitting the distribution to the parameterization of Eq. 10 is shown with dashed green lines.The prediction of the CSMS differential CC cross section are shown for neutrinos with solid blue lines and antineutrinos with dashed blue lines.The total CC charm contribution is shown in magenta, illustrating its flatter inelasticity distribution.The best-fit flux models of Tab.II are assumed for all predictions.
parameters.These parameters can be compared with the Standard Model distribution, and also used to test alternative theories.
To motivate the parameterization, recall the differential CC cross section can be written schematically at leading order as where q and q represent sums of quark PDFs [16].High-energy neutrinos probe the PDFs at low values of Bjorken−x ∼ 3 × 10 −3 (E ν /PeV), where sea quarks dominate, and they should have a power-law behavior, xq(x, Q 2 ) ∼ A(Q 2 )x −λ with λ ∼ 0.4.Following Ref. [58], when transforming variables from (x, y) to (Q 2 , y), the Q 2 -dependence of Eq. 9 can be separated from the ydependence and integrated out to give a two-parameter function, where the parameter indicates relative importance of the term proportional to (1 − y) 2 in Eq. 10.This parameterization also works for antineutrinos, but takes on a different value since q and q are interchanged.Our measurement represents an average over neutrinos and  antineutrinos.The normalized inelasticity distribution can then be written as where N is the normalization This simple parameterization can accurately represent sophisticated calculations of inelasticity distributions.For example, a fit of Eq. 10 to the full CSMS calculation produces no more than a 1% root mean square deviation (averaged over y) for neutrino energies from 1 TeV to 10 PeV.
In practice, the parameters and λ are highly correlated when fitting Eq. 11 to realistic inelasticity distributions.To avoid this correlation, it is convenient to fit for the mean of the distribution, y , and λ instead, which show far less correlation.The mean inelasticity can be found analytically, It is then straightforward to substitute into Eq.11 so that dp/dy can be found as a function of y and λ only.
The visible inelasticity distribution is shown in Fig. 7 for four half-decade energy bins, from 1 TeV to 100 TeV; a 5th bin is used for energies above 100 TeV.The data are compared to predictions based on the CSMS cross section calculation, weighted by the expected neutrino and antineutrino fluxes.The flux models used are the bestfit atmospheric and astrophysical models to be described in Sec.VIII.The data are in good agreement with the predictions.
The visible inelasticity distribution in each energy range can then be fit to the parametrization of Eq. 11 using a binned Poisson likelihood fit of the 10 bins.The goodness-of-fit test statistic is where µ i (θ) is the expected event count in each bin depending on parameters θ and n i is the observed event count per bin [59].The second sum accounts for a Gaussian prior distribution on a parameter θ j with mean θ * j and standard deviation σ j .The expected event rate is derived from weighted Monte Carlo simulations.To account for the parametrized inelasticity distribution from Eq. 11, each simulated event receives a reweighting factor, dp dy (y; y , λ)/ dp dy CSMS , where dp dy CSMS is the inelasticity distribution calculated by CSMS that is used in the simulation.A total event rate scaling factor is also included to account for uncertainties in the flux normalization.The neutrino flux is assumed to follow the best-fit flux models in Sec.VIII, but the flux model and its uncertainties have negligible effect since the size of each energy range is comparable to the energy resolution.Detector systematic uncertainties on bulk ice scattering and absorption, DOM optical efficiency, and hole ice scattering are included through the use of 4 additional nuisance parameters in the fit.They are constrained by Gaussian priors and are further described in Sec.VIII.
The best-fit parameters for describing these inelasticity distributions are shown in Tab.I. Figure 8  CSMS calculation for neutrinos and antineutrinos.The measured values of y agree well with the flux-weighted average of neutrinos and antineutrinos.The downward trend in y is due to the W -boson propagator.

VIII. LIKELIHOOD FIT RESULTS AND STARTING TRACK/CASCADE COMPARISON
Inelasticity introduces a new dimension into the studies of high-energy astrophysical neutrinos, providing sensitivity to a number of new phenomena.In addition, the much more precise measurement of the starting tracks allows new tests, such as comparing the energy spectra of astrophysical ν µ with that from cascades, a mixture that is mostly ν e and ν τ .The inelasticity distribution is also sensitive to the ν/ν ratio and to neutrino interactions that produce charm quarks.In this section, we will present a baseline maximum likelihood fit and compare it with previous analyses.
The fit is done jointly over both cascades and starting tracks.For cascades, data is binned in two dimensions with half decade bins in energy ranging from 10 2.5 GeV to 10 7 GeV and 5 bins in cosine zenith angle.For tracks, data is binned in three dimensions with the same energy and zenith binning but additionally 5 bins in visible inelasticity.The same binned Poisson test statistic in Eq. 15 is used.
We first fit the data to a model that is similar to previous IceCube analyses [1][2][3]5].It includes three components: the conventional atmospheric neutrino flux, the prompt atmospheric flux and the astrophysical flux.After describing this standard fit here, Sections VIII.A-VIII.D discuss some additional fits which each add one additional degree of freedom to explore additional aspects of the physics.
The conventional atmospheric flux is based on the HKKMS calculation, extrapolated in energy to above 10 TeV and modified to include the knee of the cosmic-ray spectrum following the H3a cosmic-ray model.To account for the uncertainties in this flux model, we include several nuisance parameters in the fit.The first is the overall normalization, Φ conv .The second, ∆γ cr , accounts for uncertainty in the energy spectrum by allowing the spectral index to vary with a prior.A third parameter, R K/π , accounts for uncertainties in the kaon to pion ratio in cosmic-ray air showers [25].Neutrinos from kaons have a somewhat different zenith angle distribution than those from pions; R K/π accounts for that possible variation.The prompt atmospheric flux follows the BERSS calculation, an update of the ERS calculation [60] used in many previous IceCube works.It is incorporated into the analysis with a single parameter, the normalization for the overall amplitude.The self-veto probability is included for both atmospheric flux calculations.
Astrophysical neutrinos are initially assumed to be isotropic.In this section the ν e : ν µ : ν τ ratio is taken to be 1  3 : 1 3 : 1 3 ⊕ , an approximation expected from almost any conventional source model, after accounting for oscillation en-route to Earth.The ν : ν ratio is taken to be 1 : 1.The flux per flavor is assumed to follow a single power-law: where f α ≈ 1/3 is the fraction of each flavor at Earth, γ is the power-law index, and Φ 0 is a normalization factor that corresponds to the average flux of ν and ν per flavor at 100 TeV.
Detector systematic uncertainties are incorporated in all results through the use of four more nuisance parameters that describe uncertainties in the detection and transmission of light through ice.The first, DOM , accounts for uncertainties in the overall optical sensitivity of the DOMs; the prior on this is ±10%.Two parameters, α Scat.and α Abs., account for uncertainties of optical scattering and absorption in the bulk ice.These parameters linearly scale the inverse scattering and absorption lengths uniformly over all ice layers.Finally, α Hole Ice accounts for uncertainties on the overall scattering in the hole ice, the columns of refrozen drill water that the strings are emplaced in.Because of the presence of visible air bubbles [44] and possible impurities, the optical quality of this ice is expected to be much worse than that of the rest of the ice.The baseline ice model assumes a 50 cm scattering length in hole ice, but calibration data are also consistent with a scattering length of 30 cm, or 1.67 times more scattering.Uncertainties in both the hole ice and bulk ice scattering can lead to a bias in the zenith angle distribution of cascades as shown in Fig. 6.The fit finds values for these parameters that are in line with expectations.The larger value of α Hole Ice is comparable with other recent IceCube measurements [4].
The fit results are shown in Tab.II, and the fit is compared with the data in Fig. 9. Based on a set of computed pseudo-experiments, the goodness-of-fit test statistic, −2 ln Λ = 175.54,corresponds to a probability (p-value) of 0.04.This value indicates the fit model may not be a complete description of the data, perhaps due to inadequately simulated systematic uncertatinties.Still, the p-value is acceptable, and there are no obvious problem areas visible in the distributions in Fig. 9.We also fit the data with the optical properties of the ice fixed to their default values.The test statistic worsened signifi-cantly, but, except for R K/π , the other fit parameters did not change significantly.The inelasticity measurements are quite robust against systematic errors.

A. Astrophysical neutrino energy spectrum
The fit finds an astrophysical power law index of γ = 2.62 ± 0.07, with a normalization Φ 0 = 2.04 +0. 23 −0.21 × 10 −18 GeV −1 s −1 cm −2 sr −1 .This is in agreement with earlier IceCube studies of contained events [1], but softer than the most recent measurements of contained events [4], and considerably steeper than the measurement using  10.Confidence regions for the astrophysical power-law index, γ and flux normalization, Φ0.The blue contours show the confidence region for the joint fit of the cascade and starting track samples, which is the main result obtained here.The red contours show the confidence region for a fit of starting tracks only, and the green contours show the confidence region for a fit of cascades only, and all are consistent with each other.The confidence region from the IceCube analysis of through-going tracks [24] is shown in orange, which is in tension with the cascade confidence region.The contours from the starting track sample are consistent with both cascade and through-going samples.through-going tracks, which found γ = 2.13 ± 0.13 [24].There are a couple of possible explanations for the latter tension.First, the through-going tracks have considerably higher neutrino energies than the current sample, with 90% of the sample sensitivity in the energy range from 194 TeV to 7.8 PeV.Using the same method as in Ref. [24], the current sample has a 90% central range of 3.5 TeV to 2.6 PeV.For reference, 90% of the selected contained events are in the energy range from 3.3 TeV to 220 TeV; this range is much lower than for the sensitivity because the most energetic events have the largest effect on the astrophysical flux measurement.If the astrophysical flux is not a single power law, then one might measure different spectral indices in different energy ranges.Or, with difficulty, one might find different spectral indices for tracks and cascades.One way to test this hypothesis is to repeat the fit, allowing the astrophysical spectral indices and normalizations to vary between the starting tracks and cascades.When this is done, we find a power-law index of γ = 2.43 +0. 28 −0.30 for starting tracks and γ = 2.62 ± 0.08 for the cascades.The two indices are compatible within uncertainties. Figure 10 shows the two-dimensional confidence regions for the cascade and track measurements, the combined measurement, and the previous through-going track fit.Confidence regions are derived from the profile likelihood over all nuisance parameters, and it is assumed the test statistic follows a χ 2 distribution throughout.The cascade sample drives The ≈ 1 3 : 1 3 : 1 3 ⊕ composition at Earth, resulting from a 1 3 : 2 3 : 0 S source composition, is marked with a blue circle.The compositions at Earth resulting from source compositions of (0 : 1 : 0) S and (1 : 0 : 0) S are marked with a red triangle and a green square, respectively.The updated best-fit neutrino oscillation parameters from [61] are used here.Though the best-fit composition at Earth (black cross) is (0 : 0.21 : 0.79) ⊕ , the limits are consistent with all compositions possible under averaged oscillations.
the combined-sample index of 2.62 ± 0.07, by virtue of the much better energy resolution and lower atmospheric background compared to tracks.Within the uncertainty, the starting track power-law index is also compatible with that from the through-going tracks.We considered alternate scenarios with a double power-law or a power law with a cutoff that could explain a harder power-law index found at high energies, but we found no evidence for either when fitting our sample alone.All later results will continue the assumption of an unbroken power-law spectrum.
The other parameters in the fit in Tab.II are in line with expectations.The best-fit prompt flux is zero, in agreement with many previous IceCube studies [1,6,24,51], but the 1σ upper limit is compatible with the expected BERSS flux.The conventional flux, cosmic-ray spectral index and R K/π are all in line with expectations.

B. Astrophysical neutrino flavor composition
A related test of the astrophysical flux is to measure the flavor composition of the contained event sample.
Compared with the previous contained event analysis [5], this analysis benefits from much better track energy resolution and also the presence of the inelasticity distribution; the inelasticity distribution has some sensitivity to the presence of ν τ , since ν τ interactions, followed by τ → µν µ ν τ decays, will lead to events with larger visible inelasticity than ν µ interactions of the same energy.A global fit combining results from contained events and through-going muons found tighter limits that enabled constraints on the source flavor composition, however it compared cascades and tracks in different energy ranges where the energy spectrum may differ [6].
Figure 11 shows confidence levels for various ν e : ν µ : ν τ ratios.Each point in the large triangle corresponds to a specific flavor composition at Earth, with its associated confidence level obtained by fitting the data with the same parameters in Tab.II.The lines and points show the expectation from different production models and standard neutrino oscillations, including conventional pion decay with source flavor composition 1  3 : 2 3 : 0 S , neutron decay with (1 : 0 : 0) S , and a model where muons lose their energy via synchrotron radiation before they decay with (0 : 1 : 0) S [62].All of these conventional scenarios are along a narrow line, and, unfortunately, all are within the 68% confidence range for the analysis.The best-fit composition, (0 : 0.21 : 0.79) ⊕ , is on the left side of the triangle.However, because most ν τ produce cascades, there is a relatively high degeneracy between ν τ and ν e , so the confidence levels nearer the middle of the triangle are high; the break in the degeneracy from the inelasticity distribution of starting tracks is inadequate to statistically separate the ν τ and ν e components.In this analysis, astrophysical cascades and tracks have a similar energy range, with 68% central ranges of 11 TeV to 410 TeV and 8.6 TeV to 207 TeV, respectively.This is the tightest limit using samples with a similar energy range; a previous global fit [6] used contained cascades and through-going muons; the latter had a much higher energy range for through-going tracks (330 TeV to 1.4 PeV) than for cascades, where energies above 30 TeV were probed [63].With track and cascade samples having different energies, if the astrophysical spectrum is not a perfect power law, then the confidence regions on the flavor triangle will be shifted.
The extreme compositions of 100% ν e and 100% ν µ are ruled out with high confidence, at 5.8σ and 7.4σ, respectively.In fact, any composition with more than 2/3 ν µ is ruled out at 95% confidence level.These constraints can be used to put limits on exotic models.

C. Atmospheric neutrino/antineutrino ratio
At energies below about 10 TeV, neutrinos and antineutrinos have substantial differences in their inelasticity distributions, since neutrino interactions are more sensitive to quarks, while antineutrinos are more attuned to antiquarks.At large Bjorken−x values where valence 12.The predicted fraction of νµ contributing to the total νµ + νµ event rate in bins of reconstructed energy and inelasticity.At energies below ∼ 10 TeV, there are substantial differences in the inelasticity distribution that enable the atmospheric neutrino to antineutrino ratio to be measured.The bottom panel shows the fraction of atmospheric neutrinos contributing to the total event rate in bins of reconstructed energy.At energies above ∼ 100 TeV where astrophysical neutrinos begin to dominate the event rate, differences in the inelasticity distribution vanish, and it is not possible to measure the neutrino to antineutrino ratio for the astrophysical flux.An equal neutrino and antineutrino composition is assumed for the astrophysical flux here.
quarks dominate, the differences are substantial, leading to roughly a factor of 2 difference in cross section [64,65] as well as a difference in inelasticity distribution.Unfortunately, as Fig. 12 shows, the difference slowly disappears above 10 TeV.So, there is little sensitivity to the ν : ν ratio of astrophysical neutrinos.However, inelasticity can be used to measure the ν : ν ratio for atmospheric neutrinos.We do this by adding a scaling parameter to the list of parameters in Tab.II: R νµ/νµ scales the conventional atmospheric ν µ /ν µ ratio up or down.The best fit value is R νµ/νµ = 0.77 +0.44 −0.25 .A flux of 100% neutrinos (no ν) is disfavored at 3.8σ, while the reverse is excluded at 5.4σ.It should be noted that these limits are also dependent on the calculated angular distributions of ν and ν in addition to inelasticity.The sensitive range for this analysis is 770 GeV to 21 TeV; at higher energies, there is little ν : ν discrimination.This is the first ν : ν measurement in this energy range.Along with measurements of the atmospheric muon charge ratio [66][67][68][69], it is a useful diagnostic of particle production in cosmic-ray interactions [54].However, since this is an overall scaling factor, the same for all energies and zenith angles, it should not be directly interpreted in terms of hadronic interaction models since they also change the energy and zenith distribution assumed.

D. Neutrino charm production
Inelasticity measurements can also be used to probe charm production in neutrino interactions.The fraction of CC neutrino interactions that produce charm quarks rises slowly, from 10% at 100 GeV to 20% at 100 TeV.Charm quarks are produced primarily when a neutrino interacts with a strange sea quark; these sea quarks have lower mean Bjorken−x values than valence quarks, so the interactions have flatter inelasticity distributions than neutrino interactions with valence quarks.There is also a roughly 10% chance for a charm quark to decay semileptonically and produce an extra muon.These muons will not be distinguishable from the primary muon from the CC interaction; the energy loss from the two muons will add, and they will be reconstructed as a single, higherenergy muon, leading to a lower apparent inelasticity.For a ν e or ν τ interaction, the presence of a muon from semileptonic charm decay may lead to a track event with an apparent high inelasticity.
The contribution of charm production to different inelasticity events at different energies is shown in Fig. 13.Charm is most visible at energies above 100 TeV and at high inelasticity.For energies between 100 TeV and 1 PeV, more than 1/3 of the events with reconstructed y vis.> 0.8 produce charm.This shape difference can be used to search for charm production in a maximum likelhood fit.
An additional parameter, R CC,charm , scaling the CC charm production event rate is added to the parameters in Tab.II.Fitting cascades and tracks jointly, we find the 68% interval, R CC,charm = 0.93 +0. 73 −0.59 .The test statistic for a null hypothesis with zero charm production is −2∆ ln(L) = 2.8, so zero charm production is excluded at 91% confidence level.The 90% upper limit is 2.3 times the leading-order HERAPDF1.5 prediction.
The central 90% of neutrino energies contributing to this test statistic is from 1.5 to 340 TeV.This is a wider energy range than the central 90% of charm production events, 1.3 TeV to 44 TeV, because charm production is larger at higher energy.This is the highest energy measurement of charm production yet, and the upper end of the energy range is above the critical energies of charm hadrons where interactions in ice must occur.Similar methods could be used to search for other special types of neutrino interactions beyond charm production, including those beyond the Standard Model.

IX. CONCLUSIONS
We have developed a tool to measure neutrino inelasticity in Gigaton-scale H 2 O based detectors and presented the first measurements of neutrino inelasticity in very high-energy (above 1 TeV) neutrino interactions, using a sample of starting track events collected by the Ice-Cube Neutrino Observatory.The measured inelasticity distributions are in good agreement with the predictions   13.The predicted fractional contribution of all-flavor neutrino CC charm production events to the total event rate in bins of reconstructed visible energy and inelasticity.The increased charm production fraction at high visible inelasticity and high energy (up to 36%) provides a shape difference that allows the presence of charm production to be identified in a likelihood fit to the data.
of the CSMS calculation.
We have made a global fit to these neutrino data, fitting cascades in two dimensions: energy and zenith angle, and starting tracks in three dimensions: energy, zenith angle and inelasticity, to extract information about the astrophysical and atmospheric neutrino fluxes.This fit finds an astrophysical power-law spectral index of γ = 2.62 ± 0.07, in good agreement with previous fits to contained events and cascades, but in tension with previous results based on through-going muons, a sample that is generally higher in energy than the contained event samples.To explore this tension, we performed a fit where we allowed the astrophysical flux to float separately for cascades and starting tracks, with different spectral indices.Unfortunately, this leads to a spectral index for the tracks, γ = 2.43 +0.28 −0.30 , intermediate between the combined result and that for through-going tracks, with an error that is consistent with either.
We then relaxed the requirement that the astrophysical neutrino flavor ratio be 1  3 : 1 3 : 1 3 ⊕ , and calculated the confidence level for other compositions.We found a best fit point consisting of 79% ν τ , 21% ν µ , but with a broad allowed contour that encompasses all of the models that invoke conventional acceleration mechanisms and standard neutrino oscillations.More exotic models may be ruled out.
We also set limits on the ν µ : νµ ratio in atmospheric neutrinos and exclude zero production of charm quarks in neutrino interactions at 91% confidence level.This is only the second study, after measurements of the cross section using neutrinos with energies above 1 TeV.
Using the indirect signature in the inelasticity distribution, we observe, at greater than 90% confidence level, CC charm production in neutrino interactions, at energies between 1.5 and 340 TeV, more than an order of magnitude higher in energy than accelerator measurements.
Looking ahead, we expect that IceCube-Gen2 [70] and KM3NeT2.0 [71] will collect larger samples of contained events, which can be used to make more precise measurements of inelasticity.These detectors could collect substantial samples of events with energies above 100 TeV.With the increased precision, it will also be possible to study several new topics.Tau neutrinos are one example; ν τ interactions have a distinctive inelasticity distribution which could be used to detect a ν τ signal.Top quark production may also be accessible if enough energetic neutrinos are available.One calculation found that, for 10 PeV neutrinos, top quarks are produced in 5% of the interactions [17].It may also be possible to study other Standard Model neutrino interactions, such as diffractive production of W bosons in the Coulomb field of oxygen nuclei [21,22]; the cross section for ν +O → l+W + +X is about 8% of the charged-current cross section for 1 PeV neutrinos.Even the first phase of IceCube-Gen2 should enable improved calibrations of the existing data, reducing the systematic uncertainties.With moderately improved calibrations, the precision of the inelasticity measurements should scale as the square root of the effective volume times the live time.
With an improved surface veto to reject atmospheric neutrinos, it might also be possible to measure the ν : ν ratio of astrophysical neutrinos.If one could use the selfveto and a surface air-shower-array veto to reject atmospheric neutrinos with energies in the 1-10 TeV energy range, the inelasticity distribution could be used to determine the ν : ν ratio of astrophysical neutrinos.
The data could also be used to search for beyondstandard-model (BSM) physics, such as supersymmetry [75], leptoquarks [73] or quantum gravity with a relatively low scale [74].These phenomena also produce cross section enhancements which could be visible via increased neutrino absorption in the Earth, but the inelasticity distribution has a higher diagnostic utility than a simple increase in neutrino absorption.The use of inelasticity allows for a more sensitive search than by merely counting cascades and tracks [72].A combined fit to cross section and inelasticity measurements would provide even better constraints on new physics.
For most of these phenomena, the LHC provides better limits compared to IceCube Gen2 and KM3NeT 2.0.However, experiments that aim to record the coherent radio Cherenkov emission from ultra-energetic neutrinos with energies above 10 17 eV can reach supra-LHC energies.The ARA [76] and ARIANNA [77] collaborations both propose to deploy large (> 100 km 3 ) arrays that will, unless ultra-high energy cosmic-rays are primarily iron, collect useful (order 100 events) samples of cosmogenic neutrinos.The challenge here is that these experiments are primarily sensitive to cascades, while the energy deposition from tracks is too diffuse to be observ-able.However, it may be possible to take advantage of the Landau-Pomeranchuk-Migdal (LPM) effect to separate electromagnetic showers, which at energies above 10 20 eV are elongated, from the hadronic showers from the target nucleus, which are less subject to the LPM effect.This leads to a moderately elongated electromagnetic shower following a compact hadronic shower [78].With this, it might be possible to measure the inelasticity of charged-current ν e interactions [79].

6 EventsFIG. 1 .
FIG.1.Two-dimensional distributions of atmospheric muon BDT score and number of photoelectrons (NPE), after the outer-layer veto, for atmospheric muons (top left), cascadelike neutrinos (top right), track-like neutrinos (bottom left) and 10% of the data (bottom right).Events above the black lines were accepted as signal events.To access ∼ 1 TeV neutrinos that typically produce ∼ 100 PE, one must overcome muon background with a rate 1,800 times higher.
FIG.2.The total event rate for data (black), atmospheric muons (red), cascade-like neutrinos (green) and track-like neutrinos (orange) in the 10% testing sample as a function of the parameter a appearing in Eq. 6 describing the elliptical cut.The chosen value of a = 0.75 is shown with a vertical black line.

FIG. 3 .
FIG. 3. Top:The track purity and efficiency as a function of normalized BDT score ŝtrack .Bottom: The signal-to-noise ratio as a function of ŝtrack .A BDT track score cut ŝtrack ≥ 0.52 (black line) maximized the signal-to-noise ratio, yielding a 92% purity and 98% efficiency.The fluctuations in the purity curve are due to low statistics in that region.

4 FIG. 4 .
FIG.4.Top: An event display showing the most energetic starting track found in the data.DOMs are represented by colored spheres with a radius corresponding to the number of photoelectrons detected and with a color showing the first photoelectron time going from red (earliest) to green (latest).The larger blue spheres show the reconstructed sequence of electromagnetic cascades along the track, and their size is proportional to the reconstructed energy of each cascade on a logarithmic scale.The event originates in the top half of the detector.Middle: The reconstructed energy loss profile as a function of distance along the reconstructed track.The detector boundaries are shown by green and red dashed lines.Energy losses outside the detector, shown in grey, are excluded from the energy and inelasticity reconstruction.Bottom: The cumulative fraction of the total deposited energy within the detector.Percentiles of the energy loss distribution, shown with blue points, are features for the random forest regression of cascade and muon energy.The cascade and muon energies are estimated to be Ecasc.= 64 TeV and Eµ = 724 TeV, respectively, leading to Evis.= 788 TeV and yvis.= 0.08 for the visible energy and inelasticity, respectively.The total deposited energy is 135 TeV, and the muon escapes the detector with most of the neutrino's energy.

1 EventsFIG. 5 .
FIG. 5. Joint distributions of reconstructed quantity versus true quantity for cascade energy, Ecasc., muon energy, Eµ, total visible energy, Evis., and visible inelasticity, yvis..The root mean square (RMS) error for each quantity is shown at the top of each panel.
FIG.6.The distribution of the difference in zenith angle between cascade and track reconstructions for tracks with yvis.> 0.75.All reconstructions use the baseline SPICE Mie ice model.192 events meeting this condition were found in the full 5-year data sample.The data (black) show a mismatch between cascade and track zenith angles, with cascades being reconstructed with a median angle of 8.1 degrees more downgoing than the corresponding tracks.When events simulated using the SPICE Mie model are reconstructed (red), the median zenith angle difference is only 1.3 degrees.This discrepancy may be caused by systematic errors in the ice model since the cascade reconstruction is much more sensitive to the ice model than the track reconstruction.Distributions for two separate simulations increasing the bulk ice scattering globally by +10% (blue) and increasing hole ice scattering by +67% (green) are shown.Increased bulk ice scattering and hole ice scattering produce a shift that can explain the downgoing bias seen in the data.

5 5 FIG. 7 .
FIG.7.The reconstructed visible inelasticity distribution in five different bins of reconstructed energy.Observed data are shown in black.Error bars show 68% Feldman-Cousins confidence intervals for the event rate in each bin[57].The result of fitting the distribution to the parameterization of Eq. 10 is shown with dashed green lines.The prediction of the CSMS differential CC cross section are shown for neutrinos with solid blue lines and antineutrinos with dashed blue lines.The total CC charm contribution is shown in magenta, illustrating its flatter inelasticity distribution.The best-fit flux models of Tab.II are assumed for all predictions.

FIG. 8 .
FIG.8.The mean inelasticity obtained from the fit to Eq. 13 in five bins of reconstructed energy.Vertical error bars indicate the 68% confidence interval for the mean inelasticity, and horizontal error bars indicate the expected central 68% of neutrino energies in each bin.The predicted mean inelasticity from CSMS is shown in blue for neutrinos and in green for antineutrinos.A flux-averaged mean inelasticity per the HKKMS calculation is shown in red.
FIG.9.Best-fit distributions with all neutrino flux and detector parameters.Top: The distribution of reconstructed visible energy, visible inelasticity, and cosine zenith for the sample of starting tracks.Bottom: the distribution of reconstructed energy and cosine zenith for the sample of cascades.Contributions of conventional atmospheric and astrophysical neutrinos are shown in blue and red, respectively, and the total predicted distribution is shown in maroon.The prompt atmospheric neutrino contribution is not shown since its best-fit value is zero.The black error bars show the data.The inclusion of detector parameters describing bulk and hole ice scattering substantially improves the model description of the cascade cosine zenith distribution.The best-fit parameters are shown in Tab.II.
FIG.10.Confidence regions for the astrophysical power-law index, γ and flux normalization, Φ0.The blue contours show the confidence region for the joint fit of the cascade and starting track samples, which is the main result obtained here.The red contours show the confidence region for a fit of starting tracks only, and the green contours show the confidence region for a fit of cascades only, and all are consistent with each other.The confidence region from the IceCube analysis of through-going tracks[24] is shown in orange, which is in tension with the cascade confidence region.The contours from the starting track sample are consistent with both cascade and through-going samples.

FIG. 11 .
FIG. 11.Confidence regions for astrophysical flavor ratios (fe : fµ : fτ )⊕ at Earth.The labels for each flavor refer to the correspondingly tilted lines of the triangle.Averaged neutrino oscillations map the flavor ratio at sources to points within the extremely narrow blue triangle diagonally across the center.The ≈ 1 3 : 1 3 : 1 3 ⊕ composition at Earth, resulting from a FIG.12.The predicted fraction of νµ contributing to the total νµ + νµ event rate in bins of reconstructed energy and inelasticity.At energies below ∼ 10 TeV, there are substantial differences in the inelasticity distribution that enable the atmospheric neutrino to antineutrino ratio to be measured.The bottom panel shows the fraction of atmospheric neutrinos contributing to the total event rate in bins of reconstructed energy.At energies above ∼ 100 TeV where astrophysical neutrinos begin to dominate the event rate, differences in the inelasticity distribution vanish, and it is not possible to measure the neutrino to antineutrino ratio for the astrophysical flux.An equal neutrino and antineutrino composition is assumed for the astrophysical flux here.
FIG.13.The predicted fractional contribution of all-flavor neutrino CC charm production events to the total event rate in bins of reconstructed visible energy and inelasticity.The increased charm production fraction at high visible inelasticity and high energy (up to 36%) provides a shape difference that allows the presence of charm production to be identified in a likelihood fit to the data.

TABLE II .
compares y as a function of energy with the predictions of the Φ0/10 −18 GeV −1 s −1 cm −2 sr −1 2.04 +0.23The best-fit parameters including all neutrino flux and detector systematic uncertainties.68% confidence intervals are shown as calculated by the profile likelihood method.The prior distribution for each parameter is shown in the last column.The last row is the goodness-of-fit test statistic.