Magnetar-powered Neutrinos and Magnetic Moment Signatures at IceCube

: The IceCube collaboration pioneered the detection of O ( PeV ) neutrino events and the identification of astrophysical sources of high-energy neutrinos. In this study, we explore scenarios in which high-energy neutrinos are produced in the vicinity of astrophysical objects with strong magnetic field, such as magnetars. While propagating through such magnetic field, neutrinos experience spin precession induced by their magnetic moments, and this impacts their helicity and flavor composition at Earth. Considering both flavor composition of high-energy neutrinos and Glashow resonance events we find that detectable signatures may arise at neutrino telescopes, such as IceCube, for presently unconstrained neutrino magnetic moments in the range between O (10 − 15 ) µ B and O (10 − 12 ) µ B .

In the context of high-energy neutrinos, previous studies of νMM were chiefly focused on neutrino propagation through the intergalactic medium [22][23][24][25].There, the idea is that despite the small magnetic field involved, neutrinos carrying large νMM can undergo efficient flavor conversion across large propagation distances.While the results of these analyses imply potentially observable effects, it should still be noted that the intergalactic magnetic fields presently come with rather large uncertainties.An alternative to studying the appearance of νMM effects over the cosmic distances is to investigate potential impact across much smaller scales, namely in the vicinity of neutrino production.Clearly, this would require involvement of magnetic fields that are orders of magnitude larger compared to the intergalactic ones.There are several types of astrophysical objects that source such magnetic fields [26,27] and, if they also produce neutrinos, νMM effects may manifest.This was sketched in [28] for Dirac neutrinos passing through white dwarfs.
In fact, astrophysical environment with large magnetic field is a prerequisite for the acceleration of charged particles which source high-energy neutrinos [29][30][31].In this work, we assume Majorana nature of neutrinos and adopt magnetars [32] as a case study.We are motivated by the fact that (i) magnetars source extremely strong magnetic fields and (ii) it has been demonstrated in the literature that such objects produce high-energy neutrinos [33][34][35][36][37]. Therefore, considering high-energy neutrinos with large νMM in such astrophysical environments appears particularly motivating.We focus on two key phenomenological features at neutrino telescopes -the flavor composition of high-energy neutrinos with energy E ≳ 100 TeV and the Glashow resonance events.IceCube already pioneered in measuring both [38][39][40].
The paper is organized as follows.In Section 2, we setup the framework for the flavor evolution of neutrinos carrying νMM in the magnetic field.In Section 3 we focus on the magnetic field and its roles in neutrino production as well as estimate the impact of νMM in realistic astrophysical scenarios of interest.In Section 4, focusing on magnetar systems, we start by specifying representative magnetic field profile used in the analysis.Subsequently, we discuss details of our numerical simulations and show respective results.In Section 5 we utilize these results in order to assess νMM prospects at neutrino telescopes, namely its impact both to the flavor composition of high-energy neutrinos and the number of Glashow resonance events.We summarize in Section 6.Unless stated otherwise, all quantities are given in units G = ℏ = c = 1.Indices a, b denote spacetime components, α, β denote neutrino flavors and j, k denote neutrino mass eigenstates.

The Formalism
We base this study on the following effective interaction term where µ αβ is the component of the νMM matrix µ, ν L is the left-handed neutrino field, and F ab is the electromagnetic field strength tensor.In terms of density matrix ρ(t, p) = |ν(t, p)⟩⟨ν(t, p)|, the dynamics of the neutrino state can be described by where the inelastic collision term C parametrizes neutrino flux attenuation, and the Hamiltonian H reads Here, H vac , H MSW and H νMM quantify neutrino oscillations in vacuum, the Mikheyev-Smirnov-Wolfenstein (MSW) matter terms [41][42][43], and the νMM effects arising from Eq. (2.1), respectively.On top of the flavor mixing induced by vacuum oscillations and MSW effects, additional helicity1 and flavor mixing will be induced by nonzero µ αβ when neutrinos propagate through the background magnetic field.Hence, the basis state |ν(t, p)⟩ encodes the irreducible number of neutrino flavors and helicities.Notice that in building Eq. ( 2.1) we used only left-handed neutrino fields which already indicates that in this study we will be chiefly focused on Majorana neutrinos.All the expressions below are valid for Majorana neutrinos; the Dirac neutrino case follows straightforwardly and, for completeness, we outline it in Appendix A. For the Majorana case, we have three flavors with {α, β} = {e, µ, τ } and two helicities denoted as {h, h ′ } = {1 (for ν L ), 2 (for νL )}.The basis state |ν(t, p)⟩ is therefore a dimension-6 vector and H can be expressed as a 6 × 6 matrix with the entries given by The left term in Eq. 2.4 is responsible for standard neutrino oscillations.Here, m j is the neutrino mass corresponding to the mass eigenstate j, U is the leptonic mixing matrix, and V is the MSW matter potential (minus sign for antineutrinos) which is diagonal in the flavor basis.The right term in Eq. 2.4 is the helicity-flipping (HF) term induced by νMM, where B ⊥ is the strength of the magnetic field projected to the plane perpendicular to the propagation direction of neutrinos, and the angle ϕ parameterizes the orientation of B ⊥ within that plane.We note that µ αβ is asymmetric due to CPT symmetry (diagonal terms vanish) indicating that a neutrino with a flavor α will be converted to an antineutrino with a different flavor β ̸ = α.The neutrino flavor evolution can only be performed numerically due to the complexity of H. Nevertheless, different components of H may dominate at different characteristic length scales which can simplify the overall treatment, see discussions below.

