Constraints on pseudo-Dirac neutrinos using high-energy neutrinos from NGC 1068

Neutrinos can be pseudo-Dirac in Nature - they can be Majorana fermions while behaving effectively as Dirac fermions. Such scenarios predict active-sterile neutrino oscillations driven by a tiny mass-squared difference $(\delta m^2)$, which is an outcome of soft lepton number violation. Oscillations due to tiny $\delta m^2$ can only take place over astrophysical baselines and hence are not accessible in terrestrial neutrino oscillation experiments. This implies that high-energy neutrinos coming from large distances can be naturally used to test this scenario. We use the recent observation of high-energy neutrinos from the active galactic nuclei NGC 1068 by the IceCube collaboration to rule out $\delta m^2$ in the region $[1.4 \times 10^{-18}, 10^{-17}]\, {\rm eV}^2$ at more than $90\%$ confidence level - one of the strongest limits to date on the values of $\delta m^2$. We also discuss possible uncertainties which can reduce the sensitivity of these results.


Introduction -
The advent of multi-messenger astronomy in the last decade has been one of the biggest success stories of particle physics [1].The observation of high-energy neutrino events, in direct correlation with the gamma-ray flares from the blazer TXS 0506 + 056, identified blazars as powerful sources of these neutrinos [2,3].Since then, additional searches have been performed to locate more of these sources.The observations of these high-energy neutrinos, which can point back to their sources, indicate the presence of some of the most powerful cosmic accelerators in the Universe.
Very recently, the IceCube collaboration released a search result for neutrino emission from a list of 110 gamma-ray sources during the period 2011-2020 [4].An excess of 79 +22 −20 neutrino events, with energies of few TeV, was reported to be coming from the direction of a nearby active galactic nucleus (AGN) -the NGC 1068 -with a global significance of 4.2σ, which was a major improvement over a previous result of 2.9σ.The neutrino spectra can be fit using a power-law, Φ νµ+νµ = Φ 0 (E ν /1 TeV) −γ , with Φ 0 being the overall flux normalization at neutrino energy E ν = 1 TeV and γ the spectral index.Using this, the collaboration reported Φ 1TeV 0 = (5.0 ± 1.5 stat.± 0.6 sys. ) × 10 −11 TeV −1 cm −2 s −1 , with γ = 3.2 ± 0.2 stat.± 0.07 sys. .The distance to the AGN is estimated to be d = 14.4 Mpc, although it has been quoted previously in the literature to lie between 10.3 ± 3 Mpc [5] and d = 16.5 Mpc [6].
The discovery of these high-energy neutrinos has opened up new avenues of testing exotic physics, which was otherwise inaccessible in terrestrial laboratories.The arrival of these neutrinos clearly indicates that the Universe is not opaque to neutrinos at this energy range.This can be used to put tight constraints on different kinds of non-standard physics in the new sector (see [7] and references therein for a detailed discussion).
The naturally long baseline offered by these neutrinos allows us to test the violation of lepton number in the Standard Model (SM).If lepton number is not considered a symmetry of the SM, and we know that it is already broken by quantum effects even within the SM, neutrinos can be Majorana fermions.The extent of lepton number violation in the SM is quantified by the Majorana mass term of neutrinos in the Lagrangian.However, if lepton number is violated softly -measurable through a tiny Majorana mass term -neutrinos can be pseudo-Dirac, where they behave effectively as Dirac neutrinos [8][9][10][11][12][13][14] 1 .In this case, the neutrino mass-eigenstates corresponding to the active and sterile states develop a tiny mass-squared difference, proportional to the extent of lepton number violation.Oscillations driven by this tiny mass-squared difference (δm 2 ) will only be relevant for long baselines, inversely proportional to δm 2 .
In this work, we utilize the high-energy neutrino events observed from NGC 1068 to probe neutrino oscillations due to a small δm 2 .The advantages are two-fold: firstly, there exists a precise measurement of the location of the source d ∼ 14.4 Mpc, RA (right ascension) and DEC (declination) = (40.667,−0.0069); secondly, the neutrino energies are well measured to be around a few TeV to tens of TeV.The combination of these allows us to precisely pin down the oscillation length, and study the sensitivity of neutrino oscillations to tiny mass-squared differences.This allows us to rule out δm 2 = 1.4 × 10 −18 eV 2 at more than 90% confidence level.We emphasize that this is the smallest value of δm 2 constrained with such a high significance.An earlier work by one of the authors used data from SN1987A to probe even tinier values of δm 2 ≲ 10 −20 eV 2 [21], however, the significance was lower due to the small number of events observed.This observation by IceCube, along with the data from SN1987A, provide two of the tightest constraints on the strength of lepton number violation in the SM.Pseudo-Dirac neutrinos -To explain the masses of the neutrinos, the SM can be expanded to include three additional sterile neutrinos with at least two of them having Majorana masses M N , which characterize the extent of lepton number violation.Apart from the Majorana mass term, the neutrino mass matrix also consists of Dirac-type masses given by M D = Y v, where Y is the Yukawa coupling and v denotes the vacuum expectation value of the Higgs field.In the limit of soft lepton number violation, |M N | ≪ |Y v|, the generic neutrino mass matrix can be diagonalized by a 6×6 block-diagonal unitary matrix U , consisting of the PMNS matrix, and another 3 × 3 unitary matrix [11].A generalization using the spectral density function was recently worked out in [28].
For a tiny |M N |, the mixing between the active and sterile states becomes maximal.In this case, a neutrino flavor state can be written as a superposition of two almost degenerate mass states, with masses m 2 k,± = m 2 k ± δm 2 k /2.Here ν ± k correspond to the mass eigenstates associated with the flavor state.This scenario leads to active-sterile oscillations driven by δm 2 k /2, in addition to oscillations due to solar (∆m 2 sol ) and atmospheric (∆m 2 atm ) mass-squared differences [29].For simplicity, we will assume that δm 2 is the same for all states, and hence drop the index k.The analysis with different δm 2 becomes considerably more complicated because of multiple parameters.In that case, we expect degeneracies to arise among different δm 2 when trying to fit the data within the (3 + 3) model.In particular, for IceCube-like experiments, which are mostly sensitive to muon neutrinos, we expect the sensitivity to be largest to δm 2  2 and δm 2 3 , whereas for experiments sensitive to solar and supernova neutrinos (which mainly detect electron neutrinos and antineutrinos), the sensitivity should be mostly to δm 2  1 .Detailed investigation on this is need in the future.
For neutrinos travelling over astrophysical baselines, flavor oscillations induced by (∆m 2 sol,atm ) average out due to rapid oscillations, causing decoherence [29].However, for small enough δm 2 , active-sterile oscillations can persist.The corresponding active neutrino survival probability can be computed as [20] Here L osc denotes the oscillation length and can be computed for a given neutrino energy E and δm 2 as, (3) However, even these oscillations driven by δm 2 can get washed out due to the separation of wave packets over long distances, leading to decoherence.This is measured through the coherence length L coh , which gives the length scale up to which flavor coherence is expected to be maintained.This depends also on the width of the neutrino wave-packet σ x , which is usually determined from the process of neutrino production [30].The coherence length is estimated as, Taking into account the active sterile oscillation probability, and the flavor averaging due to the solar and atmospheric terms, the probability of obtaining a neutrino flavor ν β , starting from ν α is given by This is the new physics that we will constrain through our analysis in the next section.The long baseline (d = 14.4 Mpc) offered by the high-energy neutrinos observed by IceCube can be used to constrain tiny values of δm 2 ∼ 10 −18 eV 2 , which are otherwise difficult to access in other sources.
Analysis -The investigation performed in this work closely follows the procedure applied in Ref. [4].Thus, we use an unbinned likelihood ratio method with a corresponding likelihood function defined as [31] log with S i and B i being the probability density functions for signal and background neutrino events, respectively.Further, n s denotes the number of signal events and N is the total number of recorded events.In Ref. [4] a dataset of track-like events created by muon (anti-)neutrinos with an exposure of 3186 days has been investigated.The analysis routines provided by the Ice-Cube collaboration [32] have been modified to account for the physics scenario of this work by adjustments of the likelihood function and the inclusion of a profile likelihood ratio test.In particular, the oscillation probability of pseudo-Dirac neutrinos, cf.Eq. 5, has been incorporated in the calculation of neutrino signal events n s for the likelihood function in Eq. 6.In addition, the likelihood function itself is added with a set of Gaussian pull terms in order to account for parameter correlations of the original analysis as well as experimental uncertainties of neutrino mixing parameters.For the latter, we took the latest global fit results [33,34] with their corresponding uncertainties.We assume normal neutrino mass ordering and for simplicity set δ CP ≃ 0 and σ x ∼ 10 −10 m -the latter is just representative of the width of the neutrino wave packet.From Eq. 4, it is clear that this choice does not lead to additional decoherence due to wavepackage separation on the Mpc-scale.Smaller values of σ x would imply a smaller coherence length, and hence averaging of the spectra.This would result in a loss of sensitivity.We chose a typical value of σ x , characteristic of charged-current production of neutrinos in AGNs, core-collapse supernovae (CCSNe) and other astrophysical objects [30].Marginalizing over σ x would reduce the sensitivity slightly since smaller values of σ would lead to decoherence.However, it is reasonable to expect that σ x cannot be smaller than 1 fm -typical scales associated with the size of a proton.The neutrino flux of NGC 1068 in muon and the antimuon flavor is assumed to follow a power-law of the form with Φ 0 being the overall flux normalization associated with an energy of 1 TeV and γ is the spectral index.At this point, the authors would like to emphasise that for the case of NGC 1068, the neutrino emission spectrum might be more complex [35]. 2 However, in order to follow the original analysis we apply the same unbroken power-law shape as the IceCube collaboration.A different spectrum, with bump/dip-like features might reduce the sensitivity of our results.Furthermore, a contribution from muons created by tau neutrino interactions is automatically taken into account by the provided Ice-Cube routines.For reasons of simplicity, we leave details of the neutrino production mechanisms aside and assume a generic neutrino flavor composition of (ν e : ν µ : ν τ ) = (1 : 1 : 1) arriving at Earth, which corresponds to a neu-trino admixture of (ν e : ν µ : ν τ ) = (1 : 2 : 0) at the source stemming from pion decay for standard neutrino oscillations.Deviating from these assumptions can lead to different compositions of neutrino mass states arriving at the Earth, and hence change the contribution to the survival probability in Eq. 5.This will have an impact on our results and require a dedicated discussion which is beyond the scope of this work.
The corresponding number of signal events is then given by [37] with the detector's lifetime t and its effective area A eff being a function of the neutrino energy E ν and the solid angle Ω.The overall neutrino oscillation probability P depends on the neutrino energy E ν , the active-sterile mass splitting δm 2 and the standard neutrino mixing angles s i ≡ sin θ i for i = {12, 23, 13}.Hence, in our analysis the number of signal events depends on the six parameters: n s ≡ n s (δm 2 , Φ 0 , γ, s 12 , s 23 , c 13 ).Existing correlations between Φ 0 and γ are accounted for by reproducing the covariance matrix of the original IceCube analysis and incorporating it with a two-dimensional Gaussian pull term to the likelihood function in Eq. 6.In doing so, we also ensure that the fit preserves the original best-fit values of Φ 0 and γ, taking into account a normalization correction due to the introduced oscillation probability in Eq. 5. Further, current knowledge of standard neutrino mixing parameters is included with individual Gaussian pull terms with uncertainties taken from Refs.[33,34].Note that for tiny mass-squared differences δm 2 ≲ 10 −22 eV 2 , the effects of active-sterile mixing become negligible for the energy region of interest, such that we coincide with the usual no additional oscillation case.In order to determine a limit on the mass-squared difference of pseudo-Dirac neutrinos, we determine a profile likelihood ratio as [38] where the numerator represents the likelihood function with a fixed δm 2 , while in the denominator all six quantities are tuned in the minimization procedure.In particular, double-hatted parameters like Φ 0 and γ are the conditional (maximum likelihood) estimators that depend on the chosen value of δm 2 , whereas the parameters occurring in the denominator are the estimators that maximize the unconstrained likelihood function.
For each δm 2 value under study, we determine the likelihood ratio in Eq. 9 and with asymptotic sampling distributions given in Ref. [38] we can assign each so- determined ratio a corresponding p-value, i.e. p δm 2 = ∞ q δm 2 ,obs dq δm 2 f (q δm 2 |δm 2 ), which we use as an exclusion criterion for the tested δm 2 values.Since the introduced oscillation probability in Eq. 5 affects the expected signal only in one direction, i.e. a decrease of the expected neutrino flux, we have to perform a one-sided statistical test.Following the IceCube collaboration, we assume that the test statistic f (q δm 2 |δm 2 ) used has approximately χ 2 -like behavior [32], except for modifications due to the onesided test case.For details about the modified profile likelihood ratio and the corresponding sampling distribution, we refer to Ref. [38].In this work, we are interested in a limit on δm 2 at 90% confidence level (C.L.) (1 − p δm 2 ).Thus, our results are determined by finding the value of δm 2 , which yields p δm 2 ≤ 0.1.The whole analysis has been performed within the Anaconda/SciPy framework [39,40], while the iminuit package [41] was used to minimize the likelihood function in Eq. 6.
Results -The fact that IceCube observed these neutrinos from NGC1068 implies that these did not oscillate into their sterile counterparts.In the presence of activesterile oscillations, the expected power-law flux will undergo spectral distortions, as shown in Fig. 1  straints on possible values of δm 2 .In our analysis, we assume that only muon (anti-)neutrinos produced at the source are detected as muon (anti-)neutrinos, i.e. we set α = β = µ in Eq. 5.This implies that we only detect ∼ 41% of all neutrinos.This leads to a slightly different flux normalization from the original best-fit values observed by IceCube: Φ 0 ∼ 1.2 × 10 −10 TeV −1 cm −2 s −1 , while γ basically remains the same.In principle, the contribution of a ν e produced at source and detected as ν µ is expected to be few tens of percentage of that of the original ν µ contribution.However, our results do not change much if we include this contribution.Our main result, depicted in Fig. 2, shows that the data from IceCube can be used to rule out δm 2 = 1.4 × 10 −18 eV 2 at more than 90% confidence level (p δm 2 ≤ 0.1).There is some sensitivity to even lower values of δm 2 , albeit at lower confidence.This agrees with the analytical estimate in Eq. 3, which shows the mass-squared difference sensitive to this oscillation length for TeV energy neutrinos.To the best of our knowledge, this is the strongest constraint to date on the smallness of δm 2 with such a high significance.
For smaller values of δm 2 , the oscillation length becomes larger than the distance to the source, hence sensitivity is reduced.
However, note that the situation is actually more complicated.Fig. 2 seems to indicate that data from NGC 1068 can be used to set bounds on arbitrarily large values of δm 2 .This is certainly not feasible since, for large δm 2 , rapid variations of the oscillation probability, and hence the flux, would be averaged out by the finite energy resolution of the detector.As a result, any reduction in the number of events would be compensated by a corresponding increase in the normalisation Φ 0 .This is the major drawback of a counting analysis with the number of events as in Eq. 8.A spectral analysis, on the other hand, would have been more sensitive to the upper limit on δm 2 .
In this analysis, Φ 0 and γ are allowed to vary within their correlation given by the covariance matrix corresponding to the original fit.Hence, although a reduction in the number of events by averaged oscillations can be compensated by an increase in Φ 0 , the correlation between Φ 0 and γ does not allow the Φ 0 to obtain values needed for a complete compensation.This is clear from the plot in Fig. 3, where we show the covariance matrix to pull Φ 0 and γ to the best-fit values corresponding to the original IceCube analysis and the Φ 0 variations during the fit procedure to compensate the reduced the number of events due to oscillations.The change in Φ 0 is not enough to compensate for the averaging effect for high δm 2 .
To simulate the loss of sensitivity at large masssplittings, we have performed the following test.We check for values of the mass-splitting where the survival probability, and hence the total number of events, flattens out, i.e. does not change the number of events below 1%.In this regime, the effect can be mimicked by a change in the normalization of the spectra only.As a result, sensitivity would be lost in this regime.This yields the sensitivity region to be 10 −19.1 ≲ δm 2 ≲ 10 −17.0 .This is explained through Fig. 4, where the left panel shows the neutrino flux variations at the Earth for two different choices of mass-splittings within which we expect the flux variations to be relevant.The right-hand panel shows the variations of the total events, indicating the values of δm 2 for which the event counts do not change by more than 1%.
Fig. 5 shows the correlation between δm 2 and γ for a fixed flux normalization Φ 0 = 1.2 × 10 −10 TeV −1 cm −2 s −1 .Analogous to our determined limits, we indicate the 90% C.L. contour for both parameters.Correlations between Φ 0 and γ are considered via an appropriate pull term in the applied likelihood function.The figure indicates that the data is clearly consistent with the "no-oscillation" hypothesis, which corresponds physically to oscillations with longer baselines.This shows that the addition of these oscillations does not help improve the fit, rather the data can be used to put constraints on these oscillations.
Finally, note that uncertainties of the source distance have not been accounted for in this analysis.Since we are varying δm 2 , there will always be a certain range of δm 2 for which the neutrino oscillation probability will undergo a few oscillation cycles.This will not get averaged out due to the uncertainty in the distance.Of course, this might lead to a slight reduction in sensitivity, but that is subdominant with respect to the effect of the other uncertainties like the determination of the neutrino energy.
A similar argument can be used to put bounds on δm 2 from observations of a neutrino event of energy E ν ∼ 290 TeV from the blazar TXS 0506 + 056 at a redshift of z = 0.33 ± 0.0010, corresponding to a distance d = 1.3 Gpc [2,3].A naive estimate yields which is stronger due to the larger energy of the neutrino.Additionally, IceCube later reported a series of 13 ± 5 events from the same direction.These observations can be used, in principle, for a similar study to constrain values of δm 2 .
Conclusion -Detection of high-energy neutrino events by IceCube has opened up new frontiers in multimessenger astronomy.In particular, the naturally long baseline available to the neutrinos, coupled with their high energies, allows them to be the harbingers of new exotic physics, otherwise inaccessible to mankind.In this work, we used the IceCube observation of ∼ 79 TeV-sh neutrino events coming from the direction of the active galactic nuclei NGC 1068 to probe the extent of lepton number violation in the Standard Model.
A soft lepton number violation can be manifested through active-sterile neutrino oscillations over astronomical baselines inversely proportional to the masssquared difference (δm 2 ) between active and sterile neutrinos.These oscillations, if present, will lead to a distortion of the event spectra, and a reduction in the number of neutrinos observed.
This simple yet powerful idea can be utilized to probe the strength of lepton number violation, characterised by the magnitude of δm 2 .The identification of the source of these neutrinos, and precise measurement of their energies, allowed us to precisely calculate the oscillation lengths associated with δm 2 , and the corresponding spectral modifications.We used the analysis of the IceCube collaboration, outlined in Ref. [4], to rule out δm 2 = 1.4 × 10 −18 eV 2 at more than 90% confidence level, using the total number of events observed.We simulated the loss of sensitivity at large mass-squared splitting by checking where the number of events does not change by more than 1% when the masssplitting is changed.This yields the sensitivity region as 10 −19.1 ≲ δm 2 ≲ 10 −17.0 eV 2 .
Parameter correlations of the original IceCube analysis as well as uncertainties of the neutrino mixing angle have been properly treated via Gaussian pull terms.The result obtained is one of the strongest bounds on the extent of lepton number violation from experiments.However, we would like to point out that the performed investigation depends crucially on the assumptions of the neutrino energy spectrum as well as their flavor composition.We leave a detailed study on the impact of various flux assumptions as well as a more involved sterile neutrino sector, i.e. different δm 2 i , for future analyses.These extreme new physics scenarios are clearly inaccessible to any terrestrial experiments.The observations of these point-sources of astronomical neutrinos like NGC 1068 and TXS 0506+056 clearly pave the way for new laboratories for neutrino physics.

FIG. 3 .
FIG.3.Top: the covariance contour to pull Φ0 and γ to the best-fit values equivalent to the original IceCube analysis.Bottom: the Φ0 variations during the fit procedure to compensate for the reduced number of events due to oscillations.

δm 2 = 10 − 19 .1 eV 2 δm 2 = 10 − 17 .0 eV 2 FIG. 4 .
FIG. 4. Left: Neutrino flux variations at the Earth due to active-sterile oscillations as a function of the neutrino true energy Eν , ( rescaled to match the flux measured by IceCuce).The standard case (P = 1) is compared with two different cases δm 2 = 10 −19.1 eV 2 and δm 2 = 10 −17 eV 2 .Right: The corresponding events observed in IceCube.The vertical lines show the benchmark values of δm 2 from where the event counts do not change by more than 1%.
2changes, minima are reached accordingly when L osc takes values close to the distance between the Earth and NGC 1068.This is the region where we expect maximum sensitivity in our results.A profile likelihood analysis can be used to set con-