Present and future constraints on flavor-dependent long-range interactions of high-energy astrophysical neutrinos

The discovery of new, flavor-dependent neutrino interactions would provide compelling evidence of physics beyond the Standard Model. We focus on interactions generated by the anomaly-free, gauged, abelian lepton-number symmetries, specifically $L_e-L_\mu$, $L_e-L_\tau$, and $L_\mu-L_\tau$, that introduce a new matter potential sourced by electrons and neutrons, potentially impacting neutrino flavor oscillations. We revisit, revamp, and improve the constraints on these interactions that can be placed via the flavor composition of the diffuse flux of high-energy astrophysical neutrinos, with TeV-PeV energies, i.e., the proportion of $\nu_e$, $\nu_\mu$, and $\nu_\tau$ in the flux. Because we consider mediators of these new interactions to be ultra-light, lighter than $10^{-10}$ eV, the interaction range is ultra-long, from km to Gpc, allowing vast numbers of electrons and neutrons in celestial bodies and the cosmological matter distribution to contribute to this new potential. We leverage the present-day and future sensitivity of high-energy neutrino telescopes and of oscillation experiments to estimate the constraints that could be placed on the coupling strength of these interactions. We find that, already today, the IceCube neutrino telescope demonstrates potential to constrain flavor-dependent long-range interactions significantly better than existing constraints, motivating further analysis. We also estimate the improvement in the sensitivity due to the next-generation neutrino telescopes such as IceCube-Gen2, Baikal-GVD, KM3NeT, P-ONE, and TAMBO.

