1 Introduction

The detection of high energy neutrinos of astrophysical origin is providing a new insight into the sources themselves and it can also help to probe different scenarios and properties of neutrino related physics [1]. In particular, neutrino inteactions with dark matter (DM) has been considered and their effects have been investigated under different assumptions [2,3,4,5,6,7,8,9,10,11,12,13,14,15].

One interesting possibility is that dark matter can be composed of ultralight particles, as this could help to alleviate problems of cold dark matter models regarding the overproduction of both substructure in the galactic haloes and satellite dwarf galaxies that are not observed [16]. Such extremely light particles can be pseudoscalars which arise naturally in scenarios beyond the Standard Model (SM) based on string theories [17] and are known as ultralight axion-like particles (UALPs). In particular, these models predict the existence UALPs spanning a wide range of masses, \(\sim 10^{-33}-10^{-10} \,\mathrm{eV}\) [18, 19]. While the lightest UALPs could be associated to dark energy [20, 21], the appropriate mass for viable DM candidates is \(m_a\sim 10^{-22}-10^{-21}\,\,\mathrm{eV}\), which has the advantage of preserving the large-scale behavior of the standard cold dark matter models, but at the same time it suppresses the substructure at shorter scales due to a large de Broglie wavelength \(\sim 1\,\mathrm{kpc}\) [16, 22].

In the present work, we study the effects of neutrino interactions with such particles during the propagation from extragalactic sources assuming that their evolution with the redshift follows the star formation rate (SFR) accoding to Ref. [23]. We follow the procedure discussed on a previous work with a different type of interaction [4]. The interaction vertex analysed in this work was considered previously in Ref. [9] assuming that UALPs are neutrinophilic, i.e., their interactions with the charged leptons can be neglected (see also Ref. [30]). In these conditions, the stringent bounds on the coupling to electrons based on stellar evolution considarations for globular clusters \(g_{ae}\lesssim 4\times 10^{-9}\,\,\mathrm {GeV^{-1}}\) [32] do not affect the coupling to neutrinos. The corresponding astrophysical and cosmological constraints that can be applied in the present context are much less restrictive, as discussed in Ref. [9]. In their work, they also described the effects of modified neutrino oscillations in baseline experiments such as DUNE [24] due to an effective potential included in the Hamiltonian. Our approach here is complementary, since we consider a flavor diagonal coupling \(g_{\alpha \beta }=g_{\nu a}\delta _{\alpha \beta }\) and this leads only to a global phase that implies no modification to oscillations due to interactions with UALP DM, and therefore no effects to be observed in baseline experiments.

Nevertheless, neutrino-UALP scatterings can still play a role and imprint effects on a diffuse neutrino flux. We explore this possibility in the present work by solving a transport of each mass (or propagation) eigenstate, for which the neutrino mass is well defined. The surviving fluxes of flavored neutrinos at the Earth are recovered by superposition, and this can be compared with observational data. We consider the cases where the hierarchy of neutrino masses are accommodated following a normal ordering (NO) or an inverted ordering (IO), and we explore three different scenarios of dominant neutrino sources with different initial flavor compositions. We find that interactions produce a change in the flavor composition of neutrinos arriving on Earth for \(g_{\nu a} > rsim 0.5\,\,\mathrm{GeV}^{-1}\) in comparison to the outcome expected with no interactions. In particular, the electron(muon) flavor comes to dominate in the NO(IO) case for neutrino energies above \(\sim 10^5\,\,\mathrm{GeV}\). Additionally, a directional dependence on the neutrino flavor ratios can be observed, given that galactic DM is more abundant in directions closer to the galactic center.

By comparing the flavor ratios integrated on the neutrino energy with available flavor measurements by IceCube[25, 26], the realization of the interaction effects studied seems unlikely, since the best fit values of the flavor ratios are not close to the predicted in the case of interactions. However, more data with more observation time and larger detectors are necessary to achieve a much higher precission in the flavor composition measurements [27, 28] and this will help to constrain definitely scenarios as the discussed here.

The rest of this work is organized as follows. In the next section, we compute the cross sections for the neutrino-UALP interaction, and in Sect. 3, we describe the neutrino propagation first in extragalactic space and then through the galactic DM halo. In Sect. 4, we present the results obtained for the neutrino fluxes and flavor ratios on Earth, and finally in Sect. 5, we conclude with a discussion.

2 Neutrino interactions with ultralight axion-like particles

In the present work, we consider the interactions between neutrinos and UALPs corresponding to the following Lagrangian term (e.g. [9, 29]):

$$\begin{aligned} \mathcal {L}_{\nu _\alpha a}= -ig_{\alpha \beta }\partial _\mu a \bar{\nu }_{\alpha }\gamma ^\mu \gamma _5 \nu _\beta , \end{aligned}$$
(1)

where \(\alpha ,\beta =\{e,\mu ,\tau \}\). This effective interaction involves only the active neutrinos, and it is implicitly assumed that the UALP coupling to the charged leptons is not relevant, which is a characteristic feature of neutrinophilic models (e.g. [9, 30, 31]). Without going into details, we note that this kind of interactions appear in models such as the applied in Ref. [33], where a new scalar doublet and a new global symmetry U(1) are included in order to give small Dirac masses to the neutrinos. In this type of models, known as neutrinophilic two-Higgs doublet models (nu2HDM), it is possible to obtain the type interactions sought, i.e., with a strong coupling to the neutrinos. Although they were originally proposed for the QCD axion [34], a similar approach could be suitable in our context for UALPs. We also focus on a minimal scenario where the UALP has flavor-conserving couplings at tree level, as we discuss below.

The interaction can be decomposed into terms corresponding to each mass eigenstate \(\nu _i\), given that \(\nu _\alpha =\sum _{i=1}^3 U_{\alpha i}\nu _i\), where \(U_{\alpha i}\) are the elements of the unitary mixing matrix for neutrinos [37, 38]. Defining \(q:= \bar{p}+p\) as the momentum carried by the UALP a, we have that \((\partial _\mu \, a) =i q_\mu \, a\), and making use of Dirac equation, it follows that

