Resonant Neutrino Self-Interactions

If neutrinos have self-interactions, these will induce scatterings between astrophysical and cosmic neutrinos. Prior work proposed to look for possible resulting resonance features in astrophysical neutrino spectra in order to seek a neutrino self-interaction which can be either diagonal in the neutrino flavor space or couple different neutrino flavors. The calculation of the astrophysical spectra involves either a Monte Carlo simulation or a computationally intensive numerical integration of an integro-partial-differential equation. As a result only limited regions of the neutrino self-interaction parameter space have been explored, and only flavor-diagonal self-interactions have been considered. Here, we present a fully analytic form for the astrophysical neutrino spectra for arbitrary neutrino number and arbitrary self-coupling matrix that accurately obtains the resonance features in the observable neutrino spectra. The results can be applied to calculations of the diffuse supernova neutrino background and of the spectrum from high-energy astrophysical neutrino sources. We illustrate with a few examples.


I. INTRODUCTION
While the interactions between neutrinos is extraordinarily feeble in the Standard Model, there are a number of reasons to entertain the possibility that new physics may introduce stronger neutrino self-interactions [1][2][3][4]. Astrophysical neutrinos may provide a powerful probe in the search for such self-interactions (νSI) [5]. Strong features, such as dips and enhancements, can be imprinted on astrophysical-neutrino spectra, which when analyzed can yield νSI parameter values. In particular, there is the diffuse supernova neutrino background (DSNB) and a collection of high-energy astrophysical neutrino (HEAN) sources. The DSNB is the isotropic time-independent flux of neutrinos and anti-neutrinos around tens of MeVs emitted by distant core-collapse supernovae [6]. These diffuse sources come from distances around 10 Mpc [7,8] up to a redshift of 5 with a peak around a redshift of 1 [9]. In comparison, although HEAN sources have no identified production mechanism, they are observed by IceCube to follow a power law [10].
If there are neutrino self-interactions, then interactions of DSNB and/or HEAN neutrinos with low-energy (∼ 0.001 eV) cosmic-background neutrinos can appear in the obseved DSNB and/or HEAN spectra as absorption features or enhancements at lower energies. The calculation of the observed flux is straightforward but, most generally, involves solving a series of coupled integrodifferential equations that evolve the energy and flavor * Electronic address: creque@jhu.edu † Electronic address: kamion@jhu.edu ‡ Electronic address: jhyde@bowdoin.edu distribution of neutrinos. These equations describe the injection of neutrinos from sources, the redshifting of neutrinos, the absorption from self-interactions, and the reinjection of lower-energy neutrinos after such interactions. The solutions to these equations require either a Monte-Carlo simulation or a straightforward, but computationally intensive, numerical integration of the equations. This thus limits the regions of νSI parameter space that can be investigated. For example, previous work assumed a universal self-coupling matrix [11][12][13], a diagonal [14] self-coupling matrix, no flavor dynamics [15,16], or particular values for the mediator mass, coupling constant, and neutrino mass [17][18][19].
In this paper, we present an analytic approach to resonant astrophysical-cosmic neutrino scattering for arbitrary self-coupling matrix and neutrino number. This solution is built on the observation that most of observable effects explored in DSNB/HEAN studies arise from resonant neutrino self-interactions. To illustrate the utility of the approach, we then use this solution to explore the discovery space for a model of τ -self-interactions (relevant for the DSNB) and another for sterile-neutrino selfinteractions (relevant for the HEAN). This paper is organized as follows. In Section II, we review the general formalism of neutrino mixing and transport. In doing so, we present three formal solutions: first in the case of no interactions; second in the case of only absorption interactions; and third in the case of both absorption and reinjection. We formulate explicit solutions in Section III where we specify the nature of neutrino self-interaction. We begin by considering only a single species of neutrino. Then, we generalize this result and present a solution for arbitrary neutrino number and selfcoupling matrix. Finally, in Section IV we use our new solution to identify regions of parameter space that can be accessed. Specifically, we consider DSNB probes of selfinteractions with Super-Kamiokande and HEAN probes with IceCube. We discuss and conclude these results in Section V and Section VI.