F Table of limits on the long-range potential 30 1 Introduction The high-energy astrophysical neutrinos discovered by the IceCube neutrino telescope [1][2][3][4][5][6][7][8] provide unprecedented potential to explore high-energy fundamental physics [8][9][10][11][12][13].Because they reach energies in the TeV-PeV range, they are sensitive to possible new physics at energy scales beyond the reach of terrestrial colliders.Because they travel Gpc-scale distances -essentially the size of the observable Universe -minute new-physics effects, αβ , with mass m ′ αβ , that mediates flavor-dependent long-range neutrino interactions.Our results are based on estimates of the measurement of flavor composition of high-energy astrophysical neutrinos in IceCube using 8 years of data [15] (through-going muons plus HESE, see Section 3) and on present-day uncertainties in the neutrino mixing parameters [16,17], assuming normal neutrino mass ordering.Limits are on neutrino-electron interactions, under the L e −L µ and L e −L τ symmetries, and on neutrino-neutron interactions, under the L µ −L τ symmetry.For the latter, we assume a Z ′ µτ -Z mixing strength of (ξ − sin θ W χ) = 5 × 10 −24 [18].Existing limits are from a recent global oscillation fit [19], atmospheric neutrinos [20], solar and reactor neutrinos [21], and non-standard interactions (NSI) [22][23][24].For comparison, we show the proofof-principle sensitivity (1σ) based on 2015 IceCube flavor-composition measurements [25] from Ref. [26].Indirect limits [27] are from black-hole superradiance (90% C.L.) [28] and the weak gravity conjecture [29], assuming a lightest neutrino mass of 0.01 eV.Figures 8 and D3 show projections for the year 2040; Fig. E4 shows results assuming inverted mass ordering.See Section 4.2 for details.Our results show that limits obtained from present-day flavor-composition measurements may improve significantly on existing ones.ordinarily undetectable, may accumulate en route to Earth and become detectable.(They also provide unique insight into the most energetic astrophysical sources [8,13,14].)We gear this potential to far-reaching studies of new, long-range neutrino interactions.
Presently, theory and experiment leave room for the possibility that neutrinos undergo interactions beyond the Standard Model weak interactions that would have gone undetected so far.In particular, new interactions that affect different neutrino flavors dif-ferently may modify neutrino oscillations, and so may be tested in oscillation experiments that use neutrinos from a variety of sources, natural and man-made.References [20,21] pioneered the idea of introducing such flavor-dependent interactions by gauging native global lepton-number symmetries of the Standard Model, generated by L α − L β , where L α (α = e, µ, τ ) is the α-flavor lepton number.Doing so introduces new neutrino-matter interactions mediated by a boson that, if light, makes it possible for distant matter to source a long-range potential that affects neutrino oscillation probabilities.Previous works have explored this possibility and placed constraints on the strength of the new interaction using atmospheric neutrinos [20], solar and reactor neutrinos [21], accelerator neutrinos [18,30], and a global fit to oscillation data [19].
Recently, Ref. [26] proposed a new way to probe these interactions using high-energy astrophysical neutrinos.Because the relative contribution of the standard oscillation mechanism, driven by the neutrino masses, falls with neutrino energy, at high energies the contribution from new interactions may become prominent and more easily detectable.By modifying oscillations -in an extreme case, by suppressing flavor transitions -the new interactions would also modify the flavor composition of high-energy astrophysical neutrinos, i.e., the proportion of ν e , ν µ , and ν τ in their flux.IceCube measures the flavor composition of the diffuse flux of high-energy astrophysical neutrinos [25,[31][32][33] (see also Refs.[34][35][36]) and, though the measurement is challenging, it is a versatile probe of highenergy neutrino physics [8,10,11,13,[37][38][39][40][41] and astrophysics [8,13,14].In particular, the sensitivity to new neutrino interactions stems from looking for characteristic deviations in the flavor composition relative to its expectation under standard oscillations.
Based on this, Ref. [26] provided a proof of principle for using the flavor composition measured by IceCube and forecast for its planned upgrade, IceCube-Gen2 [15], to constrain flavor-dependent long-range neutrino-electron interactions generated by gauging the L e −L µ and L e − L τ symmetries.It focused on ultra-light mediators, lighter than 10 −10 eV, which makes the interaction range ultra-long, from km to Gpc, depending on the mass.As a result, vast numbers of electrons in the Earth, the Moon, the Sun, the Milky Way, and in the distant Universe could source a sizable long-range potential.Yet, while the proof-ofprinciple prospects in Ref. [26] were promising, the statistical methods used therein are not amenable to being generalized to compute sensitivity beyond the 1σ statistical significance.
To address this, we revisit, revamp, and expand the tests of flavor-dependent long-range neutrino interactions based on flavor composition.There are four main improvements in our work.First, we adopt a Bayesian approach to compute constraints, based on the methods from Ref. [42].It allows us to combine consistently the uncertainties in the measurement of the flavor composition in IceCube and other neutrino telescopes, and the uncertainties in the neutrino mixing parameters that drive standard oscillations, and to compute constraints at arbitrary statistical significance.Second, not only do we consider neutrino-electron interactions induced by the L e − L µ and L e − L τ symmetries, as in Ref. [26], but also neutrino-neutron interactions induced by the L µ − L τ symmetry.Third, we study not only the sensitivity of IceCube and IceCube-Gen2, but also the sensitivity achieved by combining them with other neutrino telescopes expected by 2040, Baikal-GVD [43], KM3NeT [44], P-ONE [45], and TAMBO [46].Fourth, we use projected measurements of the mixing parameters in upcoming neutrino oscillation experiments, the Deep Underground Neutrino Experiment (DUNE) [47], Hyper-Kamiokande (HK) [48], and the Jiangmen Underground Neutrino Observatory (JUNO) [49], that will shrink the uncertainty on their values.
Figure 1 shows our main results, in the form of upper limits on the coupling strength of the new interactions, conservatively based on approximate present-day flavor-measurement capabilities in IceCube and of uncertainties in the mixing parameters.For mediators lighter than 10 −18 eV, our estimated upper limits confirm [26] the reach of flavor-composition measurements to test ultra-long-range interactions, now at higher statistical significance.We find that, already today, high-energy astrophysical neutrinos may be able to constrain flavor-dependent long-range neutrino interactions more strongly than existing limits.By 2040, our conservative forecast limits, not shown in Fig. 1, improve marginally (Section 4.2), though more power may be reaped from improvements to our analysis that could be afforded by significantly higher detection rates of high-energy neutrinos and precision in the mixing parameters.We point these out later (Section 5).
We provide our constraints as approximate, rather than as definitive, since they rely on specific, albeit realistic assumptions on present-day and upcoming experimental capabilities.For instance, they are obtained assuming that measurements of the flavor composition are performed on an astrophysical neutrino spectrum that is a power law ∝ E −2. 5 , where E is the neutrino energy.While this is compatible with IceCube observations [25], our results do not account for deviations from this assumption that lie within experimental error, and which could make the spectrum softer (see, e.g., Ref. [7]) or harder (see, e.g., Ref. [50]).Accounting for this would require deriving the IceCube flavor sensitivity for different choices of the neutrino spectrum, as in Ref. [25], which involves combining different data sets that are not available outside the Collaboration.Doing so lies beyond the scope of this paper, but we comment on possible future developments (Section 5).
This paper is organized as follows.Section 2 introduces long-range interactions and their effect on neutrino oscillations.Section 3 provides an overview of high-energy astrophysical neutrinos and shows the effect of long-range interactions on their flavor composition.Section 4 presents out statistical methods and results.Section 5 points out future improvements of our analysis.Section 6 summarizes and concludes.
2 New long-range, lepton-number neutrino interactions The Standard Model (SM), in addition to the gauge group SU (3) C × SU (2) L × U (1) Y , contains global U (1) symmetries associated to the baryon number and the three lepton numbers, L e , L µ , and L τ .While they cannot be gauged individually without introducing anomalies, certain combinations of them can be; for an extensive list, see Ref. [19].We focus on three well-motivated, anomaly-free symmetries that gauge lepton number differences [51][52][53]: U (1) Le−Lµ , U (1) Le−Lτ , and U (1) Lµ−Lτ .Each one introduces a new neutral gauge vector boson, Z ′ eµ , Z ′ eτ , and Z ′ µτ , that mediates new neutrino interactions [52][53][54] with electrons and neutrons, as we show below.(The Higgs sector differentiates between different lepton flavors [55]; however, here we focus only on the gauge interactions through the new boson.)The L α − L β gauge symmetries and their extensions have been explored in numerous scenarios, including as possible solutions of the Hubble tension [56], of the electron and muon (g − 2) anomalies [57][58][59][60], considering the new boson to be dark matter [61][62][63] and as radiation from compact binary systems [64], leptogenesis [65], their influence on muon-beam dumps at the TeV-scale [66], possible production of dark photons at the MUonE experiment [67], and explaining the electron and positron excess in the cosmic-ray flux [68,69].However, these studies consider mediators that are heavier than the ones we consider here, and couplings (see below) that are stronger.
The effective neutrino-matter interaction receives three contributions, mediated by the SM Z boson, by the new Z ′ αβ boson, and via the mixing between Z ′ αβ and Z, i.e., The first term in Eq. (2.1), L SM , is the SM contribution, i.e., where e is the unit electric charge, θ W is the Weinberg angle, l α and ν α are, respectively, the charged lepton and neutrino of flavor α = e, µ, or τ , u and d are the up and down quarks, and P L is the left-handed projection operator.(Equation (2.2), and also Eq. (2.4) below, assumes that matter is electrically neutral, i.e., that it contains equal numbers of electrons and protons [18]; this is also what we assume later when computing the new matter potential; see Section 2.2.) The second term in Eq. (2.1), L Z ′ , describes neutrino-matter interactions via the new mediator [18,19,52,53], i.e., where g ′ αβ (α, β = e, µ, τ , α ̸ = β) is the adimensional coupling constant associated to the L α − L β gauge symmetry.Because naturally occurring muons and tauons are scarce, we consider this contribution only under L e − L µ and L e − L τ , and the corresponding interaction to be sourced only by abundant electrons.
The third term in Eq. (2.1), L mix , contains terms that mix Z ′ αβ with Z [18, 70, 71], i.e., Ẑµ , where Ẑ′ µν and Bµν are the field strength tensors for U (1) ′ and U (1) Y , respectively, Ẑ and Ẑ′ are the gauge eigenstates corresponding to the neutral massive gauge bosons of the SM and the new U (1) ′ gauge symmetry, χ is the kinetic mixing angle, and δ M 2 is the squared-mass difference between Ẑ and Ẑ′ .Diagonalizing this Lagrangian redefines the fields in terms of physical states: the photon and two massive bosons, Z and Z ′ , that are related to Ẑ and Ẑ′ through a new mixing angle ξ [70], i.e., L ZZ ′ ⊃ (ξ − sin θ W χ)Z ′ µ Z µ .This introduces a four-fermion neutrino-matter interaction term via Z-Z ′ αβ mixing, i.e., where J ′ ρ = να γ ρ P L ν α − νβ γ ρ P L ν β and J ρ 3 = − 1 2 ēγ ρ P L e+ 1 2 ūγ ρ P L u− 1 2 dγ ρ P L d.To illustrate the effect of mixing, below we fix the value of the mixing strength to (ξ − sin θ W χ) = 5 × 10 −24 [18], the upper limit for a new interaction whose range is of the order of the distance between the Earth and the Sun [18,72,73].(The upper limit on the mixing strength for an interaction range of the order of the size of the Earth is slightly weaker and allows for larger mixing [18,72,73].)Because of our analysis choices, Z-Z ′ αβ mixing affects the results only under the L µ − L τ symmetry; we explain this below.
Figure 2 shows the Feynman diagrams corresponding to the above contributions.Regarding L SM , because the Z boson is heavy, the SM neutral-current interaction is shortrange and significant only in relatively high-density regions such as inside the Earth (which we account for later).Regarding L Z ′ , it induces new interactions between neutrinos and charged leptons.However, in practice, the scarcity of muons and tauons in Nature means that interactions are sourced only by electrons.As a result, this term in the Lagrangian is active only for the L e −L µ and L e −L τ symmetries.Regarding L mix , it induces interactions between neutrinos and electrons, protons, and neutrons.Under the L µ − L τ symmetry, this is the dominant contribution; its size depends on the value of the mixing strength (see above).Further, for a macroscopic collection of ordinary matter, like the ones we consider below, the contribution of electrons is nullified by the contribution of protons, leaving only neutrons to source the new interaction; see Appendix A for details.Under the L e − L µ and L e − L τ symmetries, the relative importance of L mix vs. L Z ′ depends on the value of the mixing strength; to be conservative, we assume no mixing for these two symmetries.
In summary, under the L e − L µ and L e − L τ symmetries, the new neutrino interaction is via L Z ′ and is sourced by electrons only.Under the L µ − L τ symmetry, it is via L mix and is sourced by neutrons only.Under all symmetries, L SM acts only inside the Earth.

Long-range neutrino interactions
A neutrino separated by a distance d from a collection of N e electrons, under L e − L µ and L e − L τ , or of N n neutrons, under L µ − L τ , experiences a Yukawa potential of where m ′ αβ is the mass of the Z ′ αβ boson and the coupling strength is The range of the interaction is ∼1/m ′ αβ ; beyond this distance from the source of the potential, it is suppressed due to the mediator mass.Depending on the interaction range, the potential receives contributions from nearby bodies -if the mediator is heavy -or also from distant bodies -if the mediator is light.We consider masses of 10 −35 -10 −10 eV, for which the interaction range is km-Gpc, which encompasses the Earth (⊕), Moon ( ), Sun  (⊙), Milky Way (MW) and the cosmological matter distribution (cos).Thus, the potential experienced by neutrinos is the sum of the contributions sourced by each of them, i.e., and the relative contribution of each term depends on the value of m ′ αβ .We do not compute the changing potential along the underground trajectories of the neutrinos inside the Earth or inside the Sun; see Ref. [19] for such treatment.Instead, like Ref. [26], we compute the average potential experienced by the neutrinos at their point of detection in IceCube.This approximation is especially valid for mediators lighter than about 10 −18 eV, for which the interaction range is longer than the Earth-Sun distance (see Fig. 1), and so all of the electrons and neutrons on the Earth or the Sun contribute to the potential experienced by a neutrino regardless of its position inside them.Below 10 −18 eV is also where we place limits in an unexplored range of mediator mass.
In what follows, we use Eq.(2.5) as a basis to calculate the potential induced by these sources.We adopt the methods introduced in Ref. [26] for the L e − L µ and L e − L τ cases.Below, we revisit them, extend them to the L µ − L τ case, and introduce refinements in the computation of the potential due to solar and cosmological electrons and neutrons.Throughout, we assume that matter is electrically neutral, so that the number of electrons and protons is the same, and that matter is isoscalar, so that the number of electrons and neutrons is the same, except for the Sun and for the cosmological distribution of matter.
The Earth.Astrophysical neutrinos travel inside the Earth from its surface to IceCube, located d IC = 1.5 km underground at the South Pole.Along the way, they undergo longrange interactions with underground electrons or neutrons.Under the L α − L β symmetry, the potential sourced by them is where θ z is the zenith angle along which the neutrinos travel, measured from the South Pole, the chord length traveled along this direction is , and the radius of the Earth is R ⊕ = 6371 km.
The average densities of electrons and neutrons along direction θ z are ⟨n e,⊕ ⟩ θz = ⟨n n,⊕ ⟩ θz , assuming matter is electrically neutral and isoscalar.To compute them, we adopt the matter density profile of the Preliminary Reference Earth Model [74], and we assume an electron fraction of Y e ≡ N e /(N e + N p ) = 0.5.The total number of electrons and neutrons in the Earth is N e,⊕ ≈ N n,⊕ ∼ 4 × 10 51 .Because we do not track the propagation of neutrinos inside the Earth and the changing long-range potential that they experience at different points along their trajectories, our treatment is approximate; yet, it allows us explore efficiently the parameter space.For a detailed treatment, see Ref. [19].
The Moon and the Sun.We treat the Moon and the Sun as point sources of electrons and neutrons.Under the L α − L β symmetry, the potential sourced by the Moon is where the distance between the Earth and the Moon is d ≈ 4 • 10 5 km, and the number of electrons and neutrons in the Moon is N e, = N n, ∼ 5 • 10 49 , assuming the lunar matter is electrically neutral and isoscalar.Similarly, the potential sourced by the Sun is where the distance between the Earth and the Sun is d ⊙ = 1 A.U., the number of electrons in the Sun is N e,⊙ ∼ 10 57 , and the number of neutrons in it is N n,⊙ = N e,⊙ /4.Like for Earth, we do not compute the changing long-range matter potential inside the Moon and Sun as neutrinos propagate inside them; for the latter, see Ref. [19].
The Milky Way.A neutrino of extragalactic origin traverses the Milky Way before reaching the Earth and may be affected by the long-range potential sourced by Galactic electrons and neutrons.We do not track the propagation of neutrinos inside the Milky Way.Instead, we estimate the effect of long-range interactions by computing the potential experienced by neutrinos at the location of the Earth by integrating the Galactic electron and neutron column densities across all possible neutrino trajectories that have the Earth as the endpoint.Under the L α − L β symmetry, this is × n e,MW (r, θ, ϕ) , for α, β = e, µ or e, τ n n,MW (r, θ, ϕ) , for α, β = µ, τ , (2.11) where the coordinate system is centered at the position of the Earth, 8.33 kpc away from the Galactic Center, and the densities of electrons and neutrons are n e,MW = n n,MW , assuming matter is electrically neutral and isoscalar.In the Milky Way, N e,MW ≈ N n,MW ∼ 10 67 electrons and neutrons are contained in stars and cold gas -distributed in a central bulge, a thick disc, and a thin disc -and in hot gas -distributed in a diffuse halo.Following Ref. [26], we adopt the "conventional model" from Ref. [75] for the matter density of the central bulge, thick disc, and thin disc, and the spherical saturated matter density from Ref. [76] for the diffuse halo; Fig. A1 in Ref. [26] shows the total matter distribution.
Cosmological electrons and neutrons.For m ′ αβ ≲ 10 −25 eV, the dominant contribution to the long-range potential is from the large-scale cosmological distribution of electrons and neutrons.We follow Ref. [26] to compute the potential, including the effect of the adiabatic cosmological expansion on the densities of electrons and neutrons; we defer to it for a derivation of the potential.While Ref. [26] assumed that the cosmological distribution is isoscalar, we instead account for the fact that, after the recombination epoch, the neutronto-proton ratio saturates to about 1/7 [77].Using this and assuming that cosmological matter is electrically neutral, the number of cosmological electrons is seven times higher than that of neutrons.Under the L α − L β symmetry, the potential at redshift z is × N e,cos (z) , for α, β = e, µ or e, τ N n,cos (z) , for α, β = µ, τ , where N e,cos (z) ≃ 7M H (z)/(8m p + 7m e ) is the number of electrons, m p and m e are the proton and electron mass, respectively, N n,cos (z) = N e,cos (z)/7 is the number of neutrons, M H is the baryonic mass inside casual horizon (see Eq. (16.105) in Ref. [78]), and d H is the size of the causal horizon, i.e., the radius of the largest sphere centered on the Earth within which events can be causally connected [79].We adopt a ΛCDM cosmology with Hubble constant H 0 = 100h km s −1 Mpc −1 , where h = 0.673 [80], the vacuum energy density Ω Λ = 0.692 and the matter density Ω M = 0.308 [81].Because we use the diffuse flux of high-energy astrophysical neutrinos, we account for the evolution with redshift of the number density of neutrino sources, ρ src , by averaging the potential over z, i.e., where V c is the comoving volume [82] and ρ src follows the star formation rate [83][84][85].
Figure 3 shows the total potential, Eq. (2.7), as a function of the mediator mass and coupling, for the L e − L β (β = µ, τ ) and L µ − L τ symmetries.To illustrate its behavior, we show the isocontour for the illustrative value of V αβ = 10 −18 eV (which is close to the final constraint that we obtain in Section 4.2).As the mass shrinks, the interaction range grows.As a result, the potential isocontour undergoes several step-like transitions, each of which represents the inclusion of the contribution of a more distant collection of electrons and neutrons into the total potential.From m ′ αβ ∼ 10 −10 eV to 10 −18 eV, the potential is dominated by the Earth, with a minor contribution from the Moon.The step at m ′ αβ ∼ 10 −18 eV represents the inclusion of the potential sourced by the Sun, which becomes dominant.Because the Sun contains far more electrons and neutrons than the Earth and the Moon, a smaller coupling is enough to achieve the same potential.The step at m ′ αβ ∼ 10 −27 eV represents the onset of dominance of Milky-Way electrons and neutrons, which far outnumber those in the Sun, and are concentrated in the Galactic Center.The last step, at m ′ αβ ∼ 10 −33 eV, represents the onset of the dominance of cosmological distribution of electrons and neutrons, which far outnumber those in the Milky Way.
For the L µ − L τ case, Fig. 3 shows that smaller couplings are needed to achieve the same illustrative value of the potential because it scales ∝ g ′ µτ , rather than ∝ g ′ 2 eβ , as in the L e − L µ and L e − L τ cases; see Eq. (2.6).In our work, including in Fig. 3, for the L µ − L τ case we fixed the Z ′ µτ -Z mixing strength, (ξ − sin θ W χ), to its maximum allowed value (see Section 2.1).Decreasing or increasing the mixing strength would shift the potential isocontour in Fig. 3 up or down, respectively.

Neutrino flavor-transition probabilities
In the presence of the new L α − L β gauge symmetry, the Hamiltonian that drives neutrino propagation, written in the flavor basis, is (2.14) The contributions on the right-hand side describe, respectively, oscillations in vacuum, standard neutrino-matter interactions, and the new neutrino-matter interactions.
Oscillations in vacuum are driven by where E is the neutrino energy, ∆m 2 31 ≡ m 2 3 −m 2 1 , and ∆m 2 21 ≡ m 2 2 −m 2 1 are the atmospheric and solar mass-squared splittings, respectively, m i the mass of the ν i mass eigenstate, and U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix, parameterized as a product of three rotation matrices, i.e., R(θ 23 )R(θ 13 , δ CP )R(θ 12 ), where θ ij (i, j = 1, 2, 3) are three mixing angles and δ CP is the Dirac CP-violating phase.Under normal mass ordering (NMO), where m 1 < m 2 < m 3 , the present-day best-fit values of the mixing parameters from the global oscillation analysis of Refs.[16,17] 4 and D1, we adopt these values to illustrate the effect of the new interaction on the neutrino flavor-transition probabilities.Later, in our statistical analysis, we let their values float within their present-day and future predicted allowed ranges.
In Eq. (2.14), the contribution of the SM potential from coherent forward neutrino scattering on electrons is where V CC = √ 2G F ⟨n e,⊕ ⟩ θz is the charged-current neutrino-electron interaction potential and G F is the Fermi constant.In our calculations, this contribution is relevant only inside Earth, where matter densities are high.The potential above is for neutrinos; for antineutrinos, it flips sign, i.e., V mat → −V mat .
Finally, in Eq. (2.14) the contribution from the new matter interaction is where V αβ is the potential, calculated using Eq.(2.7).Depending on the mediator mass, it is sourced by nearby and faraway sources of electrons and neutrons.The potential above is for neutrinos; for antineutrinos, it flips sign, i.e., V αβ → −V αβ .The ν α → ν β flavor-transition probability associated to the Hamiltonian, Eq. (2.14), is where L is the distance traveled by the neutrino from its point of production to the Earth, ∆ m2 ij ≡ m2 i − m2 j , with m2 i /2E the eigenvalues of the Hamiltonian, and U m is the matrix that diagonalizes the Hamiltonian, parameterized like the PMNS matrix, but in terms of new mixing parameters θ m 23 , θ m 13 , θ m 12 , δ m CP .In Appendix B, we show the evolution of the three modified mixing angles with the long-range potential.The flavor-transition probability is oscillatory, but for high-energy neutrinos that travel cosmological-scale distances, like the ones we consider, the oscillations are rapid, i.e., ∆ m2 ij L/(2E) ≫ 1.Given the limited energy resolution of neutrino telescopes, they are only sensitive to the average probability, We use this expression to produce all our results below.Under standard oscillations, i.e., when V αβ = 0, the average oscillation probability loses its dependence on neutrino energy.In contrast, in the presence of the long-range interactions, it retains the energy dependence via the interplay of the vacuum and potential contributions to the Hamiltonian, Eq. (2.14).Figure 4 shows average oscillation probabilities as functions of the potential V αβ , computed for a fixed neutrino energy of 100 TeV, representative of high-energy astrophysical neutrinos.In what follows, we focus on the probabilities under the L e − L µ symmetry.Figure 4 shows that for low values of the potential, i.e., V αβ ≪ ∆m 2 ij /(2E), oscillations are standard.Deviations start to appear when the potential becomes comparable to the Hamiltonian in vacuum, i.e., V αβ ≈ 10 −20 eV.For neutrinos, sharp features -dips or peaks -appear at V eµ ≈ 10 −17 eV in the ν e → ν β (β = e, µ, τ ) probabilities, reflecting that θ m 13 ≈ 45 • (see Appendix B), which maximizes U m e3 and makes flavor transitions resonant.At V eµ ≈ 10 −19 eV, less prominent features reflect that θ m 12 = 45 • ; they are less prominent because in vacuum θ 12 ≈ 45 • already.Similar features appear in the ν µ → ν τ and ν τ → ν τ probabilities, for analogous reasons.For large values of the potential, i.e., V αβ ≫ 10 −17 eV, the term V αβ dominates the Hamiltonian and, because it is diagonal, it suppresses flavor transitions.For L e − L τ , shown in Fig. D1, the flavor-transition probabilities are similar to L e − L µ .Figure 4 shows that, for L µ − L τ , no resonance due to θ m 13 occurs neither for neutrinos nor antineutrinos because θ m 13 never reaches 45 • , but there is a small dip or jump in the probabilities around V µτ ≈ 10 −18 eV, due to resonant θ m 12 .In the main text, we produce results assuming normal neutrino mass ordering, motivated by recent hints from global oscillation fits [16,17,[86][87][88].For antineutrinos, the sharp features above do not appear because θ m 13 and θ m 12 never become resonant under normal neutrino ordering.However, they do appear under inverted mass ordering, because ∆m 2  31 is negative; see Appendix E.
Inside the sources, the high-energy protons interact with surrounding matter, in protonproton (pp) collisions [127][128][129], or radiation, in proton-photon (pγ) collisions [128,[130][131][132].These generate high-energy pions that, upon decaying, make neutrinos, i.e., π + → µ + + ν µ , followed by µ + → νµ + e + + ν e , and their charge-conjugated processes.(At higher energies, the decay of kaons and other production channels become important [130][131][132][133][134].) Each neutrino carries, on average, 5% of the energy of the parent proton [129,130].In production via pp collisions, the neutrino spectrum is a power law ∝ E −γ that follows that of the parent protons, with the value of γ loosely expected to be in the interval [2,3].In production via pγ collisions, its shape depends on those of the parent protons and photons, and it is bump-like, on account of the photon spectrum peaking at a characteristic energy.Because the diffuse neutrino flux results from the addition of the emissions of all neutrino sources, the above features in the neutrino spectrum are softened.So far, IceCube observations are described well by an unbroken power with the value of γ dependent on the data set used and with little to no preference for alternative spectral shapes [7,50], though this may change with further observations [135].

Neutrino detection and flavor composition
Presently, IceCube is the largest neutrino telescope [6,136]: it consists of 1 km 3 of deep Antarctic ice instrumented with photomultipliers.In it, high-energy neutrinos interact with nucleons in the ice, mainly via deep inelastic scattering [32,[137][138][139], either neutral-current (NC), i.e., ν l + N → ν l + X, where l = e, µ, τ and X are final-state hadrons, or chargedcurrent (CC), i.e., ν l +N → l+X.Final-state charged particles shower and emit Cherenkov radiation, which propagates through the ice and is detected by the photomultipliers.From the number of photons detected and from their spatial and temporal distributions, IceCube reconstructs the energy, direction, and flavor of the interacting neutrino [136,140].
Broadly stated, IceCube detects three types of neutrino-induced events: tracks, cascades, and double cascades.Tracks are made mainly by CC interactions of ν µ , where the final-state hadrons initiate a shower around the interaction vertex, and the final-state muon leaves a track of Cherenkov light in its wake, several kilometers in length and easily identifiable.Tracks are also made in 17% of ν τ CC interactions where the final-state tauon decays into a muon [141].Cascades are made by CC interactions of ν e and ν τ , where both the final-state lepton and hadrons shower around the interaction vertex, and by NC interactions of neutrinos of all flavors, where only the final-state hadrons shower.Double cascades are made by CC interactions of ν τ in which the final-state hadrons make a first shower, centered around the interaction vertex.The final-state tauon is energetic enough to decay some distance away, generating a second shower [33,142,143].
Starting events are those where the neutrino interacts within the detector volume.Because a large fraction of the neutrino energy is deposited in the ice, starting events trace neutrino energy closely [140].A subset of them above 60 TeV are known as High-Energy Starting Events (HESE) [7]: they are subjected to a self-veto [144,145] that reduces the contamination of atmospheric neutrinos and muons and, as a result, they have the highest content of neutrinos of astrophysical origin.However, its detection rate is relatively low, of approximately 10 events per year [7].Through-going tracks are those where ν µ (and ν τ ) interact outside the detector and produce muon tracks that cross part of it [50].Their detection rate is orders-of-magnitude higher than for HESE, but, because they have lower energies, they are dominated by atmospheric neutrinos and muons, made in cosmic-ray interactions in the atmosphere.Further, for tracks, because the location of the neutrino interaction is unknown, the neutrino energy can only be inferred uncertainly [140].In our analysis below, we exploit the combined capabilities of HESE and through-going events for flavor measurements, motivated by Ref. [25].
Because different event types can be made by more than one neutrino flavor, a straightforward, event-by-event correspondence between detected event type and neutrino flavor is typically unfeasible, except for double cascades, which are only made by ν τ .(In addition, showers made by the Glashow resonance [146], recently discovered by IceCube [147], are triggered exclusively by the interaction of 6.3-PeV νe with electrons [148][149][150][151].) Still, from the relative number of detected events of different types, it is possible to infer statistically the flavor composition of the neutrino flux, i.e., the proportion of ν e , ν µ , and ν τ in it, even if, because of the above limitations, the measurement uncertainties are significant.Our analysis below rests on the capabilities of IceCube and upcoming neutrino telescopes to infer the flavor composition in such a way, following the same spirit as Refs.[42,135,152] The flavor composition can be inferred either using only HESE showers, tracks, and double cascades, or by combining them with through-going muons.IceCube has reported the flavor composition in several analyses: using the first 3 years of HESE data [31], a combination of 4 years of contained events plus 2 years of through-going tracks [25], 5 years of contained events starting at lower energies [32], and, most recently, 7.5 years of HESE data, including a dedicated search for double cascades [33].Independent analyses have reported complementary results; see, e.g., Refs.[34-36, 153, 154].The precision of these analyses is limited by the relative scarcity of HESE events.Presently, the most precise measurements come from the combined analysis of Ref. [25], where the large number of through-going tracks pins down the flavor content of ν µ .Such a combined analysis has not been updated since 2015 [25] in spite of the availability of larger event samples, on account of its complexity.Our analysis below is based on realistic estimates for what such analysis could look like with present and future data; see Section 3.4.
For a given choice of flavor composition at the sources, en route to Earth neutrino oscillations modify it into where the average oscillation probability, Pβα , is computed using Eq.(2.19).During propagation, the total number of neutrinos is conserved; flavor transitions simply redistribute them into different flavors.Figure 5 shows the allowed regions of flavor composition at Earth under standard oscillations [42,159], i.e., when the potential V αβ = 0, for the three scenarios, and varying the values of the standard mixing parameters within their presentday 1σ allowed ranges [16,17].In the presence of long-range interactions, i.e., when V αβ is nonzero, because the flavor-transition probabilities are modified, so is the flavor composition at Earth.In this case, the flavor composition becomes energy-dependent, via the probability, in contrast to the flavor composition computed under standard oscillations.Figure 5 illustrates the expected flavor composition at Earth in the presence of the new matter potential, under the L e − L µ and L µ − L τ symmetries, for a representative neutrino energy of 100 TeV, and for the three benchmark scenarios of flavor composition at the sources.Results under the L e −L τ symmetry are similar to those under L e −L µ ; see Fig. D2.In Fig. 5, we compute f α,⊕ for varying values of the long-range matter potential, V αβ , and of the standard mixing parameters, within their 1σ allowed ranges [16,17].For small values of the potential, i.e., V αβ ≲ ∆m 2  21 /(2E) ∼ 10 −20 eV, the flavor composition is close to the expectation from standard oscillations.As the potential grows, the behavior of the flavor composition traces that of the flavor-transition probabilities (Figs. 4 and D1): the wiggles in f α,⊕ reflect the resonant features in the probabilities.For large values of the potential, i.e., V αβ ≳ 10 −16 eV, the new matter potential becomes the dominant contribution to the Hamiltonian, Eq. (2.14).In that case, because the contribution of the new potential is flavor-diagonal, flavor transitions are suppressed (Section 2.3).Accordingly, Fig. 5 shows that the flavor composition at Earth is the same as at the sources, i.e., f α,⊕ ≈ f α,S .

Analysis choices
Flavor composition at the sources independent of energy.Because different neutrino production channels become available at different energies [129,130,132,134], the flavor composition emitted by the sources, f α,S , may vary with energy [156,157,159,[163][164][165], contingent on key source properties like the magnetic field intensity [166,167].Nevertheless, adopting the same practice as IceCube flavor measurements [25,[31][32][33], we take f α,S to be constant within the energy range of our analysis.This is justified by the limited number of HESE events, which would be further diluted by attempting to measure the flavor composition in multiple energy bins.Under this assumption, any energy dependence that the flavor composition at the Earth, f α,⊕ , might have is due solely to the presence of the new neutrino-matter interaction.Reference [158] explores the interplay between energy dependence that stems from neutrino production and from new physics.
Averaging f α,⊕ between ν and ν.Because, in high-energy neutrino telescopes, events triggered by neutrinos and antineutrinos are so far indistinguishable from one another, we average the flavor composition at Earth between neutrinos and antineutrinos, assuming they are present in equal proportions in the flux that reaches Earth.In actuality, their proportions are practically unknown (see, however, Refs.[147,168]), but the above assumption aligns with theoretical expectations, especially from multi-pion neutrino production at high energies; see, e.g., Ref. [132].As a result of averaging, the resonant features in the oscillation probabilities, which are more prominent for neutrinos than for antineutrinos (Fig. 4), appear softened in the flavor composition at Earth that we use.
Estimates of present and future flavor sensitivity.In our analysis, we use estimates of the present-day flavor sensitivity of IceCube and of future detectors; see Fig. 5.These are based on estimated combined analyses of HESE events plus through-going tracks (see Section 3.2), motivated by Ref. [25], which offer the greatest sensitivity.We adopt the flavor measurement projections for IceCube-Gen2 from Ref. [15], which were generated assuming a plausible neutrino spectrum ∝ E −2.5 .Because these projections were only generated assuming a flavor composition at the sources of (1/3, 2/3, 0) S , later we only derive constraints on long-range neutrino interactions for that one benchmark scenario.As in Ref. [42], we isolate the contributions of IceCube and IceCube-Gen2 to the flavor sensitivity from Ref. [15].This allows us to generate results based on estimates of the present-day flavor sensitivity, for the year 2020, using 8 years of IceCube, and forecasts for the year 2040, using 15 years of IceCube plus 10 years of IceCube-Gen2, and using the additional contribution of future neutrino telescopes Baikal-GVD [43], KM3NeT [44], P-ONE [45], and TAMBO [46].Following Ref. [42], we assume that future detectors will have the same detection efficiency of HESE and through-going tracks as IceCube (see Section 4), We show results for three benchmark choices of the flavor composition at the sources, (f e,S , f µ,S , f τ,S ) = (1/3, 2/3, 0), the canonical expectation, (0, 1, 0), and (1, 0, 0), compared to the expectation from standard oscillations ("Std.osc.").The flavor composition is averaged between neutrinos and antineutrinos, assuming they exist in equal proportion in the flux.For this plot only, we fix the neutrino energy to 100 TeV, assume normal neutrino mass ordering, and vary the standard mixing parameters within their 1σ allowed ranges; our statistical methods systematically vary these choices.Flavor-composition predictions are compared against three estimates of the flavor sensitivity of neutrino telescopes using 8 years of IceCube data ("IC 8 yr"), 15 years of IceCube data plus 10 years of IceCube-Gen2 ("IC 15 yr + Gen2 10 yr"), and the combined exposure of all neutrino telescopes available by 2040 ("Combined ν telescopes").For L e − L τ , results are similar as for L e − L µ ; see Appendix D. See Section 3.3 for details.though with different rates, and rescale the IceCube flavor sensitivity to their respective expected exposures.This is admittedly a simplification, made necessary by the current absence of details on the capabilities of future detectors.

Statistical analysis and results
Our goal is to compute the bounds on the long-range-interaction parameters -the mediator mass, m ′ αβ , and the coupling, g ′ αβ -that can be placed by measuring the flavor composition.To do this, we contrast our predictions of the flavor composition obtained under long-range interactions (Section 3.3) against the estimated present and future flavormeasuring capabilities of IceCube and other upcoming detectors (Section 3.4).First, we compute bounds on the long-range matter potential, V αβ ; then, we translate them into bounds on the mediator mass and coupling.

Statistical analysis
We adopt the Bayesian statistical methods introduced in Ref. [42].For given flavor ratios at the source, f α,S , we assume one of the three benchmark scenarios outlined in Section 3.

+ JUNO + DUNE + HK
Table 1: Observation epochs of our analysis: years 2020 and 2040.For each epoch, we show the neutrino telescopes that measure the flavor composition of high-energy astrophysical neutrinos and the oscillation experiments that constrain the value of the neutrino mixing parameters.We assume that upcoming neutrinos have flavor-measuring capabilities similar to those of IceCube.In 2040, the combined telescopes include IceCube, IceCube-Gen2 [15], Baikal-GVD [43], KM3NeT [44], P-ONE [45], and TAMBO [46].See Section 3.2 for details.For our 2040 projections of the measurement of mixing parameters, we assume that their real values are the present-day best-fit values from the NuFit 5.1 global oscillation fit [16,17].See Section 2.3 for details.
ϑ ϑ ϑ ≡ (s 2 12 , s 2 23 , s 2 13 , δ CP ), with s ij ≡ sin θ ij , we compute the associated energy-averaged flavor composition at Earth, ⟨f f f ⊕ ⟩ ≡ (⟨f e,⊕ ⟩, ⟨f µ,⊕ ⟩, ⟨f τ,⊕ ⟩), using the procedure introduced in Section 3.4.We assess the compatibility of these predictions with measurements of the flavor composition in neutrino telescopes through three factors: π(V αβ ), the prior associated to the value of V αβ ; π(ϑ ϑ ϑ), the prior associated to the value of ϑ ϑ ϑ; and L(⟨f f f ⊕ ⟩), the likelihood of having measured the flavor composition ⟨f f f ⊕ ⟩.We expand on these factors below.Using them, we compute the posterior probability density of V αβ , marginalized over all possible values of the mixing parameters, i.e., This method represents an improvement over the one used in the original study of longrange interactions using the flavor composition, Ref. [26].There, the sensitivity to longrange interactions was derived from a straightforward comparison of the predicted flavor composition vs. the flavor sensitivity of IceCube, and the method was not amenable to being generalized to compute sensitivity beyond the 1σ statistical significance.Here, in contrast, the Bayesian approach allows us to derive limits at higher statistical significance and to properly account for the prior on mixing parameters.Table 1 summarizes the two epochs for which we perform the above analysis: the present, represented by the year 2020, and the future, represented by the year 2040.Each epoch has an associated precision on mixing parameters and flavor measurements associated to it, encoded in the prior π(ϑ ϑ ϑ) and the likelihood L (f f f ⊕ ).We expand on them below.
Prior on the long-range matter potential, π(V αβ ).To avoid introducing unnecessary bias, we use a uniform prior on V αβ in the interval 10 −24 -10 −16 eV.This interval covers the full range of long-range-interaction effects: towards the lower end, standard oscillations dominate, and towards the upper end, long-range interactions dominate.
Prior on the mixing parameters, π(ϑ ϑ ϑ).We take the prior of each mixing parameter, today and in the future, to be centered at its present-day best-fit value from the NuFit 5.1 [16,17] global oscillation fit.The distributions for s 2  12 and s 2 13 are each normal and uncorrelated with other parameters, while those of s 2 23 and δ CP are correlated.For our 2020 estimates, we use the NuFit 5.1 χ 2 distributions to compute the joint likelihood, −2 ln L 2020 ≡ χ 2 (s 2 12 ) + χ 2 (s 2 13 ) + χ 2 (s 2 23 , δ CP ).For our 2040 projections, we combine this with future measurements of θ 12 by JUNO, as computed in Ref. [42], and of θ 23 and δ CP by DUNE [47] and HK [48], which we generate ourselves.We do so in dedicated simulations using GLoBES [169], and adopting the detector descriptions for DUNE [170,171], using 3.5 years of runtime in neutrino and antineutrino modes each (using a revised plan of 5 years for each would not change results significantly [172]), and for HK [48], using 2.5 years of runtime in neutrino mode and 7.5 years in antineutrino mode.(For θ 23 and δ CP , our projections differ from those of Ref. [42] because we generate them using the best-fit values from NuFit 5.1 [16,17] instead of NuFit 5.0, as a result of which the value of θ 23 has shifted from the upper octant to the lower octant [173].)For s 2  13 , we keep the width of its prior fixed to its present-day size, which is dominated by measurements in Daya Bay [174], as they are not expected to be improved upon significantly by 2040.In summary, for 2040, we use the likelihood −2 ln L 2040 ≡ −2 ln L 2020 + E χ 2 E , where the contribution of DUNE, HK, and JUNO is each χ 2 E = i,j (ϑ i − θi )Σ −1 E,ij (ϑ j − θj ), Σ ij is the covariance matrix for the mixing parameters ϑ i , ϑ j , and θi , θj are their assumed real values.
Likelihood of flavor measurement, L (⟨f f f ⊕ ⟩).The likelihood of flavor measurement represents the precision with which neutrino telescopes can measure the flavor composition.For the 2020 and 2040 epochs, we adopt the same flavor likelihoods as Ref. [42].By construction, they are centered on the canonical flavor composition at Earth, near (1/3, 1/3, 1/3) ⊕ , that is expected from neutrino production via the full pion decay chain (Section 3.1).Hence, the bounds on long-range interactions that we derive later are obtained under the assumption that the true value of the measured flavor is such.For our 2020 estimates, we adopt an estimate of the IceCube flavor sensitivity that would be obtained using 8 years of HESE events and through-going tracks [15], assuming a flux with spectral index γ = 2.5.For our 2040 projections, we adopt either the expected flavor sensitivity from 15 years of IceCube plus 10 years of IceCube-Gen2 measurements, extracted from Ref. [15], or the combined sensitivity from that plus the expected 2040 sensitivity of Baikal-GVD [43], KM3NeT [44], P-ONE [45], and TAMBO [46], extracted from Ref. [42], which we defer to for details.As illustration, Fig. 5 shows 95% C.L. contours for the 2020 and 2040 flavor-measurement likelihoods.

Results
Figure 6 shows the resulting posteriors of V eµ , under the L e − L µ symmetry, and of V µτ , under the L µ − L τ symmetry, computed using Eq.(4.1), for the years 2020 and 2040.Results for V eτ , under the L e − L τ symmetry, are similar to V eµ ; see Appendix D. The posteriors are maximum at a potential of roughly 10 −19 eV.Because the posteriors plateau at lower values, we can place upper limits on the values of the potentials.Values above roughly 10 −18 eV are noticeably less favored, and they become more so when moving from 2020 to 2040 due to improvements in the precision of the mixing parameters and flavor composition, i.e., to narrower π(ϑ ϑ ϑ) and L (f f f ⊕ ), respectively.We assume a measurement of the flavor composition at Earth centered around the canonical expectation of about (1/3 : 1/3 : 1/3) ⊕ , corresponding to f S = (1/3 : 2/3 : 0) at the sources; see Fig. 5. See Fig. 7 for the limits on the potential derived from the posterior.Results for V eτ , under the L e − L τ symmetry, are similar to V eµ ; see Appendix D. We assume that all neutrino telescopes available by 2040 have detection efficiencies similar to that of IceCube (Section 3.4).Results here are obtained assuming normal neutrino mass ordering (NMO); results for inverted ordering are in Appendix E. See Section 4.1 for details.
Figure 7 (also Table F1) shows the resulting upper limits on the long-range potentials.
A limit has value of V ⋆ αβ such that the integral of the posterior, V ⋆ αβ 0 P(V αβ )dV αβ , equals a desired credible level (C.L.).We show limits at the 95% C.L in Fig. 7 (also in Figs. 1, 8, and  D3).Differences in the limits obtained under different symmetries are moderate because, under the assumption of (1/3, 2/3, 0) S , their predicted flavor compositions at Earth are only moderately different; see Figs. 5 and D2.
Figure 7 reveals modest improvement in the limits between 2020 and 2040, of a factor of 2.5-3, and marginal improvement in 2040 between using only IceCube plus IceCube-Gen2 and using all future available telescopes.Given the significant increase in the combined neutrino detection rate between 2020 and 2040 (see Fig. 1 in Ref. [42]), and between using only IceCube plus IceCube-Gen2 and all telescopes, this is unexpected at first glance.Yet, close inspection of the variation of the flavor composition at Earth with the longrange potential in Fig. 5 reveals the subtle reason.For the canonical choice of flavor composition at the sources of (1/3, 2/3, 0) S , which we have adopted to obtain the limits on the long-range potential, the flavor composition at Earth approaches the center of the flavor triangle -where the likelihood of flavor measurement is largest and from whence the limit is predominantly derived -from the direction along which the gradients of the 2020 and 2040 likelihoods resemble each other the most.For comparison, we show the standard-oscillation contribution ("Std.osc.") to the Hamiltonian, ∆m 2 21 /(2E), evaluated at the present-day best-fit value of ∆m 2 21 [16,17] and at different values of the neutrino energy, E. We assume that all neutrino telescopes available by 2040 have detection efficiencies similar to that of IceCube (Section 3.4).Figure 8 shows these limits translated into limits on the coupling, g ′ αβ .See Appendix F for the values of the limits.See Section 4 for details.
Figure 8 shows how the upper limits on the potential translate into upper limits on the coupling, g ′ αβ , as a function of the mediator mass, m ′ αβ .In Fig. 8, each limit is an isocontour of potential in the g ′ αβ -m ′ αβ plane; we use Eq.(2.7) to translate the limits on V αβ in Fig. 7 into the limits in Fig. 8. Thus, the step-like transitions in Fig. 8 have the same origin as those of the potential itself (Fig. 3): the limits on g ′ αβ are stronger for smaller mediator masses, where the potential is due to a larger number of electrons or neutrons.Figure 1 spotlights our 2020 projections combining all neutrino telescopes.
For L e − L µ (and also for L e − L τ , see Appendix D), Fig. 8 reveals that already our 2020 estimate is better, by about one order of magnitude, than the proof-of-principle sensitivity based on 2015 IceCube flavor measurements from Ref. [26].(However, a fair comparison is not possible due to differences in their statistical significance and methods.).For L µ − L τ , we estimate and forecast limits based on the flavor composition for the first time.The limits for g ′ µτ shown in Fig. 8 were obtained by fixing the mixing strength to (ξ − sin θ W χ) = 5 × 10 −24 (see Section 2.1).Because V µτ ∝ g ′ µτ (ξ − sin θ W χ) (see Eqs. (2.5) and (2.6)), given an upper limit on V µτ , a smaller value of the mixing strength, which is plausible, would entail a weaker limit on g ′ µτ .(Alternatively, the limits for L µ − L τ in Fig. 8 can be interpreted as being on the product of the coupling times the mixing strength, without any assumption on the value of the latter.) Figure 8, and also Figs. 1 and D3, show that, for L e − L µ and L e − L τ , our limits, both estimated for 2020 and forecast for 2040, improve on existing ones obtained from atmospheric neutrinos [20], solar and reactor neutrinos [21], reinterpretations [27] of bounds  1, but now showing also the projected 2040 limits using either IceCube plus IceCube-Gen2, or those combined with Baikal-GVD, KM3NeT, P-ONE, and TAMBO.Left: Limits on neutrino-electron interactions under the L e − L µ gauge symmetry.Right: Limits on neutrino-neutron interactions under the L µ − L τ gauge symmetry; existing limits from accelerator neutrinos in MINOS (95% C.L.) are from Ref. [18].See Section 4 for details.
on the coefficients of non-standard neutrino interactions (NSI) [19,[22][23][24] (see Appendix C), and the recent global oscillation fit from Ref. [19].For L µ − L τ , light mediators in the mass range that we consider are largely unconstrained, except towards high masses, from accelerator neutrinos [18] and NSI.
In summary, our results show that, already today, the estimated flavor sensitivity of IceCube and the uncertainty in the mixing parameters have the potential to significantly improve the constraints on flavor-dependent long-range neutrino interactions.

Future improvements
The goal of our analysis is to estimate the reach of flavor composition of high-energy astrophysical neutrinos to constrain flavor-dependent long-range neutrino interactions rather than providing detailed results based on real flavor measurements.In doing so, we incorporated assumptions in our analysis towards facilitating our estimates.Below, we list them and point out with envisioned improvements.
Alternative flavor composition choices.To produce our limits, we limited ourselves only to the case where the flavor composition at Earth is the canonical expectation of (1/3, 1/3, 1/3) ⊕ , from neutrino production via the full pion decay production, (1/3, 2/3, 0) S .This was due to the unavailability of the projected flavor sensitivity of neutrino telescopes to alternative choices of flavor composition (Section 3.4), but our methods are general and applicable also in those alternatives.We encourage experimental collaborations to make publicly available the flavor sensitivity of their experiments under assumptions other than the canonical one.
Flavor-measurement capabilities of upcoming detectors.Our forecasts for upcoming neutrino telescopes -Baikal-GVD, IceCube-Gen2, KM3NeT, P-ONE, TAMBO -assume that their capabilities to measure the flavor composition will be similar to those of IceCube.This is a necessary assumption given the absence of their realistic capabilities at the time of writing.However, it is foreseeable already that differences in their geometries -most notably, in the spacing between photomultiplier strings -will imply differences between their detection capabilities.Further, we have not considered the possibility of proposed improvements in flavor-tagging techniques, like those provided by muon and neutron echoes [175].As the specific flavor-measurement capabilities of different upcoming experiments become better defined, future versions of our analysis will be able to incorporate them.
Using the energy dependence of the flavor composition.Our results are based on the comparison of the energy-averaged flavor composition at Earth vs. the flavor composition measured over an energy interval, assuming that it is constant within it (Section 3.4).However, long-range interactions modify the flavor composition in an energy-dependent way, so there may be additional probing power to be reaped from measuring the flavor composition in multiple energy bins.Present-day event rates are likely insufficient to do this, but IceCube-Gen2 may be able to [15].Further, doing this could allow us to distinguish between the variation of the flavor composition with energy that stems from long-range interactions from the one that stems from the neutrino production mechanism changing with energy, e.g., from full pion decay chain at low energies to muon-damped at high energies [157][158][159][163][164][165][166][167].
Astrophysical uncertainties.To produce our results, we assumed that the energy spectrum of the high-energy astrophysical neutrinos is a power law ∝ E −γ , with fixed γ = 2.5, compatible with IceCube analyses that combine HESE events and throughgoing tracks [25].However, the value of γ is only uncertainly known and varies somewhat depending on the data set used to measure it [7,50].Further, the Ice-Cube observations admit a description with differently shaped energy spectra; e.g., Refs.[7,50,117,135,176,177].Finally, because the diffuse flux of neutrinos is the addition of contributions from sources located at different distances, it depends on how the number density of sources evolves with redshift, which we have neglected in our analysis.Incorporating uncertainties in the shape of the neutrino spectrum and the redshift distribution of sources would likely weaken the bounds that we report, though a full analysis is beyond the scope of this paper.
Computing neutrino propagation.When computing the long-range matter potential (Section 2.2), we do so at the location where the neutrinos are detected, at Ice-glance, this should reduce sensitivity to long-range interactions.Surprisingly, we find that using the estimated current flavor sensitivity of IceCube and current mixing parameter uncertainties, high-energy astrophysical neutrinos could tightly constrain long-range interactions, surpassing existing limits (see Fig. 1).
The discovery of a new fundamental interaction would mark groundbreaking progress.Our findings show that high-energy astrophysical neutrinos, already today, have the potential to probe this more rigorously than present-day constraints.

E Results assuming inverted neutrino mass ordering
In the main text, we show mainly results obtained assuming normal neutrino mass ordering.Below, we show results assuming inverted ordering.There are differences in the oscillation probabilities: under inverted ordering, resonant features are prominent for antineutrinos rather than for neutrinos.However, because we average the flavor composition at Earth between neutrinos and antineutrinos, the limits on the long-range matter potentials that we obtain (Fig. 7) are largely insensitive to the mass ordering.
Figures E1 and E2 show the oscillation probabilities as functions of the long-range matter potentials.
Figure E3 shows the corresponding neutrino flavor composition at Earth. Figure E4 shows the posterior probability distribution of the long-range potentials, obtained following the procedure described in Section 4.1.
Figure E5 shows the corresponding upper limits on the couplings, g ′ αβ .

F Table of limits on the long-range potential
Table F1 shows the upper limits on the long-range matter potentials obtained in Section 4.2.

2 2 . 1 25 A 26 B 27 C 27 D 28 E
New long-range, lepton-number neutrino interactions 4 New lepton-number symmetries: L e − L µ , L e − L τ , and L µ − L τ 4 Interaction potential under the L µ − L τ symmetry Evolution of the modified mixing angles with long-range potential Reinterpretation of Super-Kamiokande NSI limits Additional results under the L e − L τ symmetry Results assuming inverted neutrino mass ordering 30

Figure 1 :
Figure 1: Estimated present-day upper limits on the coupling strength, g ′ αβ , of the new boson, Z ′αβ , with mass m ′ αβ , that mediates flavor-dependent long-range neutrino interactions.Our results are based on estimates of the measurement of flavor composition of high-energy astrophysical neutrinos in IceCube using 8 years of data[15] (through-going muons plus HESE, see Section 3) and on present-day uncertainties in the neutrino mixing parameters[16, 17], assuming normal neutrino mass ordering.Limits are on neutrino-electron interactions, under the L e −L µ and L e −L τ symmetries, and on neutrino-neutron interactions, under the L µ −L τ symmetry.For the latter, we assume a Z ′ µτ -Z mixing strength of (ξ − sin θ W χ) = 5 × 10 −24[18].Existing limits are from a recent global oscillation fit[19], atmospheric neutrinos[20], solar and reactor neutrinos[21], and non-standard interactions (NSI)[22][23][24].For comparison, we show the proofof-principle sensitivity (1σ) based on 2015 IceCube flavor-composition measurements[25] from Ref.[26].Indirect limits[27] are from black-hole superradiance (90% C.L.)[28] and the weak gravity conjecture[29], assuming a lightest neutrino mass of 0.01 eV.Figures8 and D3show projections for the year 2040; Fig.E4shows results assuming inverted mass ordering.See Section 4.2 for details.Our results show that limits obtained from present-day flavor-composition measurements may improve significantly on existing ones.

Figure 2 :
Figure 2: Feynman diagrams for the contributions to the neutrino-matter interaction Lagrangian, Eq. (2.1).(a) Neutrino-matter interaction in the Standard Model, mediated by the neutral gauge boson Z.In our analysis, this is important only for neutrinos inside the Earth, where matter densities are high.(b) Neutrino-matter interaction in the U (1) Lα−L β model, mediated by the new Z ′ αβ gauge boson.(c) Neutrino-matter interaction via Z ′ αβ -Z mixing.In our analysis, for the L e − L µ and L e − L τ symmetries, only diagram (a) (and only inside the Earth) and diagram (b), sourced by electrons, contribute, and we do not consider diagram (c).For the L µ − L τ symmetry, only diagram (a) (and only inside the Earth) and diagram (c), sourced by neutrons, contribute.See Section 2.1 for details.

Figure 3 :
Figure 3: Long-range matter potential, V αβ , experienced by a neutrino at Earth, Eq. (2.7).The potential is sourced by electrons and neutrons in the Earth, Moon, Sun, Milky Way, and the distant Universe as a function of the mass and adimensional coupling of the new U (1) gauge symmetry.Left: For the L e − L β symmetries (β = µ, τ ), sourced by electrons.Right: For the L µ − L τ symmetry, sourced by neutrons.Isocontours are drawn at the illustrative values of V eβ = V µτ = 10 −18 eV.(For the L µ − L τ symmetry, values of the coupling run lower due to the mixing between the new mediator and the SM Z boson.)The step-like transitions in the potential signal the onset of contributions from sources of electrons and neutrons located at different distances relative to the interaction range.See Section 2.2 for details.

Figure 4 :
Figure 4: Average flavor-transition probabilities, Eq. (2.19), as functions of the new matter potential induced by the U (1) gauge symmetries L e −L µ (left column) and L µ −L τ (right column).Probabilities are for neutrinos and antineutrinos, at a fixed energy of 100 TeV, assuming normal mass ordering (NMO).The features in the probabilities are due to resonances induced by the new potential.Vertical lines mark the values of the potential for which θ m 12 and θ m 13 become resonant.Probabilities under standard oscillations ("Std.osc."), i.e., for V αβ = 0, are shown for comparison.The mixing parameters are fixed at their present-day best-fit values from NuFit 5.1 [16, 17].Appendix D contains results for L e −L τ ; they are similar to L e −L µ .Appendix E contains results assuming inverted mass ordering.See Section 2.3 for details.

Figure 5 :
Figure5: Flavor composition of high-energy astrophysical neutrinos at Earth, f α,⊕ , as a function of the long-range matter potential V eµ , under L e − L µ (left), or V µτ , under L µ − L τ (right).We show results for three benchmark choices of the flavor composition at the sources, (f e,S , f µ,S , f τ,S ) = (1/3, 2/3, 0), the canonical expectation, (0, 1, 0), and (1, 0, 0), compared to the expectation from standard oscillations ("Std.osc.").The flavor composition is averaged between neutrinos and antineutrinos, assuming they exist in equal proportion in the flux.For this plot only, we fix the neutrino energy to 100 TeV, assume normal neutrino mass ordering, and vary the standard mixing parameters within their 1σ allowed ranges; our statistical methods systematically vary these choices.Flavor-composition predictions are compared against three estimates of the flavor sensitivity of neutrino telescopes using 8 years of IceCube data ("IC 8 yr"), 15 years of IceCube data plus 10 years of IceCube-Gen2 ("IC 15 yr + Gen2 10 yr"), and the combined exposure of all neutrino telescopes available by 2040 ("Combined ν telescopes").For L e − L τ , results are similar as for L e − L µ ; see Appendix D. See Section 3.3 for details.

LFigure 6 :
Figure6: Posterior probability density of the long-range matter potentials V eµ , under the L e − L µ symmetry (left) and V µτ , under the L µ − L τ symmetry (right).We assume a measurement of the flavor composition at Earth centered around the canonical expectation of about (1/3 : 1/3 : 1/3) ⊕ , corresponding to f S = (1/3 : 2/3 : 0) at the sources; see Fig.5.See Fig.7for the limits on the potential derived from the posterior.Results for V eτ , under the L e − L τ symmetry, are similar to V eµ ; see Appendix D. We assume that all neutrino telescopes available by 2040 have detection efficiencies similar to that of IceCube (Section 3.4).Results here are obtained assuming normal neutrino mass ordering (NMO); results for inverted ordering are in Appendix E. See Section 4.1 for details.

Figure 7 :
Figure 7: Upper (95% C.L.) on the long-range potential induced by the L e − L µ , L e − L τ , and L µ − L τ symmetries, estimated for 2020 and projected for 2040.For comparison, we show the standard-oscillation contribution ("Std.osc.") to the Hamiltonian, ∆m2  21 /(2E), evaluated at the present-day best-fit value of ∆m 2 21[16, 17]  and at different values of the neutrino energy, E. We assume that all neutrino telescopes available by 2040 have detection efficiencies similar to that of IceCube (Section 3.4).Figure8shows these limits translated into limits on the coupling, g ′ αβ .See Appendix F for the values of the limits.See Section 4 for details.

Figure 8 :
Figure8: Estimated present-day (2020) and projected (2040) upper limits (95% C.L.) on the coupling strength, g ′ αβ , of the new boson, Z ′ αβ , with mass m ′ αβ , that mediates flavor-dependent long-range neutrino interactions.Same as Fig.1, but now showing also the projected 2040 limits using either IceCube plus IceCube-Gen2, or those combined with Baikal-GVD, KM3NeT, P-ONE, and TAMBO.Left: Limits on neutrino-electron interactions under the L e − L µ gauge symmetry.Right: Limits on neutrino-neutron interactions under the L µ − L τ gauge symmetry; existing limits from accelerator neutrinos in MINOS (95% C.L.) are from Ref.[18].See Section 4 for details.

Figure D3 :
Figure D3: Constraints on the long-range matter potential V eτ , under the L e − L τ symmetry.Left: Posterior probability density of V eτ .Same as Fig. 6, but for V eτ .Right: Estimated (2020) and projected (2040) upper limits (95% C.L.) on the coupling strength, g ′ eτ , of the new boson, Z ′ eτ , with mass m ′ eτ .Same as Fig. 8, but for V eτ .See Section 4.2 and Appendix D for details.

Figure 7
in the main text shows a graphical representation of them and Figs. 1, 8, D3, and E5 shows them translated into limits on the coupling, g ′ αβ .

Figure E3 :LFigure E4 :
Figure E3: Flavor composition of high-energy astrophysical neutrinos at Earth, f α,⊕ , as a function of the long-range matter potential V eµ , under L e − L µ (top left), V eτ , under L e − L τ (top right), and V µτ , under L µ − L τ (bottom).Same as Figs. 5 and D2, but for inverted mass ordering (IMO).See Section 3.3 for details.