$$\begin{aligned} (\partial _\mu a )\bar{\nu }_i\gamma ^\mu \gamma _5\nu _j= & {} - i\left( m_i + m_j\right) a\bar{\nu }_i\gamma _5 \nu _j. \end{aligned}$$
(2)

We then consider the case of a flavor-diagonal coupling, where we have \(g_{\alpha \beta }=g_{\nu a}\delta _{\alpha \beta }\) for all the flavors. Then, Eq. (1) can be rewritten as

$$\begin{aligned} \mathcal {L}_{\nu a}= -ig_{\nu a}\sum _{ij} \left( m_i + m_j\right) \sum _\alpha U^*_{\alpha i}U_{\alpha j}a\bar{\nu }_i\gamma _5 \nu _j, \end{aligned}$$
(3)

and since \(\sum _{\alpha }U^*_{\alpha i}U_{\alpha j}= \delta _{ij}\), we have the following interaction terms for each massive neutrino \(\nu _i\):

$$\begin{aligned} \mathcal {L}_{\nu _i a}= & {} -ig_{i}\, a\bar{\nu }_i\gamma _5 \nu _i, \end{aligned}$$
(4)

where \(g_{i}= 2 m_i g_{\nu a}.\)

Since the coupling is flavor-diagonal and flavor-universal, this leads to an effective potential proportional to the identity. When this is added to the Hamiltonian to study the effects of neutrino oscillations in a DM background (e.g., as in Refs. [9, 11]), it only contributes with global phase which does not modify the relative energies of the different massive neutrinos. Hence, as mentioned above, no effects due to oscillations in DM medium are expected in the present context. However, the effect of neutrino scattering with UALPs can still be relevant and lead to an absorption of the neutrino flux at a given energy, as was also considered for other interaction vertices in previous works [2, 4, 6, 39].

The fact that consider only flavor-diagonal and flavor-universal interactions also allows to avoid very strong constraints for neutrinos decays [35]. Still, these asumptions correspond to a possibility which is usually explored in similar studies (e.g. [35, 36]), and in particular, the only relevant effect is due to neutrino-DM scattering. Hence, if off-diagonal couplings are considered, this would bring into play the effects of neutrino oscillations in a DM medium, as discussed above, but the study of this possibility is left for future work.

The possible interaction channels for neutrinos and UALPs are shown in Fig. 1 and the corresponding scattering amplitudes for the left and right diagrams are

(5)
(6)
Fig. 1
figure 1

Neutrino-UALP interaction channels considered

Full size image

The total squared average amplitude is

$$\begin{aligned} \left| \bar{\mathcal {M}} \right| ^2 = \frac{1}{2}\left( \left| \mathcal {M}_\mathrm{u}\right| ^2+ \mathcal {M}_\mathrm{u}\mathcal {M}_\mathrm{s}^*+ \mathcal {M}_\mathrm{s}\mathcal {M}_\mathrm{u}^* +\left| \mathcal {M}_\mathrm{s}\right| ^2 \right) , \end{aligned}$$
(7)

where for the different terms we obtain the following expressions valid if \(m_a \ll m_i\) and \(E_3\lesssim E_1\):

$$\begin{aligned}&\left| \mathcal {M}_\mathrm{u}\right| ^2= \frac{g_i^4}{E_3} \left[ {2 E_1} + \frac{8 m_i^4}{E_3m_a^2} - \frac{8m_i^2}{ m_a} \right] , \end{aligned}$$
(8)
$$\begin{aligned}&\left| \mathcal {M}_\mathrm{s}\right| ^2= \frac{g_i^4}{E_1} \left[ 2 {E_3} + 8 \frac{m_i^2}{m_a} + 8 \frac{m_i^4}{E_1 m_a^2} \right] , \end{aligned}$$
(9)
Fig. 2
figure 2

Minimum and maximum energies of the final neutrino as a function of the initial neutrino energy \(E_1\) for \(\nu _1,\ \nu _2,\) and \(\nu _3\)

Full size image

and

$$\begin{aligned}&\left| \mathcal {M}_\mathrm{s}\mathcal {M}_\mathrm{u}^*\right| = \left| \mathcal {M}_\mathrm{u}\mathcal {M}_\mathrm{s}^*\right| \nonumber \\&\quad =4{g_i^4} \left[ 2 \frac{(E_1 - E_3) }{m_a E_1 E_3} m_i^2 + 4 \frac{m_i^4}{E_1 E_3 m_a^2}-1\right] \end{aligned}$$
(10)
Fig. 3
figure 3

Total cross sections as a function of the neutrino energy for the three massive neutrinos

Full size image

For the neutrino masses, we consider two usual assumptions for their values masses of the neutrinos [37]: in the first one, the three masses are accommodated following a normal ordering (NO), i.e. \(m_1\ll m_2< m_3\) with

$$\begin{aligned} m_1= & {} 10^{-5}\,\,\mathrm{eV} \end{aligned}$$
(11)
$$\begin{aligned} m_2= & {} 8.7 \times 10^{-3}\,\,\mathrm{eV} \end{aligned}$$
(12)
$$\begin{aligned} m_3= & {} 5 \times 10^{-2}\,\,\mathrm{eV}, \end{aligned}$$
(13)

and the other option considered follows an inverted ordering (IO), \(m_3\ll m_2< m_1\), with

$$\begin{aligned} m_1= & {} 4.92\times 10^{-2}\,\,\mathrm{eV} \end{aligned}$$
(14)
$$\begin{aligned} m_2= & {} 5 \times 10^{-2}\,\,\mathrm{eV} \end{aligned}$$
(15)
$$\begin{aligned} m_3= & {} 10^{-5}\,\,\mathrm{eV}. \end{aligned}$$
(16)