Characteristic Length Scales
The length scales associated to the terms in Eq. (2.3) are given by L jk = 2π/|E j − E k |, where E j and E k denote the eigenvalues of the corresponding Hamiltonian component.The minimal length scale (L ch ) is then determined by max{E j −E k }.A crucial feature of having a HF term is that the dimension of H νMM is such that more combinations of E j and E k are possible in comparison to the standard scenario with H vac and H MSW terms; hence, multiple length scales may arise.In the absence of νMM effect, the dimension would be reduced to standard three flavors.For neutrinos traveling as separated mass eigenstates, the dimension may reduce for Dirac neutrinos if the diagonal entries µ αα dominate.In such cases, the induced νMM effects would not be subject to mass-splitting suppression.However, for the Majorana case, since non-vanishing µ αβ terms mix different flavors, such reduction cannot be made.Majorana neutrinos also suffer a large suppression of νMM effects from the mass-splitting terms when traveling over cosmic distances through regions with a weak intergalactic magnetic field of O(µG) [25].
From the left term in Eq. (2.4) we infer that the minimal length scale for vacuum oscillations is set for the largest m 2 j − m 2 k (atmospheric mass squared difference) As far as the MSW term is concerned, the refractive length is given by 2π where n e (n n ) denotes number density of electrons (neutrons).Due to the opposite signs for neutrinos and antineutrinos in Eq. (2.4), V α ∈ {V νe , V νµ,ντ , V νe ≡ −V νe , V νµ,ντ ≡ −V νµ,ντ }.The minimal length scale associated to the matter term reads where ρ B is the density of baryon matter, and Y e is the fraction of electrons.Note the relation n e = Y e n b , where n b = n p + n n is the baryon number density and n p = n e is the proton number density.Here, the benchmark value taken for ρ B aligns with numerical simulations for magnetar born in a merger event [44].Furthermore, as elaborated in the following, larger values of ρ B would lead to significant attenuation of the neutrino flux.
For neutrinos carrying O(PeV) energy, the cross sections for neutrino scattering off nuclei, A, are of the order 10 −33 cm 2 [45,46] which is comparable to G F ; and the cross section for neutrino-electron scattering is smaller.In passing, let us stress that νMM contribution to neutrino scattering is negligible since it is proportional to µ 2 αβ ≪ G F .Therefore, the length scale of the MSW effect compared with the collision length, given by L C = 1/(σ νA n A ), can be estimated as where in making this comparison we dropped O(1) factors depending on particle number densities entering L C and L M MSW .One can thus infer that, for E ∼ O(PeV), MSW effect becomes relevant at the length scale for which neutrino absorption also becomes significant.This is the rationale behind not considering MSW effects in this work, as we are not interested in the case where total neutrino fluxes get strongly attenuated.For instance, if neutrinos are produced 10 5 km above the magnetar, the flux would be attenuated when ρ B > 0.01 g cm −3 at the production site.The inverse characteristic length scales of the νMM effect (L M, −1 νMM , red), MSW effect (L −1 MSW , blue), vacuum oscillation for 1 PeV neutrinos (L −1 vac , yellow), and the inverse length itself (grey).The regions labeled as I, II and III correspond to the νMMdominated, "quiet" and vacuum-dominated regions, respectively.In both plots, the magnetic field strength and the matter density are localized within R = 10 5 km.For all lines, the νMM length scale is fixed, such that dℓ L M, −1 νMM = 1.The length scale for MSW effect is taken for different values of ρ B , such that the transparency (1 − dℓ L −1 C ) takes the respective labeled values in the left plot, and 90% in the right plot.
Obtaining eigenvalues of H νMM for Majorona neutrinos is non-trivial due to the asymmetricity of the νMM matrix.Nonetheless, we found {E j , E k } ∈ {0, a ν B ⊥ }, where degeneracy is implied and a 2 ν = µ 2 eµ + µ 2 eτ + µ 2 µτ , see Section 2.2.Hence, the minimal length scale can be estimated as (2.8) Given that L vac is dependent on neutrino energy but uniform over space, in contrast to L M νMM that depends on the magnetic field strength, different effects may dominate in different regions.In this work, our goal is not to perform an exhaustive analysis of realistic magnetic field and matter density profiles that are extremely difficult to model.Instead, we discuss a scenario which captures the most important features of the relevant physical systems.In particular, we consider a dipolar magnetic field (see Section 4) for r ≤ R, and B ⊥ ∝ (r/R) −3 for r > R. Here, r is the distance from the magnetar as the neutrinos travel, and R corresponds to the production site.
For systems with R ≪ L vac , we find two typical distances, L 1 and L 2 , such that L M νMM ≪ L vac for l ≤ L 1 , and L M νMM ≫ L vac for ℓ > L 2 .The existence of L 1 and L 2 indicates that νMM and vacuum effects can be decoupled.Fig. 1 shows the probability of a helicity flip from ν µ → νe , and a comparison of characteristic length scales, for a simple example where PeV neutrinos are produced at r = R and then travel outward across distance ℓ = r − R. Here, the matter density profile is taken as ρ B ∝ (r/R) −2 , Y e = 1 (conservatively for the matter effect), and σ νA = 6×10 −33 cm 2 .The propagation is split by L 1 and L 2 into three regions with νMM and vacuum effects dominating in region I and III, respectively.Note that in region II, although the three effects may be comparable with one another, none of them is large enough to induce a considerable flavor/helicity transition.In other words, since all length scales exceed the size of the region (min{L ch } > L 2 −L 1 within region II), none of them can amount to any significancy.We refer to this region as a "quiet" one.By comparing results with various ρ B (r = R) in the left plot of Fig. 1, one can also infer how the MSW effect would affect the flavor evolution.In the figure, labels of 100%, 90%, 50%, 10%, correspond to ρ B (r = R) = 0, 0.02, 0.1, 0.9 g/cm −3 , respectively.For all cases, the evolution of neutrino transition probability over the region II is not inducing any significant effect.
Owing to the decoupling between νMM and vacuum effects, the neutrino transition probability, from production to detection at distance L, can thus be approximately written as (2.9) Furthermore, when L 2 − L 1 ≪ L vac , this relation effectively simplifies to Here, L 1 < L cut < L 2 and we found that results only mildly depend on the choice of L cut .