A. Neutrino Mixing
Neutrinos can be represented in either the mass basis or the flavor basis. In what follows, Greek indices are reserved for flavor states and Latin indices for mass states. In order to switch between the two bases, the neutrino mixing matrix U must be used according to where the sum is over all mass states and with U unitary. Unitarity implies that for any flavor state α, i |U iα | 2 = 1, and so |U αi | 2 is interpreted as the probability that flavor α is observed as mass state i, or vice verse.

B. Neutrino Transport
The specific flux Φ i (t, E) of astrophysical neutrinos ν i (number of astrophysical neutrinos per unit conformal time per unit comoving area per unit energy) at cosmic time t and observed energy E obeys the Boltzmann equation where H(z) = H 0 E(z) is the Hubble parameter at redshift z, with z as a proxy for t, and for ΛCDM, with Ω m the matterdensity parameter today. Moreover, S i is the production rate of astrophysical neutrinos ν i , Γ i is the absorption rate of astrophysical neutrinos due to neutrino scatterings, and S tert,i is the tertiary source term accounting for the possible reinjection of astrophysical neutrinos postscattering. Note that by this definition of the specific flux, the comoving number density of astrophysical neutrinos is (1/c) dE Φ i (t, E), different from Refs. [20,21] where they considered Φ i to be defined to reproduce the physical number density. Since neutrino decoherence time scales are smaller than any other relevant time scale, the transformation Φ α (t, E) = i |U αi | 2 Φ i (t, E) can be performed to switch back to the flavor basis at any point. Moreover, if the source terms are given in the flavor basis S α , then the mass basis source term is S i = α |U αi | 2 S α . In the absence of interactions, S tert,i = Γ i = 0, the solution of Eq. (2) is obtained by identifying the total time derivative as d/dt = ∂/∂t + (dE/dt)∂/∂E, with dE/dt = −HE, leading to with a(t) the scale factor at time t. The factor of [a(t)/a(t )] inside the source term accounts for the redshifting of the energy from t to t, while outside the source term it accounts for the redshifting of the differential energy from t to t. The addition of a nonzero sink term while neglecting any reinjection, Γ i = 0, S tert,i = 0, does not complicate things much further as then with τ i (t , t, E) the optical depth of a neutrino ν i of energy E between times t and t. As a result, astrophysical neutrinos at time t not only go through the previous redshifting, but now travel through a medium of optical depth τ i from the emission time t to the observed time t.
Formally, if neutrino reinjection is taken into account, S tert = 0, a solution is easily written down, However, since the tertiary source is a function of the specific flux itself, the solution is in general not closed. Therefore, if we are to move forward, a particular model must be specified.