We note that the actual value of lightest neutrino mass is unknown, and we set it at \(10^{-5}\) eV in order to produce our illustrative results.

The differential cross section in the laboratory (lab) frame is:

$$\begin{aligned} \frac{d\sigma _{\nu a}}{dt_M}= \frac{1}{64 \pi s_M}\frac{1}{|\mathbf {p}_{1,\mathrm cm}|^2}|\mathcal {M}|^2, \end{aligned}$$
(17)

where \(s_M=(p_1+p_2)^2\), and \( t_M= (p_3-p_1)^2 \) are the usual Mandelstam variables. The momenta of the initial and final neutrino are \(p_1\) and \(p_3\), respectively, while \(p_2\) and \(p_4\) refer to the initial and final UALP, respectively. The spatial momenta of the incident neutrino in the center of mass (CM) frame is \(\mathbf {p}_{1,\mathrm cm}\), and it satisfies that \(|\mathbf {p}_1|^2= E_{1,\mathrm cm}^2-m_i^2\), with

$$\begin{aligned} E_{1,\mathrm{cm}}= \frac{s_M+m_i^2-m_a^2}{2\sqrt{s_M}}. \end{aligned}$$

In the lab frame, in turn, \(t_M= -2m_a(E_1-E_3)\), and the differential cross section can be expressed as

$$\begin{aligned} \frac{d\sigma _{\nu a}}{dE_3}= \frac{m_a}{32 \pi s_M}\frac{1}{|\mathbf {p}_{1,\mathrm cm}|^2}|\mathcal {M}|^2. \end{aligned}$$
(18)

The total cross section is, therefore,

$$\begin{aligned} \sigma _{\nu a}=\int _{E_{3\mathrm min}}^{E_{3\mathrm max}} \frac{d\sigma _{\nu a}}{dE_3}, \end{aligned}$$
(19)

where the minimum and maximum energies of the final neutrino in the lab frame are

$$\begin{aligned} E_{3\mathrm{min}}= & {} \gamma \left( E_{3\,\mathrm{cm}}- \beta \sqrt{E_{3\mathrm{cm}}^2-m_i^2}\right) \end{aligned}$$
(20)
$$\begin{aligned} E_{3\mathrm{max}}= & {} \gamma \left( E_{3\,\mathrm{cm}}+ \beta \sqrt{E_{3\mathrm{cm}}^2-m_i^2}\right) , \end{aligned}$$
(21)

with \(\gamma =(E_1+m_a)/\sqrt{s_M}\), \(\beta =(1-\gamma ^2)^{-1/2}\), and \(E_{3\mathrm cm}=E_{1\mathrm cm}\). The latter energy values are plotted in Fig. 2 as a function of the initial neutrino energy \(E_1\). It can be seen that the energy loss per interaction is very small for \(E_1\lesssim 10^8\,\, \mathrm{GeV}\), except for the lightest neutrino, i.e., \(\nu _1\) in the NO case, and \(\nu _3\) in the IO case.

The total cross sections are shown in Fig. 3 for \(g=0.5 \,\,\mathrm{GeV}^{-1}\), and \(m_a=5\times 10^{-22}\,\,\mathrm{eV}\), where the coupling for each massive neutrino is given by \(g_i=2m_i g_{\nu a}\).

In the above calculations, we have assumed for simplicity that the UALP is at rest, as it is expected to constitute non-relativistic DM (e.g. [22] ). Even in the galactic halo, considering that its dispersion velocity can be \(\sim 100 \,\,\mathrm{km \ s^{-1}}\) [40], this leads to a Lorentz factor as low as \(\varGamma \simeq 1.0000005\) with respect to a reference frame fixed to the UALP. Since the strictly correct energies in the lab system should be multiplied by a factor \(\approx \varGamma \), a relative error of \(\epsilon _E\approx (\varGamma -1)\) is assumed, and this leads to an relative error in the kinematic variables \(s_\mathrm{M}\) and \(t_\mathrm{M}\) of \(\epsilon _{s}\approx 2(\varGamma -1)\). Hence, a simple estimate of error propagation in the calculation of \(d\sigma _{\nu a}/dE_{3}\) implies that the approximation of the UALP being at rest is accurate up to a negligible error much below \(1\%\).

Throughout this work, we shall adopt coupling values \(g_{i} \sim (0.1 - 1) \,\,\mathrm{GeV^{-1}}\), which are consistent with a cutoff for the effective field theory \(\varLambda \sim (1-12)\,\, \mathrm{GeV}\). Still, the CM energies of the \(\nu _i-a\) interactions, \(\sqrt{s_M}=( m_i^2+m_a^2+ 2m_a E_{1})^{0.5}\), are well below this cutoff for a UALP mass \(m_a=5\times 10^{-22}\,\,\mathrm{eV}\) and neutrino energies \(E_1< 10^8\,\, \mathrm{GeV}\). Therefore, the validity of the effective operator approach is not compromised in the present context.

3 Neutrino propagation including interaction effects

In this section we describe the propagation of a diffuse flux of high energy neutrinos of astrophysical origin, considering separately the extragalactic propagation and then the propagation through the DM halo of our galaxy.

3.1 Injection by astrophysical sources

The neutrino injection depends on the proton acceleration mechanism operating at the astrophysical sources, and in particular, we are interested in the corresponding initial flavor ratios \(f_{\alpha ,\mathrm s}(E)\), i.e., at the beginning of the propagation stage. If the production channel is neutron decay (\(n\rightarrow p \ \bar{\nu }_e e^-\)), then only the electron flavor is generated, and hence the flavor composition at the sources is simply \((f_{e,\mathrm s}:f_{\mu ,\mathrm s}:f_{\tau ,\mathrm s})=(1:0:0)\). However, neutrino production after pion decays is more promisory as pions can be efficiently created by pp and \(p\gamma \) interactions at the sources.Footnote 1

