Explanation of IceCube spectrum with $ν\rightarrow 3 ν$ neutrino splitting in a $ν$2HDM model
Subhendra Mohanty Soumya Sadhukhan

Towards the
Raychaudhuri Equation Beyond General Relativity

    Bret Underwood
    December 30, 2023
============================================================






§ INTRODUCTION




After seven years of IceCube measurement, a clear 6σ excess of events is observed for energies above 60 TeV which cannot be 
explained by the atmospheric neutrinos <cit.>. 
The first choice for explaining the high energetic neutrino events were different astrophyscal sources <cit.>. 
The source and spectrum of neutrinos observed at IceCube remains a puzzle
as the events do not point back to any clear identifiable Active Galactic Nuclei (AGN) or Gamma Ray Bursts (GRB) and the entire spectrum of events from 60 TeV to 10 PeV cannot be explained by a single power law of neutrino flux like Φ_ν= Φ_0 E_ν^α. 
In GRB models  <cit.>, charged pions are produced when high energetic protons interact with gamma background as 

p + γ→Δ^+→ n + π^+, followed by further decay of pions as π^+→μ^+ν_μ and π^+→ e^+ ν_eν_μ, 
producing flux of high energetic muon and electron neutrinos. IceCube is sensitive to those neutrinos whose energies are above 0.1 TeV. 
In general the energy dependence of the incoming neutrino flux from the GRBs 
at the source is given by the Waxman Bahcall spectrum <cit.> Φ_ν∝  E^-2. 
The astrophysical neutrino flux modeled in this way does not give a good fit to the IceCube event distribution at all energy bins upto ∼ 2 PeV. 
With a fit at lower energy (∼ 100 TeV) bins, this model predicts an excess of neutrino events at higher energy bins, but IceCube has not observed that effect till now. 
Even if the decrease of flux amplitude can fit neutrino event observation at some bin, there remains huge mismatch with predicted events at other energies. 
If a steeper neutrino energy spectrum is taken, IceCube event distribution in the sub-PeV bins can be explained by fixing a proper neutrino flux amplitude. 
The steepness of the spectrum will result in deficiency of predicted events at energies 1-3 PeV. Invisible and visible decays of astrophysical 
neutrinos is explored in this work as a better explanation of 6 years of IceCube event distribution data. Various models for explaining the observed features in the IceCube event spectrum include neutrinos from Pev dark matter decay or annihilation to explain the 1 PeV excess (<cit.>). 
The excess of PeV events can also come from the resonant production of leptoquarks  <cit.>. 

On the other hand, there are neutrino depletion models which try to explain the non-observation of the Glashow resonance  <cit.> i.e. absence of the resonant peak due to a real W^- production in the process ν̅_e e^- → W^-. 
The decay of the real W is expected to give hadron and lepton shower or lepton track events  <cit.>. Depletion of high-energy neutrinos can occur via oscillation to sterile neutrinos in pseudo-Dirac neutrinos <cit.> and for visible decay  <cit.>. Exotic scenarios have also been invoked to explain a cutoff at the Glashow resonance energies such as Lorentz violation  <cit.> and CPT violation <cit.>. 

Neutrino splitting ν→ 3 ν can change the neutrino spectrum at IceCube and it can occur in CPT and Lorentz violation models <cit.>. 
In this paper we consider a Lorentz invariant BSM model where the neutrino splitting ν→ 3 ν occurs through a one-loop diagram with a light mediator.
One BSM case which can give rise to the neutrino decay is neutrinophilic two Higgs doublet model <cit.> 
(ν2HDM) which was originally proposed to explain the nonzero but tiny neutrino mass. 
This kind of neutrinophilic models of neutrino mass generation was initially proposed and discussed in the References <cit.>.
In ν2HDM, neutrino masses are generated by introducing a second scalar doublet with a vacuum expectation value of the eV scale. 
The astrophysical neutrino decay rate should be large enough to have a decay before neutrinos reaching the IceCube detector. 
Larger decay rates are ensured in the ν2HDM through the presence of a eV scale neutral scalar. 
Presence of an ultra light scalar results in a stringent bound from the oblique parameters. 
This can be resolved by expanding the model with vectorlike leptons, which themselves have interesting phenomenological properties. 

Neutrino decay to visible lower energy neutrinos change the spectrum of IceCube neutrinos and give a better fit of the IceCube event spectrum in the entire range 60 TeV -10 PeV
compared to a single power law neutrino flux. 

We note that in the neutrino splitting scenario, it protects the information of astrophysical source direction, as the daughter neutrinos are collinear being very energetic ones, which is consistent with the recent observation that one single source can be traced back from the observed IceCube spectrum.


In section <ref> we have briefly explained the ν2HDM model along with the explanation of emergence of tiny neutrino mass there. 
In the next section <ref>, it is shown how active SM neutrinos can split to three active neutrinos through the box loop diagram. We have also 
shown the presence of an eV scale scalar propagator increases the decay rate and therefore ensure a finite decay lifetime to have any astrophysical effect. 
In the next section <ref>, we show the flavor ratios of the daughter neutrinos after splitting, taking into account different astrophysical 
sources with different initial flavor ratios.
In section <ref>, we discuss how the presence of an ultra light scalar in ν2HDM can have large negative contribution on S, T parameters to push 
them to the values which are not experimentally viable. In section <ref> different flavor, collider and theoretical constraints are discussed along with 
their bounds on ν2HDM parameter space. In section <ref>, we do the analysis of the effects of both the neutrino splitting (ν→ 3 ν) and invisible decays  on the 
astrophysical neutrino propagation and their event distribution at IceCube. In one subsection <ref> we concentrate on visible neutrino decay and show how its presence 
fit the IceCube event distribution along with an comparison with the case where there is no neutrino decay. In another subsection <ref>, similar studies are done for 
neutrino invisible decay. Finally we summarize and conclude in section <ref>.




 

§ THE MODEL


