Probing Secret Interactions of Astrophysical Neutrinos in the High-Statistics Era

Do neutrinos have sizable self-interactions? They might. Laboratory constraints are weak, so strong effects are possible in astrophysical environments and the early universe. Observations with neutrino telescopes can provide an independent probe of neutrino self ("secret") interactions, as the sources are distant and the cosmic neutrino background intervenes. We define a roadmap for making decisive progress on testing secret neutrino interactions governed by a light mediator. This progress will be enabled by IceCube-Gen2 observations of high-energy astrophysical neutrinos. Critical to this is our comprehensive treatment of the theory, taking into account previously neglected or overly approximated effects, as well as including realistic detection physics. We show that IceCube-Gen2 can realize the full potential of neutrino astronomy for testing neutrino self-interactions, being sensitive to cosmologically relevant interaction models. To facilitate forthcoming studies, we release nuSIProp, a code that can also be used to study neutrino self-interactions from a variety of sources.

Allowed νSI would dramatically affect the evolution of systems with high neutrino densities, such as supernovae [18] and the early universe [4,[19][20][21][22]. While cosmological measurements have robustly established the presence of a radiation background compatible with expectations for the cosmic neutrino background (CνB) [23,24], its dynamics are poorly constrained. In fact, Ref. [25] found a preference for νSI in the Planck cosmic microwave background (CMB) data that could help alleviate the H 0 and σ 8 tensions [23,26] and impact constraints on inflation [27] (but see Refs. [28][29][30][31][32][33], where it is argued that polarization data restrict this possibility). Regardless of the outcome of that debate, it illustrates that a full understanding of possible νSI is critical to * Electronic address: esteban.6@osu.edu † Electronic address: phd1501151007@iiti.ac.in ‡ Electronic address: vedran.brdar@northwestern.edu § Electronic address: beacom.7@osu.edu our understanding of the early universe, especially in the upcoming high-precision cosmology era [34][35][36][37][38]. Figure 1 shows another essential point about νSI: while the constraints are strong for ν e , they are incomplete for ν µ and nearly nonexistent for ν τ [39] (this figure is explained in detail in Sec. II). This is why large cosmological effects from νSI remain allowed.
A new opportunity to probe νSI has arisen [40][41][42] due to the detection of high-energy neutrinos by Ice-Cube. The basic idea is that scattering of astrophysical neutrinos with the CνB en route to Earth redistributes their energies [43], leading to dips and bumps in the detected spectrum of the diffuse astrophysical neutrino background. Because we now know that all mass eigenstates have substantial ν τ components [44][45][46], this provides a new tool to explore ν τ self-interactions beyond the reach of laboratory experiments. However, with the relatively low statistics of IceCube data [47,48], it is hard to realize the full potential of this technique.
In this paper, we define a roadmap for making decisive progress. The first key to this is the proposed IceCube-Gen2 (Gen2 for brevity) detector [49]. As shown below, its principal advantage relative to IceCube, beyond its better statistics, is its much wider energy range. The second key is an improved theoretical treatment. We define the observable signatures of allowed flavor-dependent νSI, correcting errors and omissions in prior work. All calculations are done without any approximation, following careful optimizations. The third key is a realistic treatment of the detector response, using code provided by IceCube. To make future studies easier, we make our code publicly available at this URL . This paper is organized as follows. In Sec. II, we present the theoretical description of νSI, the relevance of , for α ∈ {e, µ, τ }). The hatched region is the "Moderately Interacting neutrino" (MIν) solution [25], argued to affect CMB observables. The dashed purple line is the interaction strength below which cosmic neutrinos free-stream as expected at cosmologically relevant times (see text for details). Above this line, our understanding of the early universe would be affected. As shown, ντ self-interactions are the least explored, leaving room for significant cosmological neutrino effects. different effects, and the interplay with other experimental observations. In Sec. III, we simulate the prospects of Gen2 to probe νSI. We show that Gen2 can probe the ν τ and other flavor sectors with sensitivity comparable to that of laboratory studies of the ν e sector. Thus, a combination of observables can fully probe the range of νSI that would significantly affect the dynamics of the cosmic neutrino background. In Sec. IV, we conclude and highlight future directions. In Appendix A, we give the results of the full νSI cross-section calculation.

II. SECRET NEUTRINO INTERACTIONS
In this section, we define the specific νSI models we consider. Sec. II A introduces the theoretical description and flavor-dependent bounds on νSI. Sec. II B specifies the spectral signatures that we expect given our present understanding of neutrino masses and mixings. Finally, in Sec. II C, we demonstrate, by computing the full νSI cross-section, the relevance of previously ignored effects both in the theoretical and observed high-energy astrophysical neutrino spectrum.