νMM effect
Let us calculate the first term in Eq. (2.10) for which νMM effect dominates.The density matrix, introduced in Eq. (2.2), will evolve from initial time t 0 by an infinitesimal time δt as ρ(δt where the exponent of the time evolution operator reads Focusing on Majorona neutrinos we define the normalized νMM tensor for which it is convenient to introduce a unitary matrix Ũ (not to be confused with the leptonic mixing matrix U aforementioned) such that Ũ |μ| 2 Ũ † = diag(0, 1, 1).After some straightforward derivations, we find For the scenario where ϕ is adiabatically evolving, namely when we can write down the density matrix at time t as ρ(t, p) = e −iA(t) ρ(0, p) e iA(t) .Furthermore, within the νMM length scale (t < L νMM ), the time evolution operator reads with (2.17) The transition probability from ν h α to ν h ′ β can then be calculated using Given the above, for the helicity-conserving (HC) case with h = h ′ we find while for HF with h ̸ = h ′ the transition probability reads Here, 0 ≤ Φ jk ≤ 1 and Φ jj = 1, are the decoherence terms describing the overlap between j and k eigenstates in the phase space [47].In contrast to the case of vacuum oscillations, the time evolution in νMM scenario cannot be diagonalized with real eigenvalues, and is instead described through the rotation matrices Ũ .This implies that there will always be oscillation effect from the Φ jj term even if Φ jk = 0 for k ̸ = j.Furthermore, the fact that νMM of Majorona neutrinos possesses an additional asymmetricity further mitigates decoherence effects.In particular, since one of the νMM eigenvalues is zero and the other two are degenerate, only the terms corresponding to the same mass eigenstates (Φ jk = 1) will survive after inserting Eq. (2.17) into Eq.(2.20).Hence, unlike for the vacuum oscillations where decoherence arises due to the separation of wave packets for different mass eigenstates (Φ jk < 1), there will be no such decoherence effects entering in the Majorana neutrino HF transition probability.
Let us now investigate HC and HF transition probabilities for θ ν ≪ 1.In that case, the two-step transitions, e.g.ν e → νµ → ν τ , can be neglected so the flavor transitions for the HC case should vanish, unlike in the HF case.This can be checked by starting from Eqs. (2.19) and (2.20) and plugging in cos θj ≃ 1 and sin θj ≃ θj for all j.We find (2.21) The optimal HF scenario (θ ν = π/2), is shown in the left panels of Fig. 2 for various flavor structures of νMM.Similarly, in the right panels we show the case where equilibrium is reached (θ ν is averaged out).For both cases, the upper panels show flavor-conserving probability P αα = P HC αα +P HF αα and the lower ones corresponds to flavor-changing probability P αβ = P HC αβ + P HF αβ (α ̸ = β).With only one non-vanishing term μβγ , the transition occurs among flavors β and γ leaving the ν α state unaffected; this is evident from the top corner in the left panels.Namely, P αα = 1 while P αβ = P αγ = 0 for α ̸ = β ̸ = γ.In the lower left corner of panels (a) and (c), μαβ = 1 yields a 100% transition from ν α to ν β in the optimal θ ν = π/2 scenario.In contrast, in the scenario where θ ν is averaged out, ν α and ν β would each have an equal share; see the lower left corner in the right panels.In fact, such an equilibrium feature holds in general; for instance, the center of panels (b) and (d) indicates that all flavors have a 1/3 share when μeµ = μeτ = μµτ .In general, it can be shown that the oscillation probabilities are unitary (2.22) as also evident from panels (a) and (c), and panels (b) and (d), respectively.

Roles of the magnetic field
In astrophysical environments, the magnetic field makes an impact to the energy and flavor composition of produced neutrinos in a rather complicated manner.In particular, possible progenitors of high-energy neutrinos are relativistic protons and the magnetic field plays the dominant role in their acceleration (Section 3.1).Therefore, the energy budget inherited by neutrinos will be limited by the magnetic field strength.The amount of energy that will be transferred to neutrinos depends also on the sequential production processes: proton collisions first generate mesons which then decay to neutrinos.During these steps, the magnetic field works to dissipate the energy from charged particles mainly via synchrotron radiation (Section 3.2).While the aforementioned aspects are related to the energy of neutrinos, note that the strength of the νMM effect is also closely associated to the magnitude and spread of the magnetic field, as was elaborated in Section 2.2.In this section, we discuss multifaceted roles of magnetic field, and based on the parameters on the so called Hillas plot [48,49], which shows the size of the accelerator region and magnetic field strength for various astrophysical environments, we will estimate (Section 3.3) realistic strength of νMM effects.

Acceleration of protons
Relativistic protons have been long thought as triggers of processes producing high-energy neutrinos carrying O(PeV) energy through, e.g.proton-photon (pγ) and/or proton-proton (pp) collisions with the subsequent meson decays.The protons can be accelerated by In the left panels we assumed that the HF effect is maximized (i.e.θ ν = π/2) while the right panels correspond to the case where θ ν is averaged out.
magnetic field and can in principle acquire energy up to [29,48] E p max = ηeBRΓ ≃ 10 5 η 0.1 Here, R is the size of the acceleration region in the inertial frame of reference, B is the strength of magnetic field that proton experiences when it reaches the energy of E p max , Γ is the Lorentz factor and η is the acceleration efficiency (lower for shorter acceleration times) tied to the acceleration mechanism and the velocity of the moving source.A portion of the proton's energy can be handed over to the generated neutrino; in average, the attainable neutrino energy can be estimated as E p max /20 [30,31,50,51].

Cooling
The condition in Eq. (3.1) does not take into account various cooling processes experienced by protons and mesons.Cooling effects do not only affect the energy of neutrinos, but also .Adiabatic (black dashed) and synchrotron (solid) cooling timescales for secondary mesons, namely µ ± , π ± , and K ± .Their decay timescales are also shown (dash-dotted).These results are obtained following [31], with B = 10 5 G at R = 10 5 km.
alter the neutrino flavor composition since some decay processes may become irrelevant in producing neutrinos of very high energy.The timescale of cooling or dissipation, t dis , is generally jointly set by the two dominant energy loss mechanisms: the adiabatic loss due to the expansion of the shell (with the timescale t ad ) and the synchrotron loss of the protons (with the timescale t syn ) through Depending on the acceleration (t acc,p ) and cooling (t dis,p ) timescales of protons, and the cooling (t dis,m ) and decay (t dec,m ) timescales of mesons, there are in principle three options for the generation of high-energy neutrinos: (i) For t acc,p > t dis,p , the protons are unable to be sufficiently accelerated, thus no highenergy neutrinos are expected.
(ii) For t acc,p < t dis,p , and t dec,m > t dis,m , protons can be highly accelerated and subsequently produce energetic mesons through pp and/or pγ collisions, but mesons will not have enough time to pass over their energy to daughter particles before cooling down.
(iii) For t acc,p < t dis,p , and t dec,m < t dis,m , protons can yield energetic mesons, and the latter decay to high-energy neutrinos.Channels that fall in this category may determine the initial flavor composition.For magnetar-powered gamma-ray burst (GRB), this usually happens 10 4 -10 5 km away from the magnetar, where the magnetic field has reduced from O (10 16 ) G around the star's surface to O(10 5 ) G [52].PeV energy.The solid lines stem from Eq. (3.4) for the region where Eq. (2.10) is applicable.The yellow contour lines show the maximum proton energies which are able to generate neutrino energy above 100 TeV, considering cooling effects and an acceleration efficiency of η = 0.1.The brown region features dipolar magnetic field [54,55] and the blue region represents GRB case [56,57].The red triangle is the considered benchmark.
For (ii) and (iii), the comparison between meson decay and cooling timescales is shown in Fig. 3 (c.f.[31,53]).We see that muons with energy larger than O(10) TeV will be damped by synchrotron radiation before they decay, while pions with energy up to ∼ 1 PeV can decay before losing significant energy via cooling.Due to their higher decay rate and larger mass, kaons with energy up to 50 PeV can decay before experiencing significant energy loss.Therefore, since the cooling process is non-trivial and energy dependent, we will discuss all typical possibilities in Section 5.

Estimated magnitude of νMM effect
Based on Eq. (2.17), if neutrinos experience a field proportional to ℓ −α , where α > 1, then where L cut ≫ R fulfills the criteria given in Eq. (2.10).This relation is determined through B and R and thus can be related to the maximal proton energy in Eq. (3.1).Using Eq. (3.3) we can estimate a ν for which we expect appreciable νMM effects.To this end, we set θ ν = 1.
Then, for the potential reach of a ν , denoted by a ν,Hillas , we find where we have defined ξ = α−1.If the condition for decoupling of νMM and vacuum effects is no longer met, the interplay with vacuum oscillation yields a more complicated form than Eq.(3.4).In this case, we define a ν,Hillas as the νMM which could give P HF = sin 2 1 at some point.In Fig. 4, the dashed and dotted lines at R > 10 8 km show such scenario, calculated via numerical evolution.The lines saturate to the value where L νMM ∼ L vac since, otherwise, vacuum oscillation would suppress P HF too much for νMM to induce an O(1) effect.Furthermore, the fact that the dashed lines overlap with the solid lines at R ≪ 10 9 km, confirms that our analytic expression from Section 2.2 agrees with the numerical result, as long as the decoupling condition is satisfied.The benchmark values of B = 10 5 G and R = 10 5 km yield a ν,Hillas of O(10 −12 )µ B and such values of νMM are presently unconstrained.These values of B and R are also indicated in Fig. 4 (red triangle).
The yellow contour lines are plotted by a joint consideration of the Hillas acceleration budget (3.1) and the proton cooling effect [case (i) in Section 3.2].In particular, the cooling from synchrotron radiation (adiabatic expansion) dominates at larger (smaller) B, resulting in a horizontal behavior (upward shift of the oblique line) [53,58].In Fig. 4 we also show regions in R and B populated by some known astrophysical objects and those that fit well with the parameter space discussed above.
In what follows, we will choose a specific astrophysical environment -slowly rotating magnetar -and perform a detailed numerical analysis, going beyond analytical estimates in Eq. (3.4).Magnetars yield extremely high magnetic field strength (10 15−17 G) near the surface of the star and in order to avoid significant cooling effects we will consider neutrinos generated 10 4 − 10 6 km away from the star.We will also scrutinize a realization in which such neutrinos travel towards the star and enter regions with rather large magnetic field, potentially inducing a strong νMM effect.

Magnetic Field Structure
We approximate the neutron star spacetime by the following line element in the rest-frame of the star ds 2 = −e −2Ψ dt 2 + e 2λ dr 2 + r 2 dθ 2 + r 2 sin 2 θdφ 2 . (4.1) Here, (t, r, θ, φ) represent the Schwarzschild coordinates, and Ψ is the lapse function of r and λ that is connected to the magnetar mass, M ⋆ , via Note that, in Eq. (4.1), we have neglected the influence of magnetic field in the stellar shape and the influence of the star's spin; this allows us to perform a semi-analytical calculation.
For simplicity, we employ dipolar magnetic field given that such configuration has been used to model the fueling of GRB afterglow for a post-merger magnetar (see e.g.[59]).In this work, we also ignore the azimuthal (toroidal) component of the magnetic field since such structure is highly uncertain and its treatment is rather involved.By denoting the unit vector normal to a space-like hypersphere as n and by setting the magnetic axis along z-direction as shown in Fig. 5, the dipolar magnetic field sourced by a neutron star can be expressed as [60,61] Here, B ⋆ is the magnetic strength at the magnetar's pole and vanishing toroidal component of the magnetic field is achieved by setting ζ(ψ 2 ) = 0. We take B ⋆ = 3 × 10 16 G which could yield B = 10 5 G at R = 10 5 km.The stream function ψ 2 for dipolar field has the form [54,62] with f (r) being a function related to certain boundary conditions, and Y 20 being the spherical harmonic function of degree 2 and order 0. The exterior part is characterized by with R ⋆ being the radius of the remnant neutron star which is set as 15 km here.This form admits that the magnetic field is force-free outside of the star and has zero-current on the stellar surface.Although the expression represents a static field, we emphasize that we are considering the slow-rotation limit of the magnetic configuration, and thus this expression is valid only within the light cylinder of neutron stars.

Simulation framework
As illustrated in Fig. 5, after being produced at a given location on the shell with radius R, neutrinos propagate towards the Earth along ⃗ ℓ at an angle θ E relative to the axis of the magnetic field.Although there are other directions as well, we are only interested in the portion of the flux that goes towards us, i.e. neutrinos streaming in a direction perpendicular to the yellow plane determined by θ E .Neutrinos produced on the upper shell (light blue) in Fig. 5 will be in-going and those produced on the lower shell (dark blue) will be outgoing with respect to the star.Furthermore, the production site can be projected along the direction of ⃗ ℓ onto the yellow plane.The projected point on the plane is determined by the distance R i = R sin θ R to the stellar center and θ B , see Fig. 5.To summarize, for a certain magnetic field structure, the geometric parameters that enter into the computation of θ ν in Eq. (2.17) are θ E , θ B , θ R and R. The two factors determining B ⊥ (ℓ) are the distance R i and the structure of the magnetic field at R i .The former is determined by R and θ R and the latter depends on θ B and θ E .
In our Monte Carlo simulation, these four parametersθ E , θ B , θ R and R -are sampled for the in-going and out-going neutrinos, depending on the system of interest.Let us take a single source with neutrinos being produced by spherical shock wave collision.In this scenario, neutrinos are expected to be produced uniformly on the shell and emitted isotropically (see Refs. [63,64]).This implies that at fixed θ E and R, θ B and θ R are uniformly sampled between [0, 2π] and [0, π], respectively.Here, θ R ∈ [0, π/2] is for in-going neutrinos and θ R ∈ [π/2, π] for out-going neutrinos.The sampling of R and θ E is more involved since it depends on specific magnetar systems, namely those listed in Appendix B. Nevertherless, owing to the symmetry of the magnetic field with respect to the pole and its equatorial plane, we can sample θ E in the range [0, π/2].To capture essential features, we will consider the following four benchmark cases with two extreme ones for the ratio of in-going and out-going neutrinos: (i) in-going : out-going = 1 : 1 (0 < θ R < π) and (ii) in-going : out-going = 0 : 1 (θ R > π/2).As for R, we will simulate two cases: one with a fixed value R = 10 5 km (the red triangle in Fig. 4) and another one where this parameter is uniformly sampled within log 10 (R/km) = [4,6].
In Section 2.2 we have shown that the evolution operation for density matrix has a simple form, see Eq. (2.16), provided the adiabatic condition in Eq. (2.15) is satisfied.Taking a ν ∼ a ν,Hillas and R = 10 5 km, we numerically test the level of adiabaticity by having 10 × 10 × 10 samples of θ E , θ B and θ R , each over 1000 steps of ℓ ≤ 10 6 km.Among these samples and steps, we find that only O(1%) of them are not adiabatic ( φ/ϕ ≥ 0.1a ν B ⊥ ).

Numerical Results
In this section, we explicitly obtain the probability distributions of θ ν using standard Monte Carlo methods for the aforementioned four benchmark cases.We first fix R = 10 5 km and obtain predictions for dℓB ⊥ (ℓ; θ E , θ B )/(BR) as a function of θ E and θ B , for θ R = π/2 (where the in-going and out-going scenarios overlap).The result is shown in the left panel of Fig. 6, and we obtain the expected value of ξ in Eq. (3.4) to be 2.66.The probability distribution of θ ν for out-going neutrinos is obtained by scanning over θ R ∈ [π/2, π], as shown in the right panel of Fig. 6 for a fixed a ν = a ν,Hillas = 9.1 × 10 −12 µ B .On the other hand, for the in-going case, we have identified the dependence on by varying θ R within the range [0, π/2].This R −2 i dependence can be understood since magnetic field falls as R −3 i , and one power of R i is compensated by dℓ.The probability distribution for the in-going neutrino case is shown in Fig. 7, where we considered three cases with fixed R (the gray curves) and one case with R sampled in the range between between 10 4 and 10 6 km (red).
Armed with θ ν probability distributions for both in-going and out-going case, we combine them by weighting them according to the considered flux ratio.This gives us the prediction for θ ν , with which we can obtain HF probability, P HF .While in Figs. 6 and 7 we employed a ν = a ν,Hillas = 9.1 × 10 −12 µ B , in general one should consider a spectrum of unconstrained values for a ν .The result of such analysis is shown in Fig. 8 where we show HF transition probability as a function of a ν .The figure shows the averaged probability; at small a ν , the contributions with largest dℓB ⊥ are most relevant (the portion to the left of the peak in Figs. 6 and 7).In this regime, HF probability scales with a 2 ν , following the behavior in Eq. (2.21) from small θ ν expansion.For larger values of a ν we observe that the oscillation features develop.The prerequisite for such behavior is that probability distribution peaks around a particular value of θ ν ≳ 1.Finally, when a ν is large enough, the probability averages to 1/2.Note that, although Eq. (2.16) may no longer apply in this scenario, the transition probabilities will still average to 1/2.From Fig. 8 we can generally infer that appreciable HF transition probability occurs at a ν values as low as approximately a ν,Hillas and 10 −2 × a ν,Hillas for fixing R = 10 5 km and varying R, respectively, when the total relevant flux consists of out-going neutrinos.The reach in a ν improves by 2 − 4 orders of magnitude for the optimistic case where ingoing and out-going neutrinos have the same contribution in the relevant flux.Additionally, to accommodate variations of different astrophysical objects, we also show, as red line in Fig. 8, the case where we sample over the value of B × R ∈ [10 9 , 10 11 ] G km which is within the GRB region of Fig. 4. In summary, given that a ν,Hillas = 9.1 × 10 −12 µ B is our reference value, we observe that appreciable P HF can be obtained for νMM values as low as O(10 −13 ) µ B for the 0 : 1 case and O(10 −15 ) µ B for the 1 : 1 scenario.Such values are presently unconstrained and in what follows we will discuss potential signatures that they may induce at neutrino telescopes.

Signatures at Neutrino Telescopes
The IceCube collaboration has measured flavor composition of high-energy neutrinos [38,39] and has also detected a Glashow resonance event [40] that corresponds to the scattering of 6.3 PeV electron antineutrino off electron with the exchanged W boson being on shell.These two observations provide us with non-trivial information about the flavor composition of high-energy neutrinos at production sites [65][66][67].In particular, five types of flavor Figure 8. Flavor transition HF probability (applies for transitions between any two flavors).We show cases where relevant fluxes consist only of out-going neutrinos and also scenarios in which the ratio between in-going and out-going neutrinos is 1 : 1.
Since IceCube cannot distinguish between neutrinos and antineutrinos at high energies, we have instead three flavor compositions considering (ν e + νe : ν µ + νµ : ν τ + ντ ), see characteristic ternary diagram in Fig. 9.The labels (empty triangles) in the plot represent the expected flavor composition in the standard case without the νMM effect and are obtained by taking present best fit values of neutrino mixing angles [68].We show 68% and 95% CL limits from IceCube (black solid and dashed) and also present expected projections for measurement of flavor ratios with IceCube-Gen2 [69] (lighter gray one corresponds to (0 : 1 : 0) and darker brown one represents (1 : 2 : 0) at production).In what follows, we assume that future observations of high-energy neutrinos will feature plenty of events (such that robust flavor composition analysis can be performed for such subset) produced in regions with high magnetic field.Among such subset, an identification of sources would also be required for the analysis.
With the P HF from Section 5, we propagate neutrinos from the production site to Earth and, for the three considered flavor compositions, we obtain predictions at IceCube -see Fig. 9.The obtained region corresponding to (1 : 2 : 0) (red) occupies, as expected, rather narrow portion of the flavor triangle.Results are more promising for (0 : 1 : 0) case; we notice that the respective region (blue) spans over significant portion of the triangle.In particular, it populates region in which flavor ratio measurements for pion decay (1 : 2 : 0) sources are expected.That could cause a degeneracy between new physics and astrophysical parameters.Nevertheless, as the blue region extends even beyond the (1 : 2 : 0) region, discovery of the νMM signatures is in principle possible.Another scenario (yellow) that we consider is neutron decay (1 : 0 : 0).One can infer from the figure that such realization is in tension with the IceCube measurements at 68% C.L. (compare yellow empty triangle with the limit).However, we observe that if new physics in the form of νMM is at play, neutron decay production mechanism becomes more consistent in light of present data since significant portion of the yellow region falls within the limit of IceCube at 68% C.L.. We should nevertheless point out that, neutron decay production mechanism has been widely taken as less likely than the other two considered above.Now we turn our attention to the Glashow resonance.As stated above, this type of events can only be induced by electron antineutrinos and, in that sense, such event topology can lead to discriminating between neutrinos and antineutrinos.That appears especially important in light of HF that occurs due to νMM.Indeed, as we show below, the neutrinoantineutrino transitions have an intriguing interplay with certain production mechanisms.
In Fig. 10, we show the fraction of high-energy electron antineutrino flux (f νe ) as a function of a ν for several astrophysical production mechanisms.In calculating that, we employed HF transition probability which matches the one for obtaining red line in Fig. 8; specifically, only out-going neutrinos are considered.In Fig. 10, we also show bands indicating future projections at neutrino telescopes, taken from Ref. [70].Notice that, for pγ-damped scenario, f νe → 0 when a ν → 0 and this is evident from all panels in which results for different textures of νMM matrix are presented.This is simply because, as discussed above, such (0 : 1 : 0 : 0 : 0 : 0) scenario does not yield antineutrinos and hence does not lead to Glashow events.However, in the considered BSM realization, due to HF transitions, antineutrinos could get produced even for such production mechanism.This is clearly shown in Fig. 10 for a ν ̸ = 0. If, in the future, at least one Glashow event from magnetar sources will be detected, pγ-damped production mechanism in such astrophysical environments may not be excluded within the νMM framework.Since we are uncertain about the origin of the detected Glashow resonance event [40], this may be relevant already for the present data.If, on the other hand, no Glashow events will be seen from magnetar sources, it will be possible to constrain νMM, as a significant fraction of electron antineutrinos will be produced for a ν ∼ a ν,Hillas across all considered production mechanisms, seen from Fig. 10.

Summary
In conclusion, our work delved into a scenario where high-energy neutrinos are generated in regions characterized by intense magnetic fields.The presence of nonvanishing magnetic moments and the resulting interactions with the background magnetic field enable such neutrinos to undergo both flavor and helicity transitions.We have performed rigorous simulations of neutrino propagation and have identified potential avenues for testing magnetic moment effects with forthcoming high-energy neutrino data.For Majorana neutrinos, such effect can washout information of the production state, by averaging out the fractions of different helicities and flavors.Specifically, for neutrino magnetic moments of O(10 −12 ) µ B , the helicity-flipping probability could reach O(1) values for all considered cases and we found that this is testable in future analyses of high-energy neutrino flavor composition and Glashow events.In the optimal scenario, where high-energy neutrinos are produced at a certain distance away from the neutron star and propagate toward it, neutrino magnetic moments can be as small as O(10 −15 ) µ B for achieving similar effects.Given that, we have demonstrated that presently unexplored values of neutrino magnetic moment could manifest in future neutrino telescope data.Hence, this work provides a framework for experimental tests of magnetic moment-related phenomena at neutrino telescopes in the realm of high-energy neutrino interactions.

B Candidates for the magnetar system
We enumerate scenarios for high-energy neutrino production in highly magnetized astrophysical environments, focusing on magnetars as emitters of jet and/or central engines for GRBs.
(i) Young magnetars tend to bear a rapid, and differential rotation, and posses a strong magnetic field with non-trivial multipolar structure.It can therefore be envisioned that the spin axis may differ from the magnetic field's axis.The unipolar induction of a rotating, magnetized neutron star will render an electric potential, which can possibly reach a magnitude of in the vicinity of the stellar surface [71].Here, Ω denotes the magnitude of stellar spin and R ⋆ is the radius of the neutron star.The associated electromotive force then accelerates the charged particles that constitute a plasma.For the cases where ⃗ B and ⃗ Ω satisfy the condition ⃗ B • ⃗ Ω < 0, the positively charged particles, including protons, will be unleashed along the open field lines in the polar regions [72].In particular, the injection rate of protons to the stellar wind [73] may be approximated by the Goldreich-Julia rate [37].Depending on the surface temperature of the pulsar, the relativistic protons are expected to hit photons and generate high-energy neutrinos via ∆-resonance [74,75], pγ → ∆ → nπ + → nν µ µ + → nν µ e + ν e ν µ , (B.2) It should be noted that the produced pions and muons will also undergo several cooling mechanisms and will hence only be able to hand over a fraction of the energy to neutrinos [76].
(ii) Within fireballs [56] that scintillate gamma-ray bursts -both long and short onesprotons and electrons will undergo Fermi acceleration by the magnetic field.Depending on the properties of fireballs (see [77] and references therein), some protons can be accelerated sufficiently to enable pp collisions.This collision will create mesons which further decays to produce neutrino transients [78,79].These pp collisions mainly lead to pion production though other mesons (e.g.kaons) can also occur [80].The fact that pions are mostly produced, however, does not necessarily imply that they are the main source of neutrinos; the less efficient radiative cooling and the shorter lifetime of kaons arguably make kaons more important source of neutrinos for energy range above TeV (thus applies to the energy range considered in the main text).We should also stress that the pγ collisions are also at play and yield photomesons.Namely, in addition to the process shown in Eq. (B.2), kaons can be produced as well and subsequently decay to neutrinos [81] pγ → Λ 0 K + , Σ 0 K + , Σ + K 0 .(B.3) (iii) At early stages of millisecond magnetars, occurring following either a binary merger or a supernova, the relativistic wind spewed from the central remnant will be braked by its interaction with ejecta, producing shocks heating up the ponderable medium.Within these pulsar wind nebula, inelastic pp collisions are expected to effectively operate, giving rise to a copious number of mesons [36], e.g.neutral and charged pions [82].These mesons may, in principle, decay and thereby generate high-energy neutrinos; yet, no such events have been detected [83].
(iv) The pp collisions within cocoon systems, formed atop the remnant magnetar of binary merger [84] or supernovae [85], will seed middle stage mesons such as pions and kaons.Their decay into highly energetic neutrinos may, however, be stagnated considerably since the accelerated mesons will be cooled via several mechanisms (e.g., [76,81]) resulting in considerable energy loss.Owing to the longer cooling times and shorter lifetimes, there are speculations that charm contribution to the neutrino population may be more important than expected [37].
(v) On top of the aforementioned proton-involved collisions, neutrons in the relativistic outflow may activate production of metastable mesons (e.g., [86]).Such events may be detected in future analyses of the brightest GRB recorded to date (GRB 221009A) [87].

Figure 1 .
Figure 1.Left plot: Transition probability ν µ → νe obtained by numerically evolving the full Hamiltonian for 1 PeV neutrinos.Right plot: The inverse characteristic length scales of the νMM effect (L M, −1 νMM , red), MSW effect (L −1 MSW , blue), vacuum oscillation for 1 PeV neutrinos (L −1 vac , yellow), and the inverse length itself (grey).The regions labeled as I, II and III correspond to the νMMdominated, "quiet" and vacuum-dominated regions, respectively.In both plots, the magnetic field strength and the matter density are localized within R = 10 5 km.For all lines, the νMM length scale is fixed, such that dℓ L M, −1 νMM = 1.The length scale for MSW effect is taken for different values of ρ B , such that the transparency (1 − dℓ L −1 C ) takes the respective labeled values in the left plot, and 90% in the right plot.

Figure 2 .
Figure 2. The ternary diagrams show how the structure of the νMM matrix with real values (relative magnitudes of the entries shown on the axes) impacts the oscillation probabilities.The panels (a) and (b) show the survival, while the panels (c) and (d) show flavor transition probability.In the left panels we assumed that the HF effect is maximized (i.e.θ ν = π/2) while the right panels correspond to the case where θ ν is averaged out.

-Figure 4 .
Figure 4.The Hillas plot.We show lines corresponding to a ν,Hillas = 10 −11 µ B (blue) and a ν,Hillas = 10 −14 µ B (red).The dashed (dotted) lines show numeric results for α = 1(3), for neutrinos with 1 PeV energy.The solid lines stem from Eq. (3.4) for the region where Eq. (2.10) is applicable.The yellow contour lines show the maximum proton energies which are able to generate neutrino energy above 100 TeV, considering cooling effects and an acceleration efficiency of η = 0.1.The brown region features dipolar magnetic field[54,55] and the blue region represents GRB case[56,57].The red triangle is the considered benchmark.

Figure 5 .
Figure 5. Sketch showing neutrino propagation from a single source.Neutrinos are produced on the shell at the distance R above the stellar surface and this can be described by θ B and R i = R sin θ R coordinates for both in-going and out-going cases.The angle between the magnetic field axis and the Earth direction is θ E .

Figure 6 .
Figure 6.The left panel shows dℓ B ⊥ /(BR) as a function of θ E and θ B .The right panel shows the probability distribution function of θ ν for out-going neutrinos when a ν = a ν,Hillas = 9.1 × 10 −12 µ B .For both panels, we have set R = 10 5 km.

Figure 7 .
Figure 7. Normalized distribution of in-going neutrinos with a ν = a ν,Hillas = 9.1 × 10 −12 µ B .The three gray lines represent the case where R is fixed to a certain value, and the red curve is obtained by logarithmically sampling R uniformly within the region between 10 4 and 10 6 km.Solid gray lines represent the case where θ E and θ B are fixed to the values of π and π/4, respectively.The dashed lines show the case in which these parameters are scanned over.

Figure 9 .
Figure 9. Ternary diagram representing flavor composition of high-energy neutrinos.The empty triangles represent the expected flavor compositions at Earth in the SM for the given initial flavor compositions.The BSM predictions for νMM are shown in blue, red and yellow for three considered flavor compositions.In producing those, all entries of the νMM matrix were considered.Relating to Fig. 8, the gradient (from darker to lighter) for each color corresponds to HF probability in the range P HF ∈ [0, 0.2], [0.2, 0.5], [0.5, 0.7], respectively.

Figure 10 .
Figure 10.The fraction of the high-energy electron antineutrino flux at Earth for various production mechanisms where BR is uniformly sampled within the GRB region in Fig.4, and only out-going neutrinos are considered (corresponding to the red line in Fig.8).Panels (a)-(d) correspond to different structures of the νMM matrix.