We then consider two additional source benchmarks: in the first one, both pions and muons decay unaffected by synchrotron losses leading to the usually expected flavor composition at the sources \((f_{e,\mathrm s}:f_{\mu ,\mathrm s}:f_{\tau ,\mathrm s})=(1:2:0)\), and we reffer to this case as pion decay sources. In the other typical case, muons are completely depleted due to synchrotron losses and we have a proportion of flavors \((f_{e,\mathrm s}:f_{\mu ,\mathrm s}:f_{\tau ,\mathrm s})=(0:1:0)\) at these highly magnetized muon-damped sources. There are also other intermediate possibilities for which muons could be cooled significantly only above a certain energy, and in these cases the flavor at the sources would be energy dependent [27, 42]. Still, similarly to the assumed in previous works [43,44,45,46,47], here we consider the mentioned three benchmark cases of initial flavor compositions, as our main goal is to focus on the effects of the possible neutrino interactions with UALP DM. To achieve this, we adopt for the neutrino injection a canonical power-law dependence with the neutrino energy times an exponential cutoff at \(E_\mathrm{c}= 10^{7}\,\,\mathrm{GeV}\):

$$\begin{aligned} Q_{\nu _i}(E,z)= K_\nu \frac{\psi _\mathrm{SFR}(z)}{\psi _\mathrm{SFR}(z=1)} f_{i,\mathrm s} E^{-\alpha }\exp \left( -\frac{E}{E_\mathrm{c}}\right) . \end{aligned}$$
(22)

Here, \(K_\nu \) is a normalization constant, and the factors \(f_{i,\mathrm s}\) depend on the flavor composition at the sources since the neutrino injection corresponding to the flavor \(\alpha \) is given by

$$\begin{aligned} Q_{\nu _\alpha }(E,z)=\sum ^3_{i=1} \left| U_{\alpha i}\right| ^2 f_{i,\mathrm s}Q_{\nu _i}(E,z). \end{aligned}$$
(23)

The evolution with the redshift z is assumed to follow the star formation rate (SFR) as given by [23],

$$\begin{aligned} \psi _\mathrm{SFR}(z)= 0.015\frac{(1+z)^{2.7}}{1+ \left[ (1+z)/2.9\right] ^{5.6}} \ M_{\odot }\,\,\mathrm{year^{-1}\,Mpc^{-3} }.\nonumber \\ \end{aligned}$$
(24)

Here we remark that we assume that the sources are also isotropically distributed in the sky, and the approach we apply in this work is adequate to capture the possible modification the total diffuse neutrino flux by accounting for the energy loss undergone by neutrinos if they interact with UALP DM. A study of these effects in the case of individual point sources should also account for the angular deflection generated by the scatterings, as it was done in Ref. [12] for the flare neutrino event IceCube-170922A associated to the blazar TXS 0506+056 [48].

In the present context, and for illustration of the effects of the described interactions on a diffuse neutrino flux, we choose to normalize the injection as

$$\begin{aligned}&\frac{c}{4\pi }\int _0^{5}\frac{dz}{(1+z)^2} \frac{dt}{dz} Q_{\nu _\mu }\left( E(1+z),z\right) \nonumber \\&\quad \times \left| _{E=10^5\mathrm{GeV}}= \varPhi ^\mathrm{IceCube}_{\nu _\mu +\bar{\nu }_\mu }(E)\right| _{E=10^5\mathrm{GeV}}, \end{aligned}$$
(25)

so that the \(\nu _\mu +\bar{\nu }_\mu \) flux obtained in the absence of \(\nu -a\) interactions matches the best fit flux obtained by IceCube evaluated at \(10^5\)GeV [51],

$$\begin{aligned}\varPhi ^\mathrm{IceCube}_{\nu _\mu +\bar{\nu }_\mu }(E)=\varPhi _0\ \left( \frac{E}{10^5 \mathrm{GeV}}\right) ^{-2.28}, \end{aligned}$$

where \(\varPhi _0=1.44\times 10^{-18}\,\,\mathrm{GeV^{-1}\,s^{-1}\,sr^{-1}\,cm^{-2}}\), and we accordingly assume \(\alpha =2.28\) for the index of injection.

3.2 Extragalactic propagation

We treat the extragalactic propagation of astrophysical neutrinos making use of the following transport equation with continuous losses for the comoving density of each massive neutrino [39, 49]

$$\begin{aligned}&\frac{\partial N_{\nu _i}(E,t)}{\partial t}= Q_{\nu _i}(E,t)-3H(t)N_{\nu _i}(E,t)\nonumber \\&\quad +\frac{\partial \left[ H(t)E\ N_{\nu _i}(E,t)+N_{\nu _i}(E,t)b_{\nu a}(E,t)\right] }{\partial E}, \end{aligned}$$
(26)

which accounts for the effect of expansion of the universe through the terms proportional to

$$\begin{aligned}H(t)= H_0\left[ \varOmega _\mathrm{m}(1+z)^3+\varOmega _\varLambda \right] ^{-1},\end{aligned}$$

and the interactions with DM UALPs are characterized by the corresponding continuous energy loss rate [50]:

$$\begin{aligned} b_{\nu a}(E,z)= n_\mathrm{dm}c\int _{E_{3\mathrm min}}^{E_{3\mathrm max}} dE_3 (E-E_3) \frac{d\sigma _{\nu a}(E,E_3)}{dE}. \end{aligned}$$
(27)

Here, the density of extragalactic DM is given by

$$\begin{aligned} n_\mathrm{dm}(z)= \left( \frac{3H_0^2}{8\pi G}\right) \frac{\varOmega _\mathrm{m}}{m_a}(1+z)^3, \end{aligned}$$
(28)