A. νSI Model and Conceptual Framework
We consider νSI parametrized by the Lagrangian where ν α are neutrino flavor eigenstates (α ∈ {e, µ, τ }), g αβ is the interaction strength between flavors α and β, φ is the interaction mediator (that for simplicity we assume to be real), M φ is its mass, and the (1/2) prefactor is present if neutrinos are Majorana fermions. We write the coupling strength in the flavor basis, as in this basis laboratory constraints are simpler to express. This is related to the coupling in the mass basis by g ij = α, β where U is the leptonic mixing matrix and i ∈ {1, 2, 3}. We assume scalar interactions, as pseudoscalar interactions give the same results at high energies and vector interactions are strongly constrained from laboratory experiments [50,51] and Big Bang Nucleosynthesis [39,52]. Although the Lagrangian in Eq. (1) is not gauge invariant due to the new couplings affecting neutrinos but not charged leptons, in this paper we are only interested in the low-energy phenomenology effectively described by Eq. (1). For possible UV completions that do not generate large couplings to charged leptons, see Refs. [39,[53][54][55]. These UV completions might induce additional signatures such as relatively large mixings with sterile neutrinos or a modified Higgs phenomenology, but these are beyond the scope of this paper. We assume Majorana neutrinos, as for Dirac neutrinos stronger BBN constraints would apply [39] (see Sec. II B for comments on the phenomenology of Dirac neutrinos). Figure 1 shows that νSI, particularly in the ν τ sector, are poorly constrained and can affect our understanding of the early universe. The present constraints (at 2σ) are shown in shaded contours, which come from Z-boson decay [56], charged kaon decays [39,53,57], and Big Bang Nucleosynthesis (BBN; obtained with the AlterBBN 2.2 code [58,59]) . To show the scale where early-universe  effects become important, the dashed purple line indicates the interaction strength below which cosmological  neutrinos free-stream as expected at times relevant to  the evolution of the CMB multipoles observed by Planck  (  2500). We conservatively assume this to be the case if the neutrino self-interaction rate (decreasing as the universe expands) is already below 10% of the Hubble expansion rate H when the smallest Fourier modes observed by Planck are still well outside the horizon (i.e., for a scale factor a = 10 k/H, with k the comoving wavenumber corresponding to ∼ 2500).
As an example of νSI that is presently allowed, the hatched region indicates the Moderately Interacting Neutrino (MIν) solution [25,39], which has been argued to affect cosmological parameter extraction from CMB data, especially the Hubble constant H 0 and the amplitude parameter σ 8 . Even if the MIν solution fades away once more data is accumulated (there seem to be indications in that direction from Planck polarization data [28][29][30][31]33]), it will remain important to probe the full parameter space above the purple dashed line.
The weakest νSI constraints are in the ν τ sector. Because of flavor mixing, astrophysical neutrinos must always contain a large ν τ component. In what follows, we explore how νSI affect astrophysical neutrino propagation, assuming that νSI apply only in the ν τ sector, i.e., g αβ = g δ ατ δ βτ (our code nuSIProp, though, can handle new couplings in all sectors). In Sec. III, we comment how our results in fact provide comparable sensitivity to all three neutrino flavors.
The basic physics of νSI affecting astrophysical neutrino propagation is as follows [40][41][42][43]. En route to Earth, high-energy astrophysical neutrinos may scatter with neutrinos in the CνB. (In the standard model, scattering is irrelevant except perhaps at ultra-high energies through the Z resonance [60,61].). As a consequence, high-energy neutrinos are absorbed and lowerenergy neutrinos are regenerated. This leads to unique dips and bumps in the astrophysical neutrino spectrum.
To quantify this description, we follow the comoving differential neutrino number density per unit energy E ν , dn/dE ν , as a function of cosmological time t. Since neutrino oscillation lengths are much smaller than astrophysical scales, flavor oscillations are quickly averaged out; mass eigenstates rapidly decohere too [62]. We can thus follow the comoving differential density of neutrinos of mass eigenstate i ∈ {1, 2, 3}, dn i /dE ν ≡ñ i (t, E ν ). Its evolution is dictated by a transport equation [41,63] The first term on the right-hand side accounts for the energy redshift due to cosmological expansion; H is the Hubble parameter. The second term represents the production rate of high-energy neutrinos from astrophysical sources: we assume it follows the star formation rate [64,65] as a function of time (a different redshift dependence has little impact on the results as it can not induce dips and bumps in the neutrino spectrum [66]), and a power law ∝ E −γ ν with spectral index γ as a function of energy. The first two terms would thus describe the evolution of astrophysical neutrinos en route to Earth without new physics. The third term accounts for absorption due to νSI, while the fourth term accounts for the subsequent regeneration. In this equation, n t i 2 × 56 (1 + z) 3 cm −3 = 8.7 (1 + z) 3 × 10 −13 eV 3 is the CνB density of ν i , with z the cosmological redshift, σ ij (E ν ) is the absorption cross-section of an incident ν i with energy E ν on a target ν j , and σ jk→il (E ν , E ν ) is the cross-section for an incident ν j with energy E ν on a target ν k to generate a detectable ν i with energy E ν , along with a ν l . For Majorana neutrinos, all final states are detectable; for Dirac neutrinos, neutrinos must be lefthanded and antineutrinos right-handed.
We ignore the neutrino matter potential induced by νSI [67][68][69]. A priori, it could affect the flavor composition at Earth, but for scalar νSI it just generates a negligible correction to neutrino masses [70]. For vector νSI, however, this effect might be measurable at higher energies than the ones we consider if the cosmic neutrino background has a large lepton asymmetry [67]. By solving the differential equation (2) with the initial conditionñ i (z → ∞, E ν ) = 0, we evaluatẽ n i (z = 0, E ν ) ≡ φ i (E ν ), the neutrino flux at Earth with energy E ν . The details on how we numerically solve Eq. (2) are given at this URL .