A. Single Neutrino Species
The particular neutrino model we consider at first is that of a single species of self-interacting neutrinos ν of mass m ν whose self-interactions are mediated by a scalar particle φ with mass m φ and coupling strength g. We ignore the existence of other neutrino species, and as such we suppress any indices present in relevant equations. That is, we initially consider the interacting Lagrangian, If this is the case, then astrophysical neutrinos will scatter with cosmic neutrinos, causing depletion of astrophysical neutrinos at a resonant energy E R = [m 2 φ /(2m ν )]c 2 at a rate Γ(t, E) = n ν (t)σ(E)c. Here and onwards, cosmic neutrinos will refer to cosmic-background neutrinos. We define n ν (t) to be the physical number density of our single cosmic neutrino species, σ(E) the scattering cross section of the process νν → νν, and c the speed of light. After depletion, neutrinos are then re-injected at energies E < E R . We take the scattering cross-section to have a Breit-Wigner form where is Planck's constant, s = 2Em ν c 2 , and Γ φ = g 2 m φ c 2 /(4π) the decay width. If the width of the resonance is small enough, resonant scattering can be approximated by a Dirac delta function. We now quantify when this occurs. The width of the resonance is where . For a detector with resolution ∆E, a width cannot be resolved and thus is a delta function ∆E. Therefore the coupling must satisfy where ∆E ≈ 1 MeV for a detector such as Super-K [22], and E R ≈ 25 MeV for masses m φ c 2 = 1 keV, m ν c 2 = 2× 10 −2 eV. We conclude that unless the coupling is of order unity, which is most of the available parameter space [15], a detector will not be able to resolve the resonance and we approximate the cross section as a delta function. A nascent delta function in the Breit-Wigner form is δ(x) = lim →0 (1/π) /( 2 +x 2 ), so that the resulting cross section is Resonant scattering is isotropic when the mediator is a scalar field, and therefore the differential cross section dσ(E 1 , E 3 )/dE 3 , where an incoming neutrino with energy E 1 scatters to an outgoing neutrino of energy E 3 , has a flat distribution With this form we now evaluate the tertiary source for neutrino production in Eq. (2).
In our case of cosmic neutrino upscattering, two neutrinos are re-injected after an initial neutrino is taken from the sink term, and cosmic neutrinos have energies much smaller than supernova neutrinos so their relative velocity is the speed of light. Thus the tertiary term takes the following expression, converting our initial differen-tial equation into an integro-differential equation With our delta-function approximation, this term is now evaluated as with Γ R (t) = n ν (t)σ R c and Θ(x) the Heaviside function with Θ(0) = 0. Moreover, we simplify the optical depth as with z a proxy for t , Thus, for E ≥ E R the spectrum is the same as a nointeraction Boltzmann equation, and so is solved in the same manner. However, when E < E R , neutrinos are reinjected at twice the rate of their depletion at the resonant energy. As such, the expression for neutrino reinjection still requires solving for the specific flux at E = E R and plugging it back in for evaluation at lower energies, which at first makes Eq. (13) seem not closed. However, Eq.
(2) has a delta function via the absorption term Γ(t, E)Φ = n ν (t)σ(E)cΦ at the resonant energy, and so we must obey the boundary condition at this point. In order to satisfy this condition, we integrate Eq. (2) around the resonant energy from below the resonance Again, above the resonant line the optical depth of freestreaming with no interactions is zero, while below, it is As a result, Eq. (13) has a closed form expression. The purpose of retaining E − R in this expression is as a reminder that in order to evaluate τ (t , t, E − R ) one must take the left-side limit of the integral in the expression of τ . Since cosmic neutrinos are low-energy, astrophysicalcosmic neutrino scattering does not add or remove energy from the astrophysical neutrino spectra, but only redistributes it. We have checked both analytically and numerically that Eq. (13) obeys this condition. There are two differences in the expression for Φ between Eq. (13) and Eq. (15). First, is the presence of the factor 1 − e −τ rather than e −τ . This factor can be understood as follows: after an astrophysical neutrino redshifts through a resonance over a short period of time, the specific flux at the resonant energy only has a fraction e −τ remaining of the original flux. It follows then that the amount that is injected at lower energies must be the complementary fraction, 1 − e −τ . The second difference is the factor of H(t)/Γ R (t), which changes the rate of injection from Γ R (t) in Eq. (5) to H(t). This change in the rate of injection is due the resonance line redshifting in time. If the scattering rate is faster than Hubble, then the flux of neutrinos at the resonant energy is suppressed. Conversely, if the rate is slower, then the scatterings have little effect.

B. Multiple Neutrino Species
In the presence of multiple neutrino species, the previous equations do not hold. Here we present the analogous calculation with additional neutrinos, taking the mass of each neutrino species to be m j with corresponding cosmic physical number density n j . The interaction Lagrangian term for the most general mass-basis interaction with a scalar mediator φ of mass m φ is given by with g ij the self-coupling matrix. In this model, astrophysical neutrinos scatter off of one of any of the cosmic neutrino species, causing depletion of astrophysical neutrinos at corresponding resonant energies The cross section σ ij ≡ kl σ ijkl is the sum of scattering cross sections for the processes ν i ν j → ν k ν l . We take σ ijkl to have the Briet-Wigner form with s j = 2Em j c 2 , and Γ φ = ij |g ij | 2 m φ c 2 /(4π) the decay width. Note that the decay width has changed since there now exists multiple decay branches for φ. After depletion, the neutrinos are re-injected at a rate ac-cording to with δ il the Kronecker delta function that accounts for the possibility of upscattering into two astrophysical neutrinos of the same state rather than just one. That is, compared to the single neutrino species case, this expression accounts for production of neutrinos of type i from an astrophysical flux Φ j hitting a cosmic neutrino density n k , with i, j, k not necessarily all being the same. Once again, the differential cross section takes a flat distribution Moreover, using our delta function limit, the cross section takes the form . As a result the tertiary source term is In addition, the optical depth is Analogous to before, we satisfy the boundary condition around each resonance by the conditions Now, it is true that only above the highest resonant line the optical depth is zero. Thus, the general solution is Then when we want to convert back to the flavor basis we use the neutrino mixing matrix once again. Eq. (21) is our main result that describes the propagation of multiple astrophysical neutrinos species that self-interact arbitrarily with cosmic neutrinos. If a flavor self-coupling matrix is given instead, the identification g ij = αβ U iα U jβ g αβ leads to an easy substitution. We present the analogous equations with this substitution in Appendix A. Note however that in order to solve Eq. (21) in a closed manner, the highest resonant boundary condition must be solved for first, as it is a function of only the source S i and scattering rate Γ i . This is to be contrasted with the boundary conditions for lower resonant energies, which depend not only on these quantities, but also the flux at higher resonances. This dependence arises because a neutrino can be absorbed at a resonance E j , downscattered to an energy E > E k , with E k < E j some other resonance, and then be redshifted down to E k . In this way, astrophysical neutrinos may cascade down the resonance pipeline until they reach energies below the lowest resonant energy.