with \(H_0\simeq 70\,\,\mathrm{km \, s^{-1}}\) and \(\varOmega _\mathrm{m}\simeq 0.3\). For instance, the extragalactic DM column \(X_\mathrm{dm}(z)=\int _0^z n_\mathrm{dm}(z') dz'\) for sources at \(z\approx 1\) is \(X_\mathrm{dm}(1) \approx 6\times 10^{52}\,\,\mathrm{cm^{-2}} \), and for sources at a redshift \(z\approx 5\) it is \(X_\mathrm{dm}(5) \approx 6\times 10^{53}\,\,\mathrm{cm^{-2}} \). In Eq. (27), E represents the initial neutrino energy and we integrate over the final neutrino energy \(E_3\) in a range that can be seen in Fig. 2 and is determined by \(E_{3\mathrm min}\) and \(E_{3\mathrm max}\) as defined above above in Eqs. (20,21). We note that in the NO case for \(\nu _2\) and \(\nu _3\), the energy loss per interaction is very small for \(E_1\lesssim 10^8\)GeV, and hence the continuous loss approximation is justified. However, this is not the case for \(\nu _1\), given that the range of initial neutrino energies \(E_1\) for a fixed final energy \(E_3\) is much broader, and additionally, the differential cross section \(\frac{d\sigma _{\nu a}}{dE_3}\) is actually peaked at \(E_1\ll E_3\). Therefore, the present approach is accurate to obtain the neutrino fluxes for a coupling

$$\begin{aligned} g_{\nu a}\lesssim 1.5 \mathrm{GeV}^{-1} \end{aligned}$$
(29)

such that no significant \(\nu _1-a\) interactions take place, i.e.,

$$\begin{aligned} \sigma _{\nu a}< & {} \frac{1}{X_\mathrm{dm}(z=5)}\approx 1.6\times 10^{-54}\,\,\mathrm{cm^{2}}. \end{aligned}$$
(30)

A similar reasoning can be applied to the IO case, where the heavier neutrinos are \(\nu _1\) and \(\nu _2\), and the lightest one will be affected only for the mentioned range of large couplings. For higher values of \(g_{\nu a}\), the only surviving flux would be the one of \(\nu _1\) in the NO case and the one of \(\nu _3\) in the IO case, so that the corresponding flavor composition would be clearly determined as

$$\begin{aligned} (f_e:f_\mu :f_\tau )_\mathrm{NO}= & {} |(U_{e1}|^2:|U_{\mu 1}|^2:|U_{\tau 1}|^2) \end{aligned}$$
(31)
$$\begin{aligned} (f_e:f_\mu :f_\tau )_\mathrm{IO}= & {} |(U_{e3}|^2:|U_{\mu 3}|^2:|U_{\tau 3}|^2), \end{aligned}$$
(32)

if absorption is not complete, but affects only the two heavier neutrinos. More on this issue is discussed below.

In terms of redshift, considering that \(dt/dz=-H(z)(1+z)\), we can write the transport equation as

$$\begin{aligned}&\frac{\partial N_{\nu _i}(E,z)}{\partial z}=-\frac{Q_{\nu _i}(E,z)}{H(z)(1+z)} \nonumber \\&\quad +\left[ \frac{1}{1+z}\right] \left[ {2}- \frac{1}{H(z)}\frac{\partial b_{\nu a}(E,z)}{\partial E}\right] N_{\nu _i}(E,z) \nonumber \\&\quad -\left[ \frac{E}{1+z}+\frac{b_{\nu a}}{H(z)(1+z)}\right] \frac{\partial N_{\nu _i}(E,z)}{\partial E}, \end{aligned}$$
(33)

and the neutrino flux after the propagation is found as

$$\begin{aligned} {\varPhi _{\nu _\alpha }}(E)=\frac{c}{4\pi }\sum _{i=1}^3f_{i,}{\left| U_{\alpha \,i}\right| ^2N_{\nu _i} (E)}. \end{aligned}$$
(34)

We solve Eq. (33) using the method of the characteristics, i.e., with the characteristic curves satisfying

$$\begin{aligned} \frac{dE}{dz}= \frac{E}{1+z}+\frac{b_{\nu a}}{H(z)(1+z)}, \end{aligned}$$
(35)

we obtain the solutions

$$\begin{aligned}&N_{\nu _i}(z,E)= \int _z^{z_\mathrm{max}}dz'\frac{Q_{\nu _i}(z',E')}{H(z')(1+z')}\nonumber \\&\quad \times \exp {\left[ -\int _z^{z'} \frac{dz''}{1+z''}\left( 2-\frac{1}{H(z'')}\frac{\partial b_{\nu a}}{\partial E}\right) \right] }, \end{aligned}$$
(36)

where we assume \(z_\mathrm{max}=5\). The results obtained for \(z=0\) are then used as input to describe the propagation through our galaxy, as discussed below.

Fig. 4
figure 4

Side and top schematic views of the galactic DM halo, indicating the galactic coordinates l and b

Full size image

3.3 Propagation through the DM halo

After propagation outside our galaxy, a diffuse flux of neutrinos arriving to the Milky Way DM halo has to traverse it at different directions, facing different DM column depths \(X_\mathrm{dm,h}(l,b)\) depending on the galactic longitudes l and latitudes b (see Fig. 4). The DM density profile adopted is a generalized spherically symmetric Navarro, Frenk, and White (NFW) one [52, 53],

$$\begin{aligned} n_\mathrm{dm,h}(r)= \frac{\rho _\mathrm{s}}{m_a}\left( \frac{r}{r_\mathrm{s}}\right) ^{-\gamma _{0}}\left( 1+\frac{r}{r_\mathrm{s}}\right) ^{\gamma _0-3}, \end{aligned}$$
(37)