B. νSI Effects on the Spectrum
Qualitatively, the effects of νSI can be understood through the s-channel neutrino-neutrino scattering crosssection (which is dominant, but not to the exclusion of other terms; see Sec. II C), where s ≡ 2E ν m j is the center-of-momentum energy squared, m j is the mass of ν j , and Γ = g 2 M φ /16π is the total scalar decay width. (Note that there is a factor of √ 2 difference in our definition of g with respect to Ref. [41].) At s = M 2 φ , the cross-section is resonantly enhanced, , leading to an enhanced astrophysical neutrino absorption. Since U τ i ∼ O(1) for all mass eigenstates, the astrophysical spectrum will generically feature multiple absorption dips, located at E ν ∼ M 2 φ /(2m j ) with j ∈ {1, 2, 3} [53,66,71]. 2.5σ over the other two scenarios, would produce two close dips in the astrophysical neutrino spectrum.
The separation between these dips, and thus their separate observability, depends on the neutrino mass spectrum. In the last few years, there has been a lot of progress in understanding it: cosmological data bounds the total neutrino mass m ν < 0.12 eV at 95% CL [23], and neutrino oscillation data prefer the Normal mass Ordering (NO), i.e., m 3 > m 2 > m 1 , over the Inverted mass Ordering (IO), i.e., m 2 > m 1 > m 3 , with ∼ 2.5σ [44]. (As is standard in neutrino oscillation studies, we denote as ν 1 and ν 2 the eigenstates that are closest in mass, with m 2 > m 1 . The other eigenstate is denoted as ν 3 .) Admittedly, the cosmological bound is subject to some model dependence [72][73][74][75], and the oscillation preference has lately weakened [44]. However, both hints are already statistically significant and the situation should improve quickly. For simplicity, we calculate our results assuming the NO. As we show below, the main difference between the NO and the IO appears in present IceCube data, which is not our focus. For Gen2, we show that it has comparable sensitivity to νSI for both mass orderings.
In the NO, νSI in the ν τ sector should feature two absorption dips, separated in energy by an O(1) factor.
The reason for this is as follows. First, m ν < 0.12 eV, together with the measured neutrino squared mass splittings |∆m 2 32 | ∼ |∆m 2 31 | ∼ 0.05 eV, implies that the mass eigenstates are not degenerate. At the same time, for the NO the spectrum contains one heavy state and two light states (as |∆m 2 31 | ∆m 2 21 , i.e., m 3 m 1 ∼ m 2 ); whereas for the IO there are two heavy states and one light state (m 2 ∼ m 1 m 3 ). Thus, for the same total neutrino mass, the lightest mass for the NO (m 1 ) will typically be much larger than the lightest mass for the IO (m 3 ). As a consequence, the NO will generically imply that light and heavy mass eigenstates are closer to each other, i.e., for the NO the νSI absorption dips will be closer to each other. Although in principle there are three absorption dips (one for each mass eigenstate) the smallness of ∆m 2 21 implies that two of them are essentially at the same energy unless 0.06 eV m ν 0.07 eV [76]. Figure 2 shows these features. We display the allflavor high-energy astrophysical neutrino flux at Earth obtained by numerically solving Eq. (2). We plot this as E 2 ν dφ/dE ν 2.3 −1 E ν dφ/d log 10 E ν , so that the conserved area under each curve represents the conservation of energy (for the Majorana case, as discussed below). For clarity we have assumed an astrophysical neutrino spectrum ∝ E −2 ν (steeper spectra would make the bumps less visible); below we relax this assumption. We observe the two dips due to neutrino absorption, as well as the two bumps due to the scattering products moving to lower energies. Although absorption mostly takes place at E ν = M 2 φ /(2m j ), neutrino energy gets redshifted as the universe expands, extending the dips to lower energies. While the IO also predicts two absorption dips for m ν = 0.1 eV, they are separated by about two orders of magnitude in neutrino energy. As the present IceCube data covers a rather narrow energy range, this would appear as a single dip in the spectrum. Future neutrino telescopes as Gen2 (see Sec. III) should have a wider energy range and thus νSI could feature two dips even for the IO.
The fact that we are considering interactions only in the ν τ sector could have an imprint on the flavor composition of the observed high-energy astrophysical neutrino flux. However, we find the effects to be modest, as shown next. The Dirac or Majorana nature of neutrinos may also affect the results. If neutrinos are Dirac fermions, part of the final scattering states will be unobservable (c.f. Table I in Appendix A). As a consequence, the characteristic νSI bumps will be smaller. This distinction between Dirac and Majorana neutrinos is expected in any new physics scenario where neutrinos have non-chiral interactions [77]. Figure 3 (top) shows that measuring the all-flavor flux is enough to explore the flavor-dependent νSI we consider. We show the high-energy astrophysical neutrino flux at Earth for each neutrino flavor, obtained by numerically solving Eq. (2) for m ν = 0.1 eV and the NO. Assuming ν τ self-interactions implies that the scattering products will be τ neutrinos, but since U flavors, these will quickly oscillate to an almost flavoruniversal composition at Earth. A very good experimental flavor-discrimination, though, might be advantageous in the future as a signature of flavor-dependent νSI [78,79].
Figure 3 (bottom) shows the difference between Dirac and Majorana neutrinos. We assume the same physics parameters as in Fig. 3 (top), but by requiring neutrinos to be Dirac fermions, we observe smaller appearance bumps. This effect can be understood from Eqs. (1) and (3): if scattering is mediated through the s-channel resonance, it can be understood as on-shell production of φ followed by its decay. From Eq. (1), we see that the decay products must have opposite chiralities, and so for Dirac neutrinos one of them will be sterile. Since this effect is subleading and, furthermore, constraints on Dirac neutrino self-interactions are rather strong [39], we assume Majorana neutrinos for the rest of this paper.
In short, our precise understanding of neutrino properties clearly identifies the experimental signatures we must look for to make the most of the opportunity suggested by Fig. 1 and explore tau neutrino self-interactions: multiple dips in the astrophysical spectrum with an almost flavor-independent flavor composition. Top: All diagrams that contribute to neutrino-neutrino scattering at lowest order. For doublescalar production, the scalars quickly decay to neutrinos. Bottom: Total neutrino self-interaction cross-section divided by the s-channel contribution as a function of the scaled center-of-momentum energy squared. At high energies, the s-channel contribution is subdominant.