IV. NUMERICAL RESULTS
Given these analytic results, a wealth of neutrino selfinteractions can be explored and constrained. However, due to the large dimensionality of the general problem, we narrow our scope to two specific models. Moreover, many factors aside from neutrino-self interactions can affect the resulting spectrum, such as detector backgrounds and energy thresholds. Such a detailed analysis, however, is outside the scope of this paper. That is, we consider only a single source of neutrinos with shot-noise error. Specifically, first we consider the standard 3-neutrino model, adjoined with a τ self-interaction coupling constant g τ τ .
Interactions of this form have been proposed to resolve the Hubble tension [23], although our analysis does not rely on this explanation.
Second, we add a sterile neutrino to the 3-neutrino model, along with a sterile s self-interaction coupling constant g ss . This case is motivated by the LSND, Mini-BooNE and reactor anomalies which suggest mixing with eV-scale sterile neutrinos [27,31,36]. While such mixing would be in tension with Planck, self-interactions of the sterile neutrino by a mediator of mass m φ MeV would bring results back into harmony [37,38]. In both models we consider all other neutrino self-coupling constants are taken to be zero.
In order to apply Eq. (21) to the four-neutrino case, we need to choose definite values for the mixing matrix elements. We use a standard parameterization [32] with R c (a, b) a 4×4 rotation matrix with matrix elements R c (a, b) ij and a mixing angle θ ab . The matrix elements are those of the 4 × 4 identity except for the following submatrix: where s ab = sin(θ ab ) and c ab = cos(θ ab ), and δ CP ab is a complex CP violating phase. In general there are also Majorana phases associated with the mixing matrix, but since since we considering lepton-conserving processes, we neglect them [34].
In the limiting case of no mixing between active and sterile neutrino states, θ 14 = θ 24 = θ 34 = 0, each of R 1 (3, 4), R 0 (2, 4), R 1 (1, 4) is the identity and we obtain the standard 3-neutrino mixing matrix [33] plus a decoupled sterile state. Note that in this model we consider self-interactions only among sterile neutrinos, so in this no-mixing limit our astrophysical spectra will return to the standard expectation, regardless of the value of g ss . Motivated by the short-baseline anomalies, we take θ 14 = θ 34 = 0 and sin 2 (θ 24 ) = 0.1, so θ 24 = 0.161.
In addition to the mixing matrix, the neutrino mass spectrum m is also constrained. We first review the constraints on the lightest three neutrinos. Oscillation experiments give the value of two mass-squared differences [47]. As a result, it is unclear whether the neutrino mass spectrum follows a normal hierarchy (NH) m 1 < m 2 < m 3 or inverted hierarchy (IH) m 3 < m 1 < m 2 . However, a lower bound on the neutrino masses is obtained by setting m 1 = 0 in the NH and m 3 = 0 in the IH. In addition, an upper bound is obtained from Planck [35], as it constrains the sum of neutrino masses to be such that j m j c 2 < 0.12 eV. As a result, the following table of neutrino mass constraints can be made The mass of the mediator is chosen to correspond to the energy ranges dictated by the sources we choose. That is, for an experiment that measures spectra between neutrino energy ranges [E min , E max ], the range of mediators that can be be probed is E min ≤ [m 2 φ /(2m j )]c 2 ≤ E max for each neutrino mass m j ∈ m. In order to simplify our analysis we only compare the null hypothesis with our model at the resonant energies, and not the entire spectrum. Finally, we choose a fixed bin size for each constraint.
We denote the event count under the null hypothesis g ij = 0 by N null . Thus, assuming only Poisson shot noise error, we find that the number of events N events in each bin can be measured away from the null hypothesis with a signal to noise Therefore, the number of events needed to distinguish from the null hypothesis N events = N null is with (S/N ) ± = ±|(S/N )|. N ± is also known as the (S/N )-σ uncertainty in the measurement of N null , with N + the upper and N − the lower uncertainty. Again, for our analysis, we only use N − when looking for depletions.