where \(\gamma _0=1.2\), r is the distance to the galactic center, \(r_\mathrm{s}=20\,\mathrm{kpc}\), and the density there (\(\rho _\mathrm{s}\)) is obtained by assuming a local density \(\rho _0=0.4 \, \mathrm{GeV \, cm^{-3}}\) at the position of our solar system \(r_0=8.127 \,\mathrm{kpc}\). For instance, we obtain \(X_\mathrm{dm,h}(l=0,b=10^\circ )= 2\times 10^{53}\,\,\mathrm{cm^{-2}}\) and \(X_\mathrm{dm,h}(l=0,b=90^\circ )= 3\times 10^{52}\,\,\mathrm{cm^{-2}}\). The transport equation in the absence of injection and adiabatic losses can be written as

$$\begin{aligned} \frac{\partial N_{\nu _i}(E,l,b)}{\partial X}= \frac{1}{ n_\mathrm{dm,h} c }\frac{\partial \left[ b_{\nu a,\mathrm h}(E) N_{\nu i}(E,l,b)\right] }{\partial E}, \end{aligned}$$
(38)

where \(b_{\nu a,\mathrm h}(E)\) is analogous to Eq. (27) but with the DM density corresponding to the galactic halo.

Using again the method of the characteristics, we find that the solution of Eq. (38) can be expressed as

$$\begin{aligned} N_{\nu _i}(E,l,b)= N_{\nu _i}(E,z=0) \left[ \frac{b_{\nu a}\left( E'({X'=0})\right) }{b_{\nu a}(E)}\right] , \end{aligned}$$
(39)

where \(E'\) is the characteristic energy corresponding to a depth \(X'\) such that the final energy for a depth \(X_\mathrm{dm,h}(l,b)\) is E.

4 Results

Fig. 5
figure 5

Diffuse fluxes of \(\nu _\mu +\bar{\nu }_\mu \) assuming a normal ordering for the neutrino masses, for \(g_{\nu a}=0.5\,\mathrm{GeV^{-1}}\) and \(g_{\nu a}=1\,\mathrm{GeV^{-1}}\) in the left and right panels, respectively. Top panels correspond to pion decay sources, middle panels to muon damped sources, and bottom panels to neutron decay sources. We show the fluxes corresponding to the case of no interactions in dashed lines, and for reference we include the 10-year CL 68% IceCube data corresponding to muon tracks. The flux corresponding to interactions with extragalactic DM only is indicated with a black solid line in each panel, and the fluxes for the arrival directions given by \(l=0^\circ \) and \(b=\left\{ 90^\circ , 45^\circ ,0^\circ \right\} \) are marked with curves in dark gray, gray, and light gray, respectively

Full size image
Fig. 6
figure 6

Diffuse fluxes of \(\nu _e+\bar{\nu }_e\) assuming an inverted order for the neutrino masses, for \(g_{\nu a}=0.5\,\mathrm{GeV}^{-1}\) and \(g_{\nu a}=1\,\mathrm{GeV^{-1}}\) in the left and right panels, respectively. Top panels correspond to pion decay sources, middle panels to muon damped sources, and bottom panels to neutron decay sources. The flux corresponding to interactions with extragalactic DM only is indicated with a black solid line in each panel, and the fluxes for the arrival directions given by \(l=0^\circ \) and \(b=\left\{ 90^\circ , 45^\circ ,0^\circ \right\} \) are marked with curves in dark gray, gray, and light gray, respectively

Full size image

In this section, we present the obtained results for the neutrino fluxes after extragalactic and galactic propagation including the effects of \(\nu -a\) scatterings. As discussed above, we consider three different benchmark possibilities for the flavor composition at the sources.

The proportion of the different massive neutrinos emitted at the sources \((f_{1,\mathrm s}:f_{2, \mathrm s}:f_{3,\mathrm s})\) in the different cases assumed is such that