C. Importance of the Cross-Section Calculation
In Eq. (2), the neutrino self-interaction cross-section σ parametrizes all the necessary particle physics inputs. Here we show that approximations made in the literature were unjustified in some important cases. So far in the literature, only the s-channel contribution (Eq. (3)) has been taken into account [40,41,63,80,81]. Furthermore, as Eq. (2) can be computationally expensive to solve, Eq. (3) was recently approximated [17,63,81] with a δ-function at s = M 2 φ . Figure 4 shows that other contributions to the cross section (see Appendix A) are in fact dominant for s > M 2 φ . One can check that for coupling strengths g ∼ O(0.1) -according to Fig. 1, the couplings that impact our understanding of the early universe -the astrophysical neutrino interaction probability is large even for s > M 2 φ , and so other scattering channels must be taken into account. Figure 5 shows that these other scattering channels are quantitatively important. We show the all-flavor high-energy astrophysical neutrino flux at Earth, obtained by numerically solving Eq. (2) under different approximations for different lines; in all cases we assume m ν = 0.1 eV and the NO. As we see, for coupling strengths O(0.1), all contributions to the cross-section must be included to correctly predict the neutrino spectrum, and the δ-function approximation is unreliable. For g ∼ 0.01, on the other hand, all approximations are quantitatively justified. Note that, for g = 0.1, the error introduced by only including s-channel scattering is smaller than the naive expectation from Fig. 4. This is because, for non s-channel scattering, the final-state neutrinos tend to have energies relatively close to the initial neutrino, and so regeneration partially compensates absorption.
We also show that the δ-function approximation, i.e., only taking into account the scattering cross-section exactly at the resonance, gives the same result for both couplings (g = 0.01 and g = 0.1) in Fig. 5. This can be understood from the neutrino mean free path at the resonance, which, ignoring the expansion of the universe, is [63] λ res n t σ res which is much smaller than the typical distances to astrophysical sources, ∼ 1/H 0 ∼ Gpc, as long as g 10 −3 (note that oscillation data imply |U τ i | = 1 [44][45][46]). That is -assuming mediator masses O(1-10 MeV) -for all couplings g 10 −3 all astrophysical neutrinos will be scattered (at E = E res ) regardless of g, and so the δfunction approximation will always give the same result. This highlights the importance of carrying out the full theoretical calculation to obtain reliable results. Figure 6 (top) shows that the effects discussed above are also present after including detection effects. We show in the top panel the expected number of events at IceCube in each bin of deposited energy in the detector, E dep . (Note that here we plot the number of events per log bin, not the energy-weighted version of that, so now the spectrum is falling.) It is important to use E dep instead of E ν as the former is different for each flavor, and takes into account the detector energy resolution [82]. Using E ν would overestimate the detectability of energydependent effects.
To simulate the detector response, we use the official high-energy starting event (HESE) data release [48], which includes all neutrino flavors and interaction topologies and takes into account Earth [80] and detector effects. We adjust the flux normalization and exposure to predict 10 events in the lowest energy bin as approximately observed by IceCube after 7.5 years (see Sec. III A), and for νSI we assume m ν = 0.1 eV and the NO. We can see the two distinct dips, appearing at E dep ∼ 2.5 × 10 5 GeV and E dep ∼ 5.5 × 10 5 GeV, in the presence of νSI. The effect of the appearance bump is that, for E dep < 2 × 10 5 GeV, the neutrino spectrum is steeper than without νSI. From here, we already foresee a degeneracy between νSI and the spectral index. Breaking it will require observing the spectrum at a wide neutrino energy range.
The bump at E dep ∼ 6 × 10 6 GeV is the Glashow resonance [83], i.e., resonant production of a W boson in neutrino-electron scattering. This feature has been recently observed by IceCube [84] and provides a window to the high-energy spectrum. As shown in Sec. III, this can help in exploring νSI. Figure 6 (bottom) shows that, for couplings g ∼ 0.1, not including the full cross-section largely overestimates the number of events. We show there the same highenergy astrophysical neutrino flux in the presence of νSI as in the top panel, but for different levels of approximation when solving Eq. (2). Such approximations were mostly used in prior work to speed up the computations. In our numerical code they are largely unnecessary, as it has been optimized to have a low computational cost. As we make the code publicly available, it can be easily employed in further studies.