A. DSNB
The production rate of neutrinos per comoving area per unit time per unit energy from core-collapse super- [39], with R CCSN the CCSN rate per comoving volume, and dN i /dE the number spectrum of neutrinos of type i emitted by one supernova explosion. For R CCSN we use the parameterization of Ref. [9] with the lower bound of the Salpeter initial mass function. Moreover, we assume equipartition of energy among neutrino species and thus approximate the spectrum of one neutrino species by a Fermi-Dirac distribution with zero chemical potential [40] with E tot = 3 × 10 46 J the total energy in neutrinos emitted by the supernova and 4 MeV ≤ k B T SN ≤ 8 MeV the supernova temperature [6]. We plot two possible flux spectra Φ e of electron anti-neutrinos from the DSNB in Fig. 1 for T . While a supernova temperature k B T SN = 8 MeV is disfavored, it does not heavily alter our conclusions. For the heaviest normal neutrino mass hierarchy, three resonances are potentially observable when flavor self-interactions are considered. While not observed yet, the addition of gadolinium sulfate to large water Cerenkov detectors would allow for the discrimination, and thus detection, of DSNB events from spallation and atmospheric neutrino events [41,42]. Electron anti-neutrinos in the DSNB are in the correct energy regime to be detected through inverse beta decay scattering at Super-Kamiokande [22]. Specifically, through the processν e p → e + n, DSNB anti-neutrinos collide with water molecules in Super-Kamiokande, producing a positron that emits Cherenkov radiation that is detectable. As a result, the colliding anti-neutrino must have minimum energy E min ν = m e c 2 + ∆ = 1.806 MeV, with ∆ = (m n − m p )c 2 . In the following, we use Eq. (25) of Ref. [43] for the inverse beta decay cross section σ IBD . The number N events of events detected in a positron energy bin [E e + , E e + + δE] is then (27) with T the time of observation and N p the number of scattering targets. Note that in this expression E + ∆ is the neutrino energy, while E is the positron energy.
We show the event counts and uncertainties corresponding to Fig. 1 in Fig. 2. Comparing our null hypothesis to our model at the resonant energies, we obtain the forecasted 1σ constraints in Fig. 3.

B. High-Energy Astrophysical Neutrinos
The production rate of high-energy astrophysical neutrinos per comoving area per unit time per unit energy is L(z, E) = W(z)L 0 (E), with L 0 the differential number luminosity for each source and W the redshift evolution of the source density. We take the redshift evolution to follow the star-formation rate, W(z) = R CCSN (z). Moreover, following IceCube's 6-year data analysis [10], we take the differential number luminosity to be a power law L 0 ∝ (E/E 0 ) −γ . We plot two possible flux spectra Φ e of electron anti-neutrinos in Fig. 4. The number of events observed by IceCube is [44] N events = T Ecasc+δE Ecasc dEΦ e (E)A eff (E), (28) with T the time of observation and A eff (E) the IceCube effective area, which we take from Ref. [45]. Note that we are approximating the neutrino energy to be the cascade energy.
We show the event counts and uncertainties corresponding to Fig. 4 in Fig. 5. Comparing our null hypothesis to our model at the resonant energies, we obtain the forecasted 1σ constraints in Fig. 6.