The beyond the SM ν2HDM theory <cit.> discussed in this work, is based on the symmetry group SU(3)_c × SU(2)_L × U((1)_Y × Z_2. 
In addition to the usual SM fermions, we have one EW singlet right-handed (RH) neutrino, ν_R, for all flavors of SM leptons.  
The model has two Higgs doublets, Φ_1 and Φ_2.  
All the SM fermions and the Higgs doublet Φ_1, are even under the discrete symmetry, Z_2, while
the RH neutrino and the Higgs doublet Φ_2 are odd under Z_2.
This leads to Yukawa interaction of all the SM fermions except the left-handed neutrinos, through Φ_1 only. 
The SM left-handed neutrinos, together with the right-handed neutrino added here, couple to the Higgs doublet Φ_2. 
The discrete symmetry Z_2 is broken by a vev of Φ_2, and we take v_2 = ⟨Φ_2⟩∼ 0.1 eV. 
Therefore, the origin of the neutrino mass in ν2HDM is due to a spontaneous breaking of the discrete symmetry Z_2.
Through their Yukawa interactions with the Higgs field Φ_2, the neutrinos acquire orders of magnitude smaller masses 
than the SM ones even with order one Yukawa couplings, due to the smallness of v_2.

For simplicity, we consider the Higgs sector to be CP invariant here. 
The most general Higgs potential consistent with the SM × Z_2 symmetry
is <cit.>


    V =
    -μ^2_1 Φ_1^†Φ_1-μ^2_2 Φ_2^†Φ_2+λ_1(Φ_1^†Φ_1)^2+λ_2(Φ_2^†Φ_2)^2+ 
    λ_3(Φ_1^†Φ_1)(Φ_2^†Φ_2)-
    λ_4|Φ_1^†Φ_2|^2
    -1/2λ_5[(Φ_1^†Φ_2)^2+(Φ_2^†Φ_1)^2].

In the ν2HDM scalar potential of Eq. <ref>, the parameters
m_12^2 and λ_6,7 are put to zero due to the imposition of Z_2 symmetry.
The physical Higgs fields are charged fields H^±, two neutral CP even scalar fields h and H, and a neutral CP odd field A. 
In the unitary gauge, the two doublets can be written


    Φ_1 =  1/√(2)([       √(2) (v_2/v)H^+; h_0 + i (v_2/v)A +v_1;                       ]),


    Φ_2 = 1/√(2)(
                                                   [       -√(2) (v_1/v) H^+; H_0 - i (v_1/v) A + v_2;                         ]),

 where v_1 = ⟨Φ_1⟩, v_2 = ⟨Φ_2⟩, and v^2 = v^2_1 + v^2_2. 
Here the charged and CP odd interaction states mix with corresponding charged and neutral Goldstone modes through an orthogonal mixing with angle β. 
This mixing produces mass eigenstates H^±, A along with the massless Goldstone bosons that are removed in the unitary gauge. 
The mixing angle is expressed as tanβ = v_2/v_1. 

The particle masses are

 
    m^2_H^± = 1/2(λ_4 +λ_5)v^2,   m^2_A = λ_5 v^2,

 
    m^2_h,H = (λ_1 v^2_1 +λ_2 v^2_2)±√((λ_1 v^2_1 -λ_2 v^2_2)^2
     +(λ_3-λ_4-λ_5)^2 v^2_1 v^2_2).

An immediate consequence of the scenario under consideration is a
very light scalar H with mass,


    m^2_H = 2λ_2 v^2_2[1+O(v_2/v_1)].

The mass eigenstates h, H are related to the weak eigenstates h_0, H_0 by


    h_0 = c_α h + s_α H,  H_0 = -s_α h + c_α H,

where c_α = cosα, s_α = sinα, and
are given by


    c_α   = 1+O(v^2_2/v^2_1),
    
    
    s_α   = -λ_3-λ_4-λ_5/2λ_1(v_2/v_1)+O(v^2_2/v^2_1).

This mixing is extremely tiny as along with the absence of an explicit Z_2 breaking m_12^2 term, v_2∼ 0.1 eV and v_1∼ 250 GeV, 
and thus can be neglected.
The smallness of tanα can be seen from Eq. <ref>, imposing the limit v_2/v_1→ 0, which evidently results in tiny tanβ. 
Hence, the neutral scalar h essentially behaves like the SM Higgs, except of course in
its Yukawa interaction with the neutrinos. 
So we do not expect any observable deviation from the Higgs couplings to the SM particles,
except possibly in the loop induced couplings, e.g. hγγ (through H^± loop). 
The effect of H^± loop in the hγγ vertex is computed in Ref. <cit.>, 
along with the estimation that a sizable Higgs invisible decay is allowed here.
Therefore, the charged fermions obtain mass entirely from the Φ_1 vev v_1, and neutrinos get it through its coupling exclusively to Φ_2.
In the ν2HDM with an imposed lepton number conservation, the neutrino mass is of Dirac nature. 
If that symmetry is not enforced, ν2HDM allows Majorana mass generation of neutrinos with a low scale seesaw mechanism.  

In the ν2HDM with Z_2 symmetry, the second Higgs doublet Φ_2 shows some interesting features. 
The neutral scalars H, A, in their Yukawa interaction, couple almost fully to the neutrinos, while the interactions between neutrinos and 
charged leptons are mediated through H^±. 
Presence of tiny and nonzero neutrino mass can be explained even with the Yukawa couplings at O(1) values. 
Among the H, A coupling to the W and Z bosons, tiny v_2 supresses the triple gauge couplings involving only one scalar. 
The triple gauge couplings with two scalars and one gauge boson can be there and can give rise to a significant scalar pair 
production cross section at the LHC, in the processes like pp→γ^* → H^+H^-.







§ SPLITTING OF ACTIVE NEUTRINOS



In ν2HDM scenario, Φ_1 behaves like the SM scalar doublet while h resembles the SM Higgs. As explained in section <ref>, Φ_1 couples to all the SM 
fermions except the neutrinos with the second doublet Φ_2 and the scalars inside, H,A,H^± interacting in a neutrinophilic way in the Yukawa sector. 
We add a right handed neutrino ν_R, which is odd under Z_2 symmetry. This RH neutrino along with Φ_2, which is also odd under Z_2 forms the Yukawa 
interaction as


    L_Y   = y_e L_e^†Φ̃_̃2̃ν_R + y_μ L_μ^†Φ̃_̃2̃ν_R + y_τ L_τ^†Φ̃_̃2̃ν_R + h.c., 
    
       = ∑_e, μ, τy_e v_2/√(2) U_eiν_i ν_R + ∑_e, μ, τy_e/√(2) U_ei H (i A) ν_i ν_R + other terms , 
       =  m_ν_iν_i ν_R + m_ν_i/v_2 H (i A) ν_i ν_R + other terms ,

where m_ν_i=∑_e, μ, τy_e v_2/√(2) U_ℓ i are neutrino masses and U_ℓ i is the PMNS matrix.




We take the right handed neutrino ν_R mass to be of the order of the SM neutrino mass, i.e. m_N ≈ m_ν and the neutrino Yukawa coupling y_i = m_ν_i/v_2. 
The box diagram of neutrino decay ν→ 3 ν is formed through the neutral scalar (H/A) and right handed neutrino propagators as shown in Fig. <ref>. 
The ν_1 → 3 ν_3 decay through a box diagram is computed in ν2HDM and the decay rate is 

    Γ_ν = 16 π^4 m_ν_1 m_ν_3^3/m_H^4 ( y_1 y_3^3 )^2  |F( m_N^2/m_H^2)|^2  PS_1 → 3,
 
where PS_1 → 3 is the phase space factor of one particle decay to three particles, computed following the Ref. <cit.>. 
The loop factor F( m_N^2/m_H^2) is obtained as 

    F (y) = 1 + y/(1-y)^2 + 2 y log y/(1-y)^3.

Neutrino decay lifetime τ_ν is the inverse of the decay rate computed here.

In the model considered here, three light neutrinos get mass from a neutrinophilic doublet Φ_2 with y_i ∼ O(1). 
We assume a scenario where the heaviest two of these three neutrinos can split to a number of the lightest ones, 
given that the distances they travel from their sources to the place of their detection at IceCube are large enough.  
If three neutrino masses are nearly degenerate, then the daughter neutrino carries almost all of the decaying neutrino energy, and 
thus contributing to the flux at that energy. 
So there will be no depletion in neutrino count in that energy range. 

To show the decrease in neutrino flux, the nature of the mass spectrum must be hierarchical, i.e. having two mass eigenstates with mass much higher than the final state neutrino. 

In this case, the daughter neutrinos after the splitting will be much less energetic than the parent one and therefore, those are left out of the neutrino flux at the energy scale of the splitted neutrino <cit.>.




In this work we take ν_3 to be the lightest final state neutrino. The ν_3- ν_e mixing through tiny PMNS element U_e3, suppresses the ν_e number 
that come out from a ν_3 state, making the ν_e count even smaller. For ν_3 to be the lightest neutrino, inverse neutrino hierarchy is preferred 
where ν_1, ν_2 are relatively heavy (∼ 0.1 eV) and almost with similar mass and the ν_3 is extremely light (∼ 0.01 eV). 
 
The distances of the neutrino sources like AGN and GRB's are of order of 10's of MPc from the earth, where we hope to have the detection. 
For the final state neutrinos to be counted out, the decay factor e^-L/γ c τ should be negligibly small which translates to L/γ c τ≫ 1.
For the neutrinos coming with PeV scale energies,

    L/γ c τ = L/E (mc^2/c τ ) ≫ 1,
 
which corresponds to τ/m < 10^3 s/eV where τ is the neutrino rest frame lifetime.
Big Bang Nucleosynthesis (BBN) provides a lower bound on the neutrino lifetime <cit.>. 
With these limits the allowed window of neutrino decay lifetime is constrained to the range <cit.>,

    10^-6 s/eV ≤τ/m≤ 10^3  s/eV .









§ FLAVOR RATIOS IN NEUTRINO SPLITTING


 We consider three different production mechanism of astrophysical neutrinos here and explore the effects of neutrino splitting on the flavor composition of final 
state neutrinos, both for normal and inverted mass hierarchies. 
The allowed ranges of neutrino flavor ratios at the Earth for different standard and non-standard cases, with arbitrary initial flavor ratios, is 
discussed previously in Ref. <cit.>.
Flavor composition of neutrinos after neutrino oscillation and decay is discussed in Ref. <cit.> earlier. 

Assuming the mass degeneracy of ν_1 and ν_2, for the neutrino splitting rates, we can write Γ(ν_1 → 3 ν_3) ≈Γ(ν_2 → 3 ν_3) for 
inverted mass hierarchy, whereas for normal hierarchy it reads Γ(ν_3 → 3 ν_1) ≈Γ(ν_3 → 3 ν_2). 
Neutrino mixing angles are also assumed as θ_23≈ 45^∘, θ_13≈ 0^∘ and PMNS matrix elements are computed taking these in account,

    U_PMNS = 
    
    [      c_12      s_12         ϵ; s_12/√(2) c_12/√(2)    1/√(2); s_12/√(2) c_12/√(2)   1/√(2), ]

where c_12= cosθ_12, s_12= sinθ_12, with θ_12≈ 32^∘ being the solar neutrino mixing angle, ϵ≈ 0.15. 
The neutrino flux in the mass basis before the decay is 

    ϕ_ν_i = ∑_α |U_i α|^2 ϕ_ν_α.

For normal hierarchy, splitting of ν_3 happens with equal probability to three ν_1, ν_2 states each.
Assuming that the neutrino decay is complete before it reaches the detector, only stable states ν_1, ν_2 are included in the final state flux. 
It also ensures that the exponential energy dependent decay factor e^-L m/E τ does not modify the flux. After the splitting, the flux ratio takes the form 
ϕ_1^d : ϕ_2^d : ϕ_3^d = ϕ_1 + 3/2ϕ_3 : ϕ_2 + 3/2ϕ_3 : 0, as the neutrino flux ϕ_3 splits equally to produce 
three daughter ν_1 or ν_2 from one parent neutrino ν_3. 
Final flux ratio, in the flavor eigenstate takes the form 

    Φ_ν_α = ∑_i |U_α i|^2 Φ_i^d.


The three neutrino sources that are considered here consist of neutrino generation from the pion decay, muon decay and neutron decay.

Pions are the primary products expected from hadronic collisions of accelerated protons with ambient matter in AGN's. The  pions decay  to μ and ν_μ 
and μ further decays to one ν_e and ν_μ each. 
Here the initial flavor ratio of the flux is ϕ_e : ϕ_μ : ϕ_τ = 1:2:0. 

The muon damped sources  the muons from the  π decay is absorbed before its decay to provide only one initial ν_μ with initial flux flavor 
ratio ϕ_e : ϕ_μ : ϕ_τ = 0:1:0.

Another source of neutrinos is the neutron decay where one ν_e is originated in n → p e ν_e decay, with initial flavor flux 
ratio ϕ_e : ϕ_μ : ϕ_τ = 1:0:0.

For all these three cases initial neutrino flux at mass basis is computed using Eq. <ref> and Eq. <ref>. 
Neutrino flux ratio after the neutrino splitting for the case of normal hierarchy is computed using Eq. <ref> and Eq. <ref>. 
For normal hierarchy, neutrinos coming from pion decay split like ν→ 3 ν to produce flavor ratio of daughter neutrinos as 
ϕ_e : ϕ_μ : ϕ_τ≈ 2:1:1. When neutrinos come from muon damped sources with a normal hierarchy, they split to provide final flux ratio in the 
flavor basis as ϕ_e : ϕ_μ : ϕ_τ≈ 1.8:1:1. For the neutrinos originating from the neutron decay, their splitting results in the final flux 
ratio ϕ_e : ϕ_μ : ϕ_τ≈ 3:1:1.

In all these cases, for inverted mass hierarchy m_ν_1≈ m_ν_2 >>  m_ν_3 the final neutrino 
mass eigenstates are purely ν_3, which in flavor basis contribute to the neutrino flux as ϕ_μ : ϕ_τ≈ 1:1, with ν_e flux being negligible. 
Therefore with neutrinos with inverted hierarchy, the flavor ratio of the daughter neutrinos after splitting is ϕ_e : ϕ_μ : ϕ_τ≈ 0:1:1 for 
any initial flavor ratio. 





§ OBLIQUE PARAMETERS: S, T, U


The radiative corrections to the gauge boson two point functions (Π_VV) are known as oblique corrections. 
Peskin-Takeuchi oblique parameters manely S,T, U are defined by different combination of those modified two point functions.  
The S parameter measures the running of two point functions of the neutral gauge bosons (ZZ, Zγ and γγ) between zero
momentum and the Z pole with p^2= m_z^2. 
Therefore, it will be specially sensitive to new degrees of freedom at low scales, particularly below the m_Z. 
Thus, it becomes very important in the presence of an ultra light neutral scalars, which is the case for the ν2HDM model with exact Z_2 symmetry.  
The T parameter measures the amount of custodial symmetry breaking at zero momentum, i.e. the difference between the WW and the ZZ two point functions at q^2=0. 
It usually plays a significant role in constraining the parameter space of particles charged under SU(2)_L. 
Splitting the masses of particles inside a doublet breaks custodial symmetry and affects T. 
The impact of adding a second Higgs doublet in the electroweak precision tests (EWPT), that is encoded in the
oblique parameters S, T, and U, has been discussed in detail in the literature <cit.>.













The experimentally measured values of these oblique parameters from the LHC are:

    Δ S^SM = 0.05±0.11,
       Δ T^SM = 0.09±0.13,
       Δ U^SM = 0.01±0.11,











In the Z_2 symmetric ν2HDM, scalar spectrum of the model is very tightly constrained. 
Unbroken Z_2 symmetry leads to m_12^2=0 and the mixing angles are also very small with sin^2α,sin^2β≪ 1. 
While one CP-even scalar h is taken as the 125 GeV Higgs particle discovered at the LHC, and thus has fixed λ_1. 
The neutrinophilic CP-even neutral scalar H is extremely light, m_H ∼eV (O(v_2)) ≪ v. For modest values of the quartic
couplings within the perturbative limit, the charged scalars and the CP-odd scalar masses are expected to be below the TeV scale.
As the presence of a light neutral scalar is essential to explain sufficient depletion of neutrino flux and that kind of scalar is present here, 
it turns out that the oblique parameters S and T will play a decisive role in constraining the model.

The modification of the S parameter due to the presence of an ultra light neutral scalar, m_H ∼ O(v_2), pushes it to a large negative value to 
rule out the model at 2 σ confidence level.

We modify the Z_2 symmetric ν2HDM with the addition of one vectorlike lepton SU(2)_L doublet (ψ) with hypercharge Y_ψ and one 
vectorlike lepton SU(2)_L singlet χ with hypercharge Y_χ=Y_ψ - 1/2. Both the lepton doublet and singlet are Z_2 even. 
So the additional terms that will be added to the ν2HDM Lagrangian being consistent with all the symmetries of the model is

    L = m_ψψ̅ψ + m_χχ̅χ - y_χψ̅Φ_1 χ + h.c. ,
 
After electroweak symmetry breaking ψ_2 and χ mix with each other to provide two lepton mass eigenstates. 
These two along with ψ_1 contribute to the S,T parameter calculation. 
The values of these parameters for the cases with and without vectorlike leptons are shown in Fig. <ref>. 
Phenomenology of this kind of construction is discussed in detail at Ref. <cit.>.






§ OTHER CONSTRAINTS ON Ν2HDM


The charged Higgs (H^+) from 2HDM has interesting phenomenology at the LHC. 
The p p → H^+H^- pair production happens at the LHC through the exchange of off-shell γ and Z. 
This 2HDM production rate remains same for the ν2HDM charged Higgs. 
As explained before, the second scalar doublet here has a very small mixing with the scalars of the SM-like doublet. 
The quark decay modes of the charged scalar are highly suppressed with the mixing factor v_2/v. 
So the charged scalar decays mostly to the leptonic channels H^±→ lν. 
Corresponding bound on the charged Higgs mass from the LEP in this channel is given as m_H^± > 80 GeV.

The charged scalars can contribute to muon and electron g-2 calculation at one loop, but the contribution is negligible due to a 
suppression factor m_l^4/m_H^±^2 in the amplitude. 
Unlike a general 2HDM where the two loop contributions can be dominant, charged lepton couplings to H, A are suppressed by a factor tanβ in the ν2HDM, 
making g-2 constraints insignificant in this scenario. 

The charged scalar loop contribution to h →γγ is of the same sign as the W loop while top loop has a greater contribution with opposite sign. 
Therefore, the new H^± decreases the diphoton decay amplitude, though the reduced amplitude still satisfies the experimental limit. 

There are strong experimental constraints on the lepton flavor violating processes through the loop diagrams, mainly in the μ→ e
γ channel. For this process, the charged scalar mediated decay branching ratio
reads <cit.>


    BR(μ→ e γ) = BR( μ→ eν̅ν) α_EM/192 π |⟨ m_μ e^2 ⟩|^2 ρ^2


where ρ = (G_F m_H^±^2 v_2^2)^-1.
The strongest bound on this lepton flavor violating (LFV) decay comes from
the MEG-2 experiment which gives the upper limit as BR(μ→ e γ)<5.7 ×
10^-13 <cit.>, while bounds are weaker for the other LFV decay channels (τ→ e γ, τ→μγ),
that are vobtained from the BaBar Collaboration.  
In terms of ρ defined earlier, the 90% confidence level bounds read <cit.>,


    ρ≲ 1.2  eV^-2   for    μ→ e γ .

The limit on ρ translates into a limit where for v_2 ≲ 1 eV one must have the lower bound as m_H^±≳ 250 GeV.
This is the tightest limit on the v_2 and m_H^± parameter space upto now. 


With the sensitivity of the MEG expected to be improved further, the bound on ρ is expected to be improved by about one order of magnitude. 
The limits imposed due to the MEG bound on the (m_H^±, v_2) plane are shown in Fig. 1 of Ref. <cit.>. 
During the big bang nucleosynthesis (BBN) era 
the ultra light neutral scalar H cannot stay at equilibrium with thermal bath as H with m_ν > 0.1 eV decays strongly through the H → 2 ν mode in the ν2HDM. 
Therefore, the ultra light bosonic mode which will otherwise be relativistic, does not contribute to the effctive degrees of freedom N_eff in the BBN era.
Astrophysical constraints on the neutrinophilic 2HDM is discussed in the Ref <cit.>.





§ ICECUBE ANALYSIS


IceCube is one prominent neutrino telescope which can probe high energetic astrophysical 
neutrino emission from supernovae, gamma ray bursts (GRB), active 
galactic nuclei (AGN) or from other possible new cosmic sources. 
IceCube experiment, after a run of six years, has recorded a total of 82 events combining 
both the track and shower events with neutrino interaction vertex inside the detector. Track events in the IceCube detector arise from the charged current interactions of 
muon neutrinos whereas shower events happen due to charged currents of ν_e and ν_τ and neutral currents of all neutrino flavor. These neutrino events observed 
at IceCube are together called 'High Energy Starting Events' (HESE). 
We assume the presence of a single astrophysical neutrino flux along with atmospheric 
components to analyze the observed diffuse neutrino flux. 
Flavor components of isotropic astrophysical neutrino flux can be different for different astrophysical sources, as discussed in section <ref>. 
The flavor ratios of the flux after neutrino spliting are also obtained there. 
If we take only the shower events, the daughter neutrinos are dominantly of electron flavor. 
The all flavor initial astrophysical neutrino spectrum with single component is parametrized as 

    Φ (E_ν) = Φ_astro^0(E_ν/100  TeV)^-γ

When the astrophysical neutrinos after traveling a distance L from the source with energy E_ν gets splitted with decay lifetime τ, the disappearance 
of those results in the suppression of the neutrino flux. So, just after the decay the neutrino energy spectrum takes the form

    Φ_decay (E_ν) = Φ_astro^0(E_ν/100  TeV)^-γ e^-β/E_ν,

when β for a neutrino with mass m is given as 

    β = L m/τ.






 §.§ Splitting of Neutrinos



















In the neutrinophilic two Higgs doublet model (ν2HDM) the active neutrinos do not decay to the invisible final state particles. 
Instead, one heavier neutrino mass eigenstate splits to three lighter neutrino mass eigenstates with a process like ν_3 → 3 ν_1, 
creating three daughter neutrinos for each initial neutrinos. 


To account for the general time dependence in the neutrino propagation, we take the case where some of the neutrinos have decayed while rest of them have not till that time.
Part of the neutrino flux that is decayed will be increased three times after the decay here, compared to the case when neutrinos just vanish due to 
their decay to invisible particles. 
The amount of the initial flux that remains unchanged after the decay is of the form 

    Φ_nd (E_ν) = Φ_astro^0(E_ν/100  TeV)^-γ e^-β/E_ν.

The daughter neutrinos originating after the decay will have different energies compared to the parent one and we denote that 
energy as E_d. We assume all three daughter neutrinos will have same energy i.e. E_d ≈ E_ν/3 and they will interact 
with the IceCube detector with that energy.
So, the part of the flux that is decayed to three neutrinos each, can be written as 

    Φ_d (E_ν) = Φ_astro^0(E_ν/100  TeV)^-γ( 1- e^-β/E_ν).

For a particular bin of detection energy in the IceCube, the daughter neutrino energy will have a lower bound equal to the detection energy. 
With all these considerations in mind we write the total number of astrophysical neutrinos being observed at a particular energy bin as <cit.>

    N_i = 4 πΔ t n_int∫^E_i+1_E_i dE [ ∫^∞_EΦ_nd (E_ν) d σ/dE (E_ν)  dE_ν
                + ∫^∞_E 3  Φ_d (E_ν) d σ/dE (E_d) dE_d ],

when d σ/dE = ( dσ^NC/dE + dσ^CC/dE) is the total neutral current (NC) and charged current (CC) interaction.

 This represents the effective interaction area per effective interacting nucleon present inside the matter of the IceCube experiment. 
Here n_int is the total effective number of nucleon-neutrino interaction points at the IceCube, where Δ t is the total time 
exposure of the IceCube detector i.e. the duration of IceCube run of 2078 days.

 






 









The interaction cross-section can be written as <cit.>, dσ/dE = 1/E_ddσ (E_d)/dy, where y 
is the energy fraction of final state lepton and can be written here as E_l/E_d and E is the energy deposit. 
A complete shower acceptance simulation, taking into account the difference of energy transfer in CC scattering of ν_e and ν_τ along with 
unequal hadronic content in the final states, is beyond the scope of this work. Following the Ref. <cit.>, full transfer of ν_e energy into showers is 
approximated as, E_d ≈ E that can be used to write dσ/dE = σ^CC (E^e_d) δ(E-E^e_d). Taking into account 50 % energy 
transfer for ν_τ initiated CC interaction we write dσ/dE = σ^CC (E^τ_d) δ(E- 0.5 E^τ_d). For all neutrino flavors, 
neutral current interaction is written as dσ/dE = σ^NC (E_d) δ(E- 0.2 E_d). The energy dependent ν_e-nucleon charged current (CC) and 
neutral current (NC) interaction cross section σ(E_d) in the SM case can be extracted from Fig. 3 of Ref. <cit.>.   














We choose some benchmark points of the relevant parameters like the astrophysical neutrino flux amplitude Φ_astro, the spectral index 
γ and β that quantifies the amount of neutrino decay. It is shown in Fig. <ref> how an active neutrino splitting to three active neutrinos 
can have a significant impact on the shape of the distribution of IceCube events. IceCube event distribution in the energies from 250 
TeV to 1 PeV are better explained if the the neutrino spectrum is relatively steeper. A moderate astrophysical neutrino flux amplitude along with a high spectral
index results in heavy suppression of the number of expected events at energies above 1 PeV causing a large deviation from the observed IceCube neutrino events. 
If the neutrino flux amplitude is increased for a steep spectrum it can explain observed events at energy bins above 250 TeV but the deviation is there 
in the lower energy (60-250 TeV) bins, as the expected number of neutrino events increase with flux amplitude. 
This issue can be resolved introducing a split of active neutrinos and therefore suppressing the energy spectrum with the factor e^-β/E_ν with 
greater suppression at lower energy bins, which therefore can help expected neutrino events to match the observed ones. The plots in Fig. <ref> present effects of 
the introduction of a splitting of neutrinos for different benchmark points comparing with the case where there is no neutrino splitting. Amount of suppression that 
the decay brings to the neutrino event distribution is dependent on β, and increases with it. The left plot with β = 3.8 × 10^5GeV suppresses 
the event distribution of the neutrino without decay case more than that in the right plot with β = 2 × 10^4 GeV. 


We quantify the effects of neutrino decay on the χ^2 values in our statistical analysis here. As the χ^2 values by themselves 
are not too informative, we try to compare the χ^2 values with neutrino decay with its values when there way no neutrino decay, over the parameter 
region that is interesting for IceCube phenomenology. Therefore, we use a metric 
that can compare the two cases of absence and presence of neutrino splitting and probe if there 
is a better fit for the case of neutrino decay. We define the quantity as,

    Δχ^2 = 100 ×χ^2( no  decay) - χ^2( with  decay)/χ^2( no  decay),

and compute its values with variation of different model parameters.  The variation of Δχ^2 over different parameter region is presented in Fig. <ref>.







 §.§ Invisible Decay of neutrinos


Decay of active neutrinos with absence of any SM neutrinos in the final state decay products is termed as neutrino invisible decay. In the ν2HDM model discussed here,
neutrino decay is not invisible. Invisible decays are possible in different scenarios like lepton number violating triplet Majoron models <cit.> and 
singlet Majoron models <cit.>. Therefore, we like to explore the observable effects of invisible neutrino decay at the IceCube. 
We start with an isotropic astrophysical neutrino flux shown above in Eq. <ref>.   
The expression of neutrino flux in Eq. <ref> depicts the suppression of neutrino flux due its invisible decay. 
Using that expression, total number of astrophysical neutrinos being observed at a particular energy bin [E_i, E_i+1] as 

    N_i = 4 πΔ t n_int∫^E_i+1_E_i dE ∫^∞_EΦ_decay (E_ν) ( dσ^NC/dE + dσ^CC/dE) dE_ν ,

when ( dσ^NC/dE + dσ^CC/dE) is the neutrino-nucleon interaction as discussed and computed above. 


How the distribution of events are modified at IceCube in presence of an invisible neutrino decay is shown in Fig. <ref>. For a flat energy 
spectrum of astrophysical neutrinos, IceCube events in multi-PeV bins are satisfied along with those in the energies upto 250 TeV. That results in the 
intermediate energy bin event distribution being far off from the expected results. Similarly, a flat energy spectrum which satisfies events in all bins above 250 TeV,
falls short of observed events in the lower energy bin. A relatively steeper neutrino energy spectrum with higher γ can satisfy all the energy bin above 250 TeV, 
but with high spectrum amplitude that results in excess events at lower energy (60 -250 TeV) bins. Introduction of an invisible neutrino decay can suppress the flux 
at lower energy to satisfy the IceCube event distribution. 



§ SUMMARY AND CONCLUSION



The anomalies in the IceCube data have been  addressed in terms of excess at PeV energies or a cutoff at 6 PeV Glashow resonance energies. In this paper we discuss neutrino splitting which transfers neutrino energies to roughly 1/3 after decay and due to the Lorentz factor there is a larger depletion of lower energy neutrinos. 
The two effects together give an explanation of the PeV excess and 6 PeV depletion in a single power law initial neutrino flux scenario. The benchmark points that explain the IceCube spectrum put the decay parameter in the range, β=L/(τ/m)= (5 × 10^4 - 1.25 × 10^5) GeV with the correlated spectral index being in the range
γ = (2.8-3) as shown in Table <ref>. The χ^2 values for benchmark points in the case with no neutrino splitting or decay along with the case when they are present are also shown in Table <ref>. 
To achieve a decay at L=1 Mpc of m_ν=0.1 eV neutrinos we need τ= (3 × 10^-2 - 2 × 10^-1 ) second. 
To achieve this value of τ/m we need a small mass mediator. This small mass scalar occurs naturally in a ν2HDM model where the vev of the scalar is taken to O(eV) in order to explain the neutrino mass. Presence of an ultra light scalar along with the mass splitting between the neutral and charged component of the second Higgs doublet gives rise to a large contribution to the S and T oblique parameters. This problem is fixed here by adding a vector lepton doublet and singlet with mass 200 GeV. The neural component of the heavy fermions can serve as dark matter as has been studied in <cit.>. The charged and CP-odd neutral scalars and charged vector-like leptons can potentially be observed at the LHC, along with interesting phenomenological properties <cit.>.

For the neutrinos splitting ν→ 3 ν scenario the daughter neutrinos are in the ratio ϕ_e : ϕ_μ : ϕ_τ≈ 0:1 :1 for the case of inverted hierarchy and for normal hierarchy the ratio is ϕ_e : ϕ_μ : ϕ_τ≈ 2:1 : 1 with pion decay as neutrino source, ϕ_e : ϕ_μ : ϕ_τ≈ 1.8:1 : 1 for muon damped source of neutrinos and ϕ_e : ϕ_μ : ϕ_τ≈ 3: 1 : 1, when neutrinos come from the neutron decay. So the flavor ratio determination <cit.> can be a good test of the neutrino splitting scenario which can also give information of the mass hierarchy of neutrinos and distinguish the neutrino splitting scenario from other phenomenon like invisible decay of neutrinos <cit.>.


SS acknowledges the useful discussion with Bhavesh Chauhan on the effects of an ultra light scalar on the BBN constraints. 



99


Aartsen:2017mau 
  M. G. Aartsen et al. [IceCube Collaboration],
  
  arXiv:1710.01191 [astro-ph.HE].
  
  


Aartsen:2015rwa 
  M. G. Aartsen et al. [IceCube Collaboration],
  
  Phys. Rev. Lett.  115, no. 8, 081102 (2015)
  doi:10.1103/PhysRevLett.115.081102
  [arXiv:1507.04005 [astro-ph.HE]].
  
  


Cholis:2012kq 
  I. Cholis and D. Hooper,
  
  JCAP 1306, 030 (2013)
  doi:10.1088/1475-7516/2013/06/030
  [arXiv:1211.1974 [astro-ph.HE]].
  
  


Anchordoqui:2013dnh 
  L. A. Anchordoqui et al.,
  
  JHEAp 1-2, 1 (2014)
  doi:10.1016/j.jheap.2014.01.001
  [arXiv:1312.6587 [astro-ph.HE]].
  
  


Murase:2014tsa 
  K. Murase,
  
  AIP Conf. Proc.  1666, 040006 (2015)
  doi:10.1063/1.4915555
  [arXiv:1410.3680 [hep-ph]].
  
  


Sahu:2014fua 
  S. Sahu and L. S. Miranda,
  
  Eur. Phys. J. C 75, 273 (2015)
  doi:10.1140/epjc/s10052-015-3519-1
  [arXiv:1408.3664 [astro-ph.HE]].
  
  


Waxman:1998yy 
  E. Waxman and J. N. Bahcall,
  
  Phys. Rev. D 59, 023002 (1999)
  doi:10.1103/PhysRevD.59.023002
  [hep-ph/9807282].
  
  


Esmaili:2013gha 
  A. Esmaili and P. D. Serpico,
  
  JCAP 1311, 054 (2013)
  doi:10.1088/1475-7516/2013/11/054
  [arXiv:1308.1105 [hep-ph]].
  
  


Murase:2015gea 
  K. Murase, R. Laha, S. Ando and M. Ahlers,
  
  Phys. Rev. Lett.  115, no. 7, 071301 (2015)
  doi:10.1103/PhysRevLett.115.071301
  [arXiv:1503.04663 [hep-ph]].
  
  


Dev:2016qbd 
  P. S. B. Dev, D. Kazanas, R. N. Mohapatra, V. L. Teplitz and Y. Zhang,
  
  JCAP 1608, no. 08, 034 (2016)
  doi:10.1088/1475-7516/2016/08/034
  [arXiv:1606.04517 [hep-ph]].
  
  


Bhattacharya:2014yha 
  A. Bhattacharya, R. Gandhi and A. Gupta,
  
  JCAP 1503, no. 03, 027 (2015)
  doi:10.1088/1475-7516/2015/03/027
  [arXiv:1407.3280 [hep-ph]].
  
  



Barger:2013pla 
  V. Barger and W. Y. Keung,
  
  Phys. Lett. B 727, 190 (2013)
  doi:10.1016/j.physletb.2013.10.021
  [arXiv:1305.6907 [hep-ph]].
  
  



Dutta:2015dka 
  B. Dutta, Y. Gao, T. Li, C. Rott and L. E. Strigari,
  
  Phys. Rev. D 91, 125015 (2015)
  doi:10.1103/PhysRevD.91.125015
  [arXiv:1505.00028 [hep-ph]].
  
  



Dey:2015eaa 
  U. K. Dey and S. Mohanty,
  
  JHEP 1604, 187 (2016)
  doi:10.1007/JHEP04(2016)187
  [arXiv:1505.01037 [hep-ph]].
  
  


Chauhan:2017ndd 
  B. Chauhan, B. Kindra and A. Narang,
  
  arXiv:1706.04598 [hep-ph].
  
  



Glashow:1960zz 
  S. L. Glashow,
  
  Phys. Rev.  118, 316 (1960).
  doi:10.1103/PhysRev.118.316
  
  


Berezinsky:1977sf 
  V. S. Berezinsky and A. Z. Gazizov,
  
  JETP Lett.  25, 254 (1977)
  [Pisma Zh. Eksp. Teor. Fiz.  25, 276 (1977)].
  
  


Berezinsky:1981bt 
  V. S. Berezinsky and A. Z. Gazizov,
  
  Sov. J. Nucl. Phys.  33, 120 (1981)
  [Yad. Fiz.  33, 230 (1981)].
  
  


Bhattacharya:2011qu 
  A. Bhattacharya, R. Gandhi, W. Rodejohann and A. Watanabe,
  
  JCAP 1110, 017 (2011)
  doi:10.1088/1475-7516/2011/10/017
  [arXiv:1108.3163 [astro-ph.HE]].
  
  


Joshipura:2013yba 
  A. S. Joshipura, S. Mohanty and S. Pakvasa,
  
  Phys. Rev. D 89, no. 3, 033003 (2014)
  doi:10.1103/PhysRevD.89.033003
  [arXiv:1307.5712 [hep-ph]].
  
  



Pakvasa:2012db 
  S. Pakvasa, A. Joshipura and S. Mohanty,
  
  Phys. Rev. Lett.  110, 171802 (2013)
  doi:10.1103/PhysRevLett.110.171802
  [arXiv:1209.5630 [hep-ph]].
  
  


Tomar:2015fha 
  G. Tomar, S. Mohanty and S. Pakvasa,
  
  JHEP 1511, 022 (2015)
  doi:10.1007/JHEP11(2015)022
  [arXiv:1507.03193 [hep-ph]].
  
  



Anchordoqui:2014hua 
  L. A. Anchordoqui, V. Barger, H. Goldberg, J. G. Learned, D. Marfatia, S. Pakvasa, T. C. Paul and T. J. Weiler,
  
  Phys. Lett. B 739, 99 (2014)
  doi:10.1016/j.physletb.2014.10.037
  [arXiv:1404.0622 [hep-ph]].
  
  


Liao:2017yuy 
  J. Liao and D. Marfatia,
  
  arXiv:1711.09266 [hep-ph].
  


Cohen:2011hx 
  A. G. Cohen and S. L. Glashow,
  
  Phys. Rev. Lett.  107, 181803 (2011)
  doi:10.1103/PhysRevLett.107.181803
  [arXiv:1109.6562 [hep-ph]].
  
  



Machado:2015sha 
  P. A. N. Machado, Y. F. Perez, O. Sumensari, Z. Tabrizi and R. Z. Funchal,
  
  JHEP 1512, 160 (2015)
  doi:10.1007/JHEP12(2015)160
  [arXiv:1507.07550 [hep-ph]].
  
  


Gabriel:2006ns 
  S. Gabriel and S. Nandi,
  
  Phys. Lett. B 655, 141 (2007)
  doi:10.1016/j.physletb.2007.04.062
  [hep-ph/0610253].
  
  


Davidson:2009ha 
  S. M. Davidson and H. E. Logan,
  
  Phys. Rev. D 80, 095008 (2009)
  doi:10.1103/PhysRevD.80.095008
  [arXiv:0906.3335 [hep-ph]].
  
  


Seto:2015rma 
  O. Seto,
  
  Phys. Rev. D 92, no. 7, 073005 (2015)
  doi:10.1103/PhysRevD.92.073005
  [arXiv:1507.06779 [hep-ph]].
  
  



Romao:1998sr 
  J. C. Romao and S. Andringa,
  
  Eur. Phys. J. C 7, 631 (1999)
  doi:10.1007/s100529801038
  [hep-ph/9807536].
  
  


Bustamante:2015waa 
  M. Bustamante, J. F. Beacom and W. Winter,
  
  Phys. Rev. Lett.  115, no. 16, 161302 (2015)
  doi:10.1103/PhysRevLett.115.161302
  [arXiv:1506.02645 [astro-ph.HE]].
  
  


Pakvasa:2004hu 
  S. Pakvasa,
  
  Mod. Phys. Lett. A 19, 1163 (2004)
  [Yad. Fiz.  67, 1179 (2004)]
  doi:10.1134/1.1772452
  [hep-ph/0405179].
  
  



Grimus:2007if 
  W. Grimus, L. Lavoura, O. M. Ogreid and P. Osland,
  
  J. Phys. G 35, 075001 (2008)
  doi:10.1088/0954-3899/35/7/075001
  [arXiv:0711.4022 [hep-ph]].
  
  


Haber:2010bw 
  H. E. Haber and D. O'Neil,
  
  Phys. Rev. D 83, 055017 (2011)
  doi:10.1103/PhysRevD.83.055017
  [arXiv:1011.6188 [hep-ph]].
  
  


Ellis:2014dza 
  S. A. R. Ellis, R. M. Godbole, S. Gopalakrishna and J. D. Wells,
  
  JHEP 1409, 130 (2014)
  doi:10.1007/JHEP09(2014)130
  [arXiv:1404.4398 [hep-ph]].
  
  


Bertuzzo:2015ada 
  E. Bertuzzo, Y. F. Perez G., O. Sumensari and R. Zukanovich Funchal,
  
  JHEP 1601, 018 (2016)
  doi:10.1007/JHEP01(2016)018
  [arXiv:1510.04284 [hep-ph]].
  
  


Adam:2013mnn 
  J. Adam et al. [MEG Collaboration],
  
  Phys. Rev. Lett.  110, 201801 (2013)
  doi:10.1103/PhysRevLett.110.201801
  [arXiv:1303.0754 [hep-ex]].
  
  



Sher:2011mx 
  M. Sher and C. Triola,
  
  Phys. Rev. D 83, 117702 (2011)
  doi:10.1103/PhysRevD.83.117702
  [arXiv:1105.4844 [hep-ph]].
  
  


Dey:2016sht 
  U. K. Dey, S. Mohanty and G. Tomar,
  
  arXiv:1606.07903 [hep-ph].
  
  



Kistler:2013my 
  M. D. Kistler, T. Stanev and H. Yüksel,
  
  Phys. Rev. D 90, no. 12, 123006 (2014)
  doi:10.1103/PhysRevD.90.123006
  [arXiv:1301.1703 [astro-ph.HE]].
  
  


Ahlers:2006pf 
  M. Ahlers, J. Kersten and A. Ringwald,
  
  JCAP 0607, 005 (2006)
  doi:10.1088/1475-7516/2006/07/005
  [hep-ph/0604188].
  
  


Gelmini:1980re 
  G. B. Gelmini and M. Roncadelli,
  
  Phys. Lett.  99B, 411 (1981).
  doi:10.1016/0370-2693(81)90559-1
  
  


Chikashige:1980qk 
  Y. Chikashige, R. N. Mohapatra and R. D. Peccei,
  
  Phys. Rev. Lett.  45, 1926 (1980).
  doi:10.1103/PhysRevLett.45.1926
  
  


Bhattacharya:2015qpa 
  S. Bhattacharya, N. Sahoo and N. Sahu,
  
  Phys. Rev. D 93, no. 11, 115040 (2016)
  doi:10.1103/PhysRevD.93.115040
  [arXiv:1510.02760 [hep-ph]].
  
  


Kumar:2015tna 
  N. Kumar and S. P. Martin,
  
  Phys. Rev. D 92, no. 11, 115018 (2015)
  doi:10.1103/PhysRevD.92.115018
  [arXiv:1510.03456 [hep-ph]].
  
  


Wang:2016vfj 
  W. Wang and Z. L. Han,
  
  Phys. Rev. D 94, no. 5, 053015 (2016)
  doi:10.1103/PhysRevD.94.053015
  [arXiv:1605.00239 [hep-ph]].
  
  
 

Gopalakrishna:2015wwa 
  S. Gopalakrishna, T. S. Mukherjee and S. Sadhukhan,
  
  Phys. Rev. D 93, no. 5, 055004 (2016)
  doi:10.1103/PhysRevD.93.055004
  [arXiv:1504.01074 [hep-ph]].
  
  
 
  
Aartsen:2015ivb 
  M. G. Aartsen et al. [IceCube Collaboration],
  
  Phys. Rev. Lett.  114, no. 17, 171102 (2015)
  doi:10.1103/PhysRevLett.114.171102
  [arXiv:1502.03376 [astro-ph.HE]].
  
  



Maltoni:2008jr 
  M. Maltoni and W. Winter,
  
  JHEP 0807, 064 (2008)
  doi:10.1088/1126-6708/2008/07/064
  [arXiv:0803.2050 [hep-ph]].
  
  
  
  
Baerwald:2012kc 
  P. Baerwald, M. Bustamante and W. Winter,
  
  JCAP 1210, 020 (2012)
  doi:10.1088/1475-7516/2012/10/020
  [arXiv:1208.4600 [astro-ph.CO]].
  
  
  