III. ICECUBE-GEN2: THE ROAD TO PRECISION NEUTRINO ASTROPHYSICS
In this section, we show that Gen2 will be a unprecedented tool for probing νSI. In Sec. III A, we review the present IceCube data and discuss its limitations for exploring νSI. In Sec. III B, we discuss how these limitations will be overcome by Gen2. Finally, in Sec. III C, we carry out a data analysis to quantify the future sensitivity of Gen2 to νSI as well as the present exclusion by IceCube.

A. Status of IceCube Data
The IceCube observatory has firmly established the existence of high-energy astrophysical neutrinos by detecting O(100) neutrinos with O(0.1 − 1) PeV energies. As νSI introduce an energy-dependent distortion in the as-trophysical neutrino spectrum, probing them requires a good reconstruction of the neutrino energy. The HESE dataset [85] is particularly appropriate for this, as it only includes neutrinos interacting within the internal part of the detector, thus reducing the amount of energy deposited outside the instrumented volume. Figure 7 shows that, intriguingly, the IceCube HESE data suggests some non-trivial spectral features. Highenergy astrophysical acceleration mechanisms generically predict that the neutrino spectrum should follow a power law as a function of energy [86,87]. This would roughly correspond to a straight line in Fig. 7, but the data shows some deviations. In particular, there seems to be a deficit of events at E dep ∼ 5 × 10 5 GeV, and an excess at E dep ∼ 10 6 GeV. These features have been interpreted in the literature as hints for νSI [17,80,81].
To illustrate this point, we show in Fig. 7 that the spectrum can be fit by nonzero νSI and a different spectral index. While the fit is visually better than the no-νSI case, we caution that the significance is not high (see Sec. III C). This νSI interpretation relies on the presence of two separate dips (at E dep ∼ 5 × 10 5 GeV and at E dep ∼ 3 × 10 6 GeV), each with its corresponding bump at lower energies. We also note that the degeneracy between νSI parameters and the spectral index could perhaps be exploited to reduce the longstanding tension in the different values of the spectral index deduced from different IceCube datasets [85,88,89]. However, this is beyond our scope.
This discussion identifies two weaknesses of the present IceCube data to explore νSI: apart from having relatively low statistics, the data only covers about one order of magnitude in neutrino energy. As Figs. 2 to 6 show, νSIinduced spectral features cover a wide energy range. Furthermore, a wider energy coverage also entails a better understanding of the spectral index that to some extent can be degenerate with νSI-induced dips, as Fig. 7 shows.