$$\begin{aligned} f_{i,\mathrm s}\propto \left\{ \begin{array}{ll} \left( |U_{ei}|^2+2|U_{\mu i}|^2\right) &{} \quad \text{ for } \pi \text{ and } \mu \text{ decays } \\ \left( |U_{\mu i}|^2\right) &{} \quad \text{ for } \mu \text{ damped } \text{ case } \\ \left( |U_{ei}|^2\right) &{} \quad \text{ for } n \text{ decays } \end{array} \right. . \end{aligned}$$

The neutrino propagation including the effect of interactions is described by the solution of the trasport equations for each massive neutrino as mentioned above, and the flux of neutrinos of the different flavors is computed by superposition as in Eq. (34). In Fig. 5, we show the obtained fluxes of \(\nu _\mu + \bar{\nu }_\mu \) in the NO case including the effects of interactions for \(g_{\nu a}=0.5 \mathrm{GeV^{-1}}\) and \(g_{\nu a}=1 \mathrm{GeV^{-1}}\). We show the fluxes after extragalactic propagation only as well as the arriving in directions given by \(l=0^\circ \) and \(b=(20^\circ , 45^\circ , 90^\circ )\), and we compare these results with the flux corresponding to no interactions and with the flux of \(\nu _\mu + \bar{\nu }_\mu \) measured by IceCube at CL 68% based on 10 years of muon track data [51]. It can be seen that the interaction effects are more noticeable if the arrival directions are closer to the galactic center (GC) since the DM column along the neutrino path is greater. The effect of increasing the coupling strength can be appreciated by comparing the results of the left panels to those of the right panels of Fig. 5: a more significant flux attenuation takes place as the coupling is increased, as can be expected due to a larger cross section. It can also be seen that the fluxes are most affected in the case of muon damped sources and least affected for neutron decay sources. This can be understood taking into account that the cross section is higher for neutrino \(\nu _3\), followed by that for \(\nu _2\) and the weakest one is for \(\nu _1\), as can be seen in Fig. 3. Therefore, since in the case of muon-damped sources only one muon neutrino is produced, the fact that \(|U_{\mu 3}|\approx 0.55\) implies that a significant portion of the emitted flux would be more affected by interactions. Likewise, in the neutron-decay case, only the electron flavor is produced and since \(|U_{e 3}|\approx 0.022\), the net attenuation is weaker and due basically to the depletion corresponding to the neutrino \(\nu _2\), while the \(\nu _1\) remains unaffected for the coupling values considered.

In Fig. 6, we show the diffuse fluxes of the electron flavor neutrinos in the NO case, since in this flavor the fluxes are more affected. This is because in the IO case, the neutrino \(\nu _3\) is the least affected one and its contribution to the electron flavor is very small, given that \(|U_{e3}|\approx 0.022\) as in the NO case. In order to appreciate the contribution of the three flavors, we show the energy dependence of the flavor ratios on the Earth in Fig. 7 for the NO case and in Fig. 8 for the IO case. In these plots we include the averaged results over two regions in the sky: one for \((-90^\circ<l<90^\circ )\), i.e., centered around the GC, and the other for \((90^\circ< l <270^\circ )\), i.e., centered around the galactic anti-center (GAC). Clearly, neutrinos corresponding to the former hemisphere face a higher column of galactic DM than those arriving corresponding to the GAC hemisphere, and attenuation is more significant. By comparing the results of Fig. 7 with the ones of Fig. 8, it can be seen that the effects of interactions lead to very different flavor ratios if neutrino masses follow a NO or if they follow an IO. In the former case and in relation to the flavor composition expected without interactions, the electron flavor is increased and becomes dominant, followed by the tau flavor which is decreased, and the muon flavor contribution is significantly reduced. Quite the opposite is expected in the IO case if interactions are important: the muon flavor becomes dominant, followed by the tau flavor still with reduced relevance, and the electron flavor contribution becomes very small.

Fig. 7
figure 7

Energy dependent flavor ratios of neutrinos reaching the Earth in the NO case for the electron, muon, and tau flavors in the left, center, and right panels, respectively. Thick solid lines correspond to the half of the sky centered around the galactic center (GC), while thin lines correspond to the other half of the sky containing the galactic anti-center (GAC). The results for no interactions are marked using dashed lines. Top panels correspond to pion decay sources, middle panels to muon damped sources, and bottom panels to neutron decay sources

Full size image
Fig. 8
figure 8

Energy dependent flavor ratios of neutrinos reaching the Earth in the IO case for the electron, muon, and tau flavors in the left, center, and right panels, respectively. Thick solid lines correspond to the half of the sky centered around the galactic center (GC), while thin lines correspond to the other half of the sky containing the galactic anti-center (GAC). The results for no interactions are marked using dashed lines. Top panels correspond to pion decay sources, middle panels to muon damped sources, and bottom panels to neutron decay sources

Full size image

We can integrate on the neutrino energy the differential fluxes in order to obtain the flavor ratios of the integrated fluxes above \(E_\mathrm{min}\) as:

$$\begin{aligned} F_{\alpha }=\frac{\int _{E_\mathrm{min}}^\infty dE_\nu \varPhi _{\nu _\alpha }(E_\nu )}{\int _{E_\mathrm{min}}^\infty dE_\nu \left[ \varPhi _{\nu _e}(E_\nu )+ \varPhi _{\nu _\mu }(E_\nu )+ \varPhi _{\nu _\tau }(E_\nu )\right] }. \end{aligned}$$
(40)

The results for \(E_\mathrm{min}=5\times 10^5\mathrm{GeV}\) are shown in Fig. 9 as a function of the galactic latitude b for the NO case and in Fig. 10 for the IO case. In this figures, we show the flavor ratios of the integrated fluxes for the three cases of sources assumed and with different values of the coupling, \(g_{\nu a}=0.5\,\,\mathrm{GeV}^{-1}\) and \(g_{\nu a}= 0.75\,\mathrm{GeV}^{-1}\). Also, as mentioned in Sect. 2, for high values of the coupling, only the \(\nu _{1(3)}\) would survive in the NO(IO) case, and can be completely unaffected due to a much lower cross section. In these cases, the flavor ratios could also be homogeneous across the full sky, but fixed at the particular values

$$\begin{aligned} (F_e:F_\mu :F_\tau )_\mathrm{NO}= & {} \left( |U_{e1}|^2:|U_{\mu 1}|^2:|U_{\tau 1}|^2\right) \end{aligned}$$
(41)
$$\begin{aligned}\approx & {} (0.68: \, 0.075: \, 0.24)\nonumber \\ (F_e:F_\mu :F_\tau )_\mathrm{IO}= & {} \left( |U_{e3}|^2:|U_{\mu 3}|^2:|U_{\tau 3}|^2\right) \nonumber \\\approx & {} (0.022: \, 0.56: \, 0.41). \end{aligned}$$
(42)
Fig. 9
figure 9

Flavor ratios of the integrated fluxes in the NO case as a function of the galactic latitude b corresponding to the arrival direction, for a fixed galactic longitude \(l=0^\circ \). The electron, muon, and tau flavor ratios are shown in the left, middle, and right panels, respectively. Top panels correspond to pion decay sources, middle panels to muon damped sources, and bottom panels to neutron decay sources

Full size image
Fig. 10
figure 10

Flavor ratios of the integrated fluxes in the IO case as a function of the galactic latitude b corresponding to the arrival direction, for a fixed galactic longitude \(l=0^\circ \). The electron, muon, and tau flavor ratios are shown in the left, middle, and right panels, respectively. Top panels correspond to pion decay sources, middle panels to muon damped sources, and bottom panels to neutron decay sources

Full size image

In fact, this can be illustrated with a ternary plot as the one in Fig. 11, where each position represents a unique combination of the three flavor ratios obtained by integration of the neutrino fluxes over the energy and also over all incoming directions. In this plot, we mark with a white triangle the expected values in the cases of no interactions corresponding to pion decay (top panels), muon damped (middle panels), and neutron decay sources (bottom panels). If the coupling is set as \(g_{\nu a}=0.5\,\,\mathrm{GeV}^{-1}\) then the results shift to the indicated by the gray squares, to the dark gray ones if we set \(g_{\nu a}=0.75\,\,\mathrm{GeV}^{-1}\), and to the black ones for \(g_{\nu a}=1\,\,\mathrm{GeV}^{-1}\). We also include the \(68\%\) CL region obtained with the sample of high energy starting events (HESE) by IceCube [51], and the projected \(68\%\) CL constrained ragion with 8 yr of IceCube if the composition were that for pion decay sources, i.e., \((f_{e,\mathrm s}{:}f_{\mu ,\mathrm s}{:}f_{\tau ,\mathrm s})=(2{:}1{:}0)\) (see Ref. [27]).

Fig. 11
figure 11

Ternary plots of the flavor ratios for neutrinos above \(E_\mathrm{min}=5\times 10^5\mathrm{GeV}\) in the NO case (left plots) and in the IO case (right plots) under the assumptions of pion decay sources (upper panels), muon damped sources (middle panels), and neutron decay sources (bottom panels). The outcomes expected in the absence of neutrino interactions are marked by white triangles, while the results with interactions enabled are shown in gray squares for \(g_{\nu a}=0.5\,\,\mathrm{GeV}^{-1}\), in dark gray squares for \(g_{\nu a}=0.75\,\,\mathrm{GeV}^{-1}\), and in black for \(g_{\nu a}=1\,\,\mathrm{GeV}^{-1}\). For reference, we show with black-dashed lines the contour corresponding to the \(68\%\) C.L. obtained by IceCube with 7.5 years of accumulated data of the HESE sample [51], and we show in gray the projected region expected to be obtained with IceCube-gen2 if the sources were dominated by pion decay

Full size image

5 Discussion

In this work, we have analyzed the impact of possible interactions of astrophysical neutrinos with UALP DM particles a, with \(m_a\approx 5\times 10^{-22}\,\,\mathrm{eV}\). Taking into account a diagonal coupling with the same value for the three flavors (\(g_{\nu _\alpha a}=g_{\nu a}\)), we treated the extragalactic and galactic propagation using transport equations for the density of each massive neutrino. Considering as benchmarks three usually adopted initial flavor compositions for pion decay, muon-damped , and neutron decay sources, we found that interactions can cause an important change in the flavor composition to be observed on Earth with neutrino telescopes, and this effect is also sensible to the the mass ordering of the massive neutrinos. As the coupling is increased gradually from \(g_{\nu a}=0.5\,\,\mathrm{GeV}^{-1}\) to \(g_{\nu a}=1\,\,\mathrm{GeV}^{-1}\), the flavor composition gradually departs from the result corresponding to no interactions to \((F_e:F_\mu :F_\tau )_\mathrm{NO}= (|U_{e1}|^2:|U_{\mu 1}|^2:|U_{\tau 1}|^2)\) in the NO case, and to \((F_e:F_\mu :F_\tau )_\mathrm{IO}= (|U_{e3}|^2:|U_{\mu 3}|^2:|U_{\tau 3}|^2)\) in the IO case. We also found that for somewhat smaller couplings \(0.5\, \mathrm{GeV}^{-1}\lesssim g_{\nu a}\lesssim 1\, \mathrm{GeV}^{-1}\), neutrinos are affected differently at different directions in the sky because the column of galactic DM will be higher for neutrino path passing closer to the galactic center. Therefore, if the measured neutrino flavor composition is homogeneously distributed across the whole sky, this will be evidence against neutrino interactions with UALP DM in the mentioned range of the coupling \(g_{\nu a}\). Considering the current data from IceCube regarding the flavor composition [51], the neutrino interactions discussed in this work with a coupling \(g_{\nu a} > rsim 1\, \mathrm{GeV}^{-1}\) appear to be disfavored, as can be seen in the ternary plots of Fig. 11. Still, more statistics is necessary to obtain a higher precision in the experimental determination of the flavor ratios and this will be achieved with longer time exposure and larger detectors such as IceCube-gen2 [27].

These conclusions are in agreement with the constraint estimated in Ref. [9] based on the detection of the neutrino event IceCube-170922A, leading to \(g_{\nu a} > rsim 0.8\mathrm{GeV}^{-1}\), and in particular our analysis shows that a lower coupling \(g_{\nu a}=0.5\mathrm{GeV}^{-1}\) interactions leave their signature by causing a change in the flavor ratios for energies greater than \(\sim 10^5\mathrm{GeV}\), with respect to the result corresponding to no interactions, as shown in the left panels of Figs. 7 and 8.

Given that the energy dependence of the flavor ratios is expected to be experimentally determined by future measurements [27, 28], the analysis here presented can be useful to probe a diagonal coupling with a value \(g_{\nu a}\approx 0.5\mathrm{GeV}^{-1}\) for all the flavors. This is below the ranges that, as mentioned in Ref. [9], can be excluded by avoiding a too fast cooling in supernovae (\(g_{\nu a}\lesssim 10^4\mathrm{GeV^{-1}}\)) [54] and also by allowing the free streaming of neutrinos in the early universe before photon decoupling (\(g_{\nu a}\lesssim 10^3\mathrm{GeV^{-1}}\)) [35].Footnote 2 In the cases of different couplings for each flavor or even non-diagonal couplings, observations with upcoming underground neutrino detectors such as JUNO [55] and DUNE [24] are expected to explore values as \(g_{\nu _\alpha a}\sim 0.1-0.01\,\,\mathrm{GeV}^{-1}\) searching for oscilation effects [9].

We leave for a future work the treatment of higher energy neutrino fluxes such as the so-called comogenic neutrinos produced by cosmic ray interactions with the cosmic microwave background. Given that the energy loss per interaction becomes higher at higher energies, the interactions of these neutrinos is therefore expected to alter the flavor composition for lower \(g_{\nu a}\) values than the explored in the present work.