V. DISCUSSION
Several points are worth examining in further depth. First, Eq. (21) only holds when each cosmic neutrino species is cold. However, we know there exists at least one cold neutrino species. Therefore, in the case where one or more cosmic neutrino species are not cold, this equation is modified so that any sum over neutrino scattering cross sections is only over all cold species. Moreover, interactions with the lower-mass neutrinos should be suppressed relative to that of the heavier cold species due to thermal broadening.
Second, the spectra shown all have three resonances, and this not need be the case. The heaviest allowed normal neutrino mass hierarchy is special in this case, since the nearly-degenerate pair have masses much larger than any neutrino mass splitting. In the inverted scenario, this cannot be the case and so at most two resonances could be seen in any spectra that does not have a large energy range. If a single resonance is seen it is unclear how to disentangle the two scenarios, but such distinction is outside the scope of this work.
Third, while we made our constraints by looking for absorption features, one in principle could also look for enhancements in spectra. In experiments, it is simple as one needs to look for when the signal surpasses some threshold for statistical significance, N events > N + . However, it is less obvious theoretically what bins or how many bins one should look at in order to create a forecasted constraint from enhancements in a time-efficient manner. It  depends on the number of resonances, the shape of the null hypothesis spectrum, and the detection method. For example, in the DSNB the enhancements are much more pronounced at low energy compared to HEAN sources, since at low energies the DSNB spectrum falls off while the HEAN source spectrum grows.
Fourth, we wrote down our formulas assuming a single scalar φ, however it is straightforward to generalize to multiple scalars φ k with self-coupling matrices g k ij . The only possible subtlety is if degeneracies in the resonances occur, in which case the resonant condition needs to be altered accordingly by a sum over degenerate resonances.
Fifth, while this paper is focused on neutrino selfinteractions, it is also straightforward to incorporate arbi-trary resonant scattering between any species. The most obvious other cold species to generalize to would be cold dark matter.
Finally, we took the noise to be only Poissonian and assumed fiducial astrophysical parameters. In a realistic experiment, other backgrounds must be taken into account as well as degeneracies with their parameters. However, such a proper treatment, similar to Refs. [46,48], is outside the scope of this work.

VI. CONCLUSION
In this paper we have have considered the consequence of beyond the Standard Model neutrino self-interactions on various astrophysical neutrino spectra. We began by presenting the necessary formalism for neutrino mixing and transport. We did this not only to establish notation, but also in order to demonstrate that neutrino reinjection is a problem that is generally not closed.
In order to overcome this hurdle, we then took the limit where the scattering cross section goes to a delta function, and found that the former partial integrodifferential equations turn into a standard partial differential equation with simple boundary conditions. As a result, we then presented the solution for astrophysical neutrino spectra for a single neutrino species, following with one for an arbitrary number of neutrino species. These solutions were specified in either the mass basis or the flavor basis.
From this, we then demonstrated the utility of the analytic solution by considering our astrophysical sources to be either the diffuse supernovae background or highenergy astrophysical neutrinos. From there, we established forecasts and constraints on a normal 3-neutrino hierarchy with τ self-interactions, as well as a 4-neutrino hierarchy with sterile self-interactions. None of these calculations took a significant amount of time, and were routine in their implementation.
It will be interesting to implement this calculation in future work to explore the effects of neutrino selfinteractions on DSNB and HEAN spectra for a wider range of models that involve new neutrino interactions.
We consider the most general flavor interaction for a single scalar mediator φ of mass m φ L flavor The identification of g ij = αβ U iα U jβ g αβ allows us to use Eq. (21). We reparameterize the scattering rate with Γ k,αβγδ R = |U βk | 2 n k (t)σ αβγδ R c and σ αβγδ R = |g αβ | 2 |g γδ | 2 /[4(m φ c 2 )Γ φ ]. We choose such a reparameterization in order to separate the neutrino conversion probabilities from the scattering cross sections. In doing so, and invoking unitarity, we obtain the result dt [a(t)/a(t )]e −τi(t ,t,E) S i {t , [a(t)/a(t )E]}, with Γ ij R (t) = α |U αi | 2 βγδ Γ j,αβγδ R . While we have presented here these equations in the flavor basis for analytic insight, we note that in general it is easier numerically to use Eq. (21) with the appropriate substitution, as it contains less summations.