B. Prospects for Gen2
Fortunately, Gen2 [49] will overcome the issues noted above for IceCube. Its effective volume will be about one order of magnitude larger, so Fig. 7 already indicates that Gen2 should observe events at energies up to at least E dep ∼ 10 7 GeV. (Other neutrino telescopes [90][91][92] could have comparable sensitivity as long as their effective volume is as large as that of Gen2 and their energy resolution is good enough.) To quantify the potential of this future observatory, we carry out a simplified but realistic simulation of Gen2. We take the Gen2 optical effective area from Fig. 25 in Ref. [49] that takes into account Earth attenuation, which is significant at Gen2 energies; for simplicity we neglect ν τ regeneration that is less important, specially for steep power laws. We compute the deposited energy distribution following Refs. [82,93], and we assume an energy resolution of 10% for the deposited energy. Following Ref. [49], we assume a minimum detectable neutrino energy of 200 TeV, below which the data may contain too many atmospheric background events. Figure 8 shows that Gen2 will indeed be very sensitive to νSI. We show the expected numbers of events after 10 years of data taking as a function of deposited energy, using the same physics parameters as in Fig. 7. The bottom panel shows the difference between the expected number of events with and without νSI, divided by the square root of the former, which gives a conservative underestimate of the expected statistical uncertainties. As we see, unlike presently at IceCube, the much wider energy range, together with the increased statistics, will undoubtedly disentangle a power-law spectrum from νSI with a modified spectral index. To better understand the sensitivity of Gen2, though, we must explore the impact of weaker νSI. Figure 9 quantitatively shows that even relatively weak νSI would profoundly impact the Gen2 observations. We show the expected spectrum in Gen2, assuming a normalization compatible with current IceCube observations. In the bottom panel, we show the difference between the number of events with and without νSI, divided by the expected statistical uncertainties.
The dotted purple line is of particular physical relevance. It corresponds to the coupling at which, for M φ = 10 MeV, the neutrino mean free path (c.f. Eq. (4)) is 1 Gpc, roughly the distance to astrophysical neutrino sources. We expect the effects shown by that line to be close to the weakest detectable νSI propagation effects: for smaller couplings, the neutrino flux attenuation where L is the distance to sources, is hardly different from 1, and uncertainties on the underlying astrophysical flux could start playing an important role. According to the bottom panel, Gen2 should have enough statistical power to be sensitive to the dotted purple line. It will thus realize the full potential of high-energy astrophysical neutrino propagation for testing νSI.
The results above also illustrate that the unique spectral features cannot be accommodated by a modified power law and so would be a smoking gun for νSI. Note that it is not easy for random statistical fluctuations to emulate νSI, because (1) the separation between the absorption dips is precisely determined by external neutrino oscillation and cosmological data, and (2) νSI predict a bump with a fixed amplitude just below each energy dip. These features suggest that, should spectral anomalies appear at Gen2, it could be relatively easy to identify them as being due to νSI or not.

C. νSI sensitivity
To quantify the future sensitivity of Gen2, along with the νSI presently allowed by IceCube, we carry out a data analysis using the theoretical framework discussed in Sec. II. Starting with the present IceCube HESE data, we follow the procedure in the official data release [85,94]. We assume an initial neutrino flux following a power law with spectral index γ ∈ [2, 3] and a normalization φ 6ν (100 TeV) ∈ [10 −16 , 10 −20 ] GeV −1 cm −2 s −1 sr −1 . In addition to the standard analysis parameters, we include ν τ self-interactions parametrized by the coupling strength g, the mediator mass M φ , and the total neutrino mass m ν . For m ν , we include as a prior the cosmological bound coming from CMB and Baryon Acoustic Oscillation data [23], m ν < 0.12 eV (95% CL). We assume the NO and the neutrino squared mass splittings and mixings from Ref. [44].
To forecast the expected sensitivity of Gen2, we generate mock data corresponding to 10 years of exposure at Gen2. As IceCube will keep accumulating data, we also include 15 years of IceCube exposure for deposited energies E dep ∈ [6 × 10 4 , 2 × 10 5 ] GeV that are below the range of Gen2 [49]. This increases the sensitivity to lower mediator masses. We assume an astrophysical flux following a power law, φ 6ν (E ν ) = 5 × 10 −18 (E ν /100 TeV) −2.5 GeV −1 cm −2 s −1 sr −1 , compatible with present IceCube data. We analyze the mock data assuming νSI as outlined in the paragraph above (including marginalizing over the spectral index).  Fig. 1). The dark green region, including part of the MIν region, is excluded by IceCube data, and the dashed brown line shows the Gen2 sensitivity (2σ). Gen2 will exploit the full potential of testing νSI with high-energy astrophysical neutrino propagation (see text), being sensitive to a large parameter space where neutrinos have a non-trivial cosmological behavior. The sensitivity to other flavors is comparable.
We use only the HESE spectrum, and do not exploit flavor information. As Ref. [79] shows, Gen2 will have very good flavor discrimination. Including it could enhance the sensitivity to νSI and probe its flavor structure. Figure 10 shows that Gen2 will have superb sensitivity, covering a huge range of cosmologically relevant νSI parameters. The shape of the Gen2 sensitivity curve can be easily understood. Its best sensitivity to the coupling strength g covers about one order of magnitude in M φ which, for resonant scattering (i.e., M 2 φ ∼ 2E ν m ν ) corresponds to two orders of magnitude in E ν . As we see from Fig. 8, this is roughly the energy range of Gen2.
The abrupt sensitivity decrease for M φ < 2 MeV is because, at these mediator masses and for m ν = 0.1 eV, the highest-energy νSI resonant absorption dip is below the smallest detectable neutrino energy. Gen2 will only be sensitive to such small mediator masses if neutrinoneutrino scattering outside the resonance is relevant. This, as Fig. 5 shows, only happens for coupling strengths g ∼ 0.1. These results also illustrate that, since the lowenergy statistics is very large (c.f. Fig. 9), the presence of a single spectral feature is enough to probe νSI.
For mediator masses between 2 MeV and 20 MeV, Gen2 has superb sensitivity. The sensitivity there can be approximated as g/M φ ∼ several times 10 −4 MeV −1 . This corresponds to a neutrino mean free path at the resonance (c.f. Eq. (4)) of λ res ∼ Gpc ∼ H −1 0 , the typical distance to astrophysical neutrino sources. For lower couplings, the neutrino flux attenuation ∼ e −g 2 /M 2 φ n t L ∼ 1 − g 2 /M 2 φ n t L, where L is the distance to sources, is hardly different from 1. Because Gen2 will have good statistics over a wide energy range, it will be hard to improve upon this sensitivity for mediator masses between 2 MeV and 20 MeV. For example, increasing the statistics by a factor N would only improve the sensitivity by g ∼ N 1/2 , even in the ideal case of only having statistical uncertainties.
There are also some local sensitivity improvements for mediator masses of ∼ 4 MeV, ∼ 20 MeV, and ∼ 30 MeV. These correspond to the mediator masses at which νSI spectral features either enter the Gen2 spectrum; or have the same energy as the Glashow resonance, which increases their statistics.
The sensitivity decreases again for M φ 20 MeV, the mediator masses above which the highest-energy νSI dip is at energies E dep 10 7 GeV. Although there are still other νSI features in the spectrum, the poor statistics at high energies dramatically reduces the sensitivity. For these mediator masses, in contrast to 2 MeV M φ 20 MeV, additional statistics from higher energy detectors could significantly increase the sensitivity to νSI. Note that this particular mediator mass above which Gen2 loses sensitivity depends on the spectral index, which sets the highest energy neutrinos that can be detected.
We have systematically explored how our results depend on the input choices. In short, changing the inputs within reasonable ranges leads to results that are similar enough; hence, Fig. 10 is sufficient to represent the Gen2 sensitivity. Assuming the IO instead of the NO would keep one of the spectral features unaffected (c.f. Fig. 2), which, as discussed above, is enough to probe νSI. The assumed total neutrino mass or squared mass splittings do not change the results either, as present oscillation and cosmological data are precise enough to fix the location of the relevant νSI spectral features [76]. Finally, the assumed spectral index sets the largest energy detectable by Gen2 and thus the largest value of M φ that can be explored, but we have checked that changing it to 2.0 or 3.0 only changes the higher end of the Gen2 sensitivity by about a factor of 2 in M φ . As we discuss below, the higher end of the Gen2 sensitivity can be complemented with higher-energy observatories. Fig. 10 also shows the present IceCube HESE exclusion region. It is more powerful than BBN and Z decay bounds in some regions of the νSI parameter space, and complementary in others. The explored coupling range already indicates the relevance of the theoretical effects we discuss in Sec. II: we have checked that, if we ignored them, the IceCube excluded region would be quantitatively different.
Finally, we note that, unlike in prior work [17,80,81], we do not find any statistically significant indication in favor of νSI in present IceCube data. In particular, despite being described by more parameters, νSI are only preferred over a simple power law neutrino spectrum at the ∼ 1σ level. This difference may stem from our more precise theoretical treatment described in Sec. II, our use of a more recent data set, our precise treatment of detector effects that matches that used by IceCube, or a combination of these effects. Both Figs. 9 and 10 illustrate that Gen2 will be a very powerful tool to explore νSI. We expect it to be sensitive to the full range of the presently allowed MIν solution and to almost the full parameter space for which cosmological neutrinos do not free-stream. The remaining allowed region could be explored with higher-energy observatories including Gen2-radio [49,[95][96][97][98][99][100][101]. Furthermore, for mediator masses between 2 and 20 MeV, Gen2 would realize the full potential of neutrino astronomy to test νSI, even competing with the strongest laboratory probes (see meson-decay limits in Fig. 1) and setting the strongest bound on νSI. Since the presence of a single spectral feature is enough to probe νSI at Gen2, we expect the constraints from Fig. 10 to apply with comparable strength to ν e and ν µ . This future observatory will be a unique bridge between the cosmological and laboratory quests to understand if neutrinos have large self-interactions.

IV. CONCLUSIONS AND OUTLOOK
We return to our opening question: Do neutrinos have sizable self-interactions? In the laboratory, this is impossible to answer through scattering and is not adequately constrained through particle decays. But this question is of central importance to particle theory, as neutrinos allow unique tests of new physics. And it is of central importance to cosmology, as allowed νSI parameters would indicate the presence of strong self-scattering in the early universe. New techniques to probe νSI are needed, especially in the ν τ sector. These may lead to discoveries that should be tested in multiple ways, or to limits that will improve our abilities to search for physics beyond the standard models of particle physics and cosmology.
A way forward is to look for signatures of scattering with neutrinos in the CνB, which leads to characteristic features in the spectrum of astrophysical neutrinos at Earth [40][41][42][43]. This has become newly promising now that IceCube has detected O(100) events at energies larger than 60 TeV [47,48].
In this paper, we take advantage of the proposed Gen2 detector, develop a comprehensive theoretical treatment, and make predictions that include realistic experimental effects. We benefit from the improved knowledge of the neutrino mass spectrum: the measurements of the total neutrino mass in cosmology and of the neutrino mass ordering in the laboratory shape the experimental signatures that should be looked for. We also demonstrate that common approximations in the literature were unjustified in some important cases.
Our primary result is in Fig. 10: Gen2 will significantly improve the sensitivity to νSI, realizing the full potential of high-energy neutrino astronomy for testing νSI in propagation. It can probe a huge range of parameters where neutrinos do not free-stream as expected in the early universe. This includes the full parameter space relevant for the Moderately Interacting neutrino solution [25,39]. And, as discussed above, modifying the analysis inputs only induces slight sensitivity changes. At its best sensitivity, Gen2 will overcome laboratory constraints and become the strongest probe of neutrino selfinteractions between any flavors. The future quest for νSI discovery will not remain bounded to τ neutrinos as it is today.
We make our code publicly available so that our formalism can be straightforwardly applied to explore different parts of the parameter space with various sources. Some examples include the diffuse supernova neutrino background or cosmogenic neutrinos.
Should a signal appear at Gen2, there will be plenty of opportunities to test it. A key observable is the flavor composition (see, e.g., Ref. [79]). Another is exploiting point sources. The main purpose of Gen2 is to resolve individual neutrino sources [49]; any detection would be highly valuable to explore νSI [102]. The reason is that nearby sources should not be affected by νSI and could provide a better understanding of the high-energy astrophysical neutrino spectrum. The appearance of spectral signatures in far but not near sources would be a smoking gun for νSI. In addition, the scattering of neutrinos en route to the Earth could introduce measurable time delays [41]. In addition to the Gen2 optical array that we considered, the Gen2 radio array and other neutrino observatories can greatly extend the reach of our work by detecting higher-energy neutrinos [49,[95][96][97][98][99][100][101]. Finally, hints for νSI at Gen2 could leave signatures in astrophysical and cosmological observables.
As neutrino physics enters the precision era, the properties of these ghostly particles will be scrutinized better than ever. In this paper, we provide the necessary particle-physics framework as well as the phenomenological roadmap to make the most out of neutrino selfinteraction measurements in present and next-generation neutrino telescopes. Improvements in understanding high-energy astrophysical sources and further experimental sensitivity studies will enhance this progress. This will open a window into understanding whether neutrinos have sizable self-interactions, providing insight about physics beyond the standard model and the evolution of the early universe.
Under these assumptions, the differential scattering cross-section is given by [104]

Dirac neutrinos
If neutrinos are Dirac fermions, and assuming that there is the same amount of target neutrinos and antineutrinos, the total inclusive cross-section for initial mass eigenstate i and final mass eigenstate j is given by The first sum is the cross-section with two neutrinos in the final state, and the last term is the cross-section for double-scalar production, present only if s > 4M 2 φ . For double-neutrino production, we can classify the differential cross-section in terms of the visibility of the final state. Here, visible means that neutrinos must be left-handed and antineutrinos right-handed. We show the results in Table I. dσ dt k visible, l invisible 1 32π k visible, l visible 0  The differential cross-section is simpler than in the Dirac neutrino case, as all final states are observable. For