Activating the 4th Neutrino of the 3+1 Scheme

Non-Standard Interactions (NSI) of neutrinos with matter has received renewed interest in recent years. In particular, it has been shown that NSI can reconcile the $3+1$ solution with IceCube atmospheric data with $E_\nu>500$ GeV, provided that the effective coupling of NSI is large, e.g. $\sim 6 G_F$. The main goal of the present paper is to show that contrary to intuition, it is possible to build viable models with large NSI by invoking a new $U(1)$ gauge symmetry with gauge boson of mass $\sim 10$ eV. We refer to these new constructions as $3+1+ U(1)$ models. In the framework of a $3+1$ solution to LSND and MiniBooNE anomalies, we show that this novel NSI can help to solve the tension with cosmological bounds and constraints from IceCube atmospheric data with $E_\nu>500$ GeV. We then discuss the implications of the MINOS and MINOS+ results for the 3+1+$U(1)$ scenario.


Introduction
Adding a light sterile neutrino mixed with active neutrinos to the standard model is one of the most economic extensions which leads to a rich phenomenology. The longstanding LSND anomaly [1] has been largely confirmed by MiniBooNE [2] recently reaching over 6σ evidence by combining the two data sets. Along with the reactor [3] and Gallium [4] anomalies there is a simple solution to these anomalies within the so-called 3+1 scheme which requires a sterile neutrino of mass ∼ 1 eV with a mixing with active neutrinos, θ ∼ O(0.1). This scheme is however in serious tension with the observation of atmospheric neutrinos by IceCube and with cosmological constraints on the presence of new light neutrinos in the early universe.
Within the standard 3+1 scheme, the propagation in matter is governed by the following Hamiltonian where V e = √ 2G F N e − ( √ 2/2)G F N n and V µ = V τ = −( √ 2/2)G F N n and U is the 4 × 4 unitary mixing matrix including the sterile neutrino. If sterile neutrinos have no interaction with matter, the corresponding matter potential vanishes, V s = 0. For high energy neutrinos, ∆m 2 41 /(2E) can be of order of V µ , leading to a resonant enhancement of the ν µ → ν s oscillation which can be probed by atmospheric neutrinos at IceCube [5][6][7][8][9][10][11][12], see fig. 2. Null results from IceCube on the deviation of P (ν µ → ν µ ) from the standard 3ν scheme prediction set a bound on active sterile mixing which is in tension with the value derived from LSND [13]. Ref. [14] suggested turning on NSI for ν µ and ν τ to suppress the effective mixing for E ν > 500 GeV. However, as shown in [14], this requires values of NSI which are larger than the standard weak coupling. It is very challenging to build viable models with such large NSI couplings that satisfy various bounds on the couplings of neutrinos and active neutrinos. However, the bounds on the coupling of sterile neutrinos are more relaxed so it is intriguing to entertain the possibility of a fourth neutrino with sizeable interaction with matter fields even if MiniBooNE and LSND anomalies are refuted by the upcoming SNB experiment.
Let us suppose the sterile neutrino (in the flavor basis) has an effective interaction of the following form with matter fields (f = u, d) where P L ≡ (1 − γ 5 )/2 is the left-handed projector. The value of the sterile matter potential in Earth can be then estimated as V s = 2.5 × 10 −12 eV ρ 5 gr cm −3 , which should be compared to 10 −12 eV (∆m 2 /eV 2 )(TeV/E). The active sterile mixing in Earth will be suppressed for neutrinos with energy larger than TeV if 10. On the other hand, to reproduce LSND and MiniBooNE we need unsuppressed mixing for energies smaller than few hundred MeV so we obtain an upper bound on < 10 4 . Intermediate values of can dramatically affect the long baseline and atmospheric neutrino data in the energy range 10 GeV < E < 100 GeV which has been detected by IceCube. Thus, to solve the tension, we will focus on | | 10. According to the MiniBooNE collaboration, the data can be better explained by an agnostic introduction of "new" matter effects with a resonance energy of 300 MeV 1 . We shall comment on this possibility within the framework that we are discussing.
Another tension shows up between the 3 + 1 solution to the short baseline neutrino observation and cosmological bounds. It has been shown in a series of papers [16][17][18][19][20][21][22][23][24][25][26] that self-interaction of sterile neutrinos can ease this tension. We shall discuss under what conditions the interactions that we are discussing can solve the tension. Long baseline NOvA [27] and MINOS and MINOS+ [28] experiments can also constrain the 3+1 scenario. In fact, according to [28], the constraints from MINOS+ can rule out a significant part of the 3+1 solution to LSND. However, the strength of these bounds are debated in [29]. We study the bounds that MINOS+ can set on various combinations of θ 14 , θ 34 , and θ 24 with and without interaction of sterile neutrinos with matter.
The paper is organized as follows. In section 2, we enumerate the various bounds that already exist on the new couplings and show that despite these bounds, it is still possible to obtain 1. In section 3, we describe a model that can lead to the couplings we are interested in. We show that the region of the parameter space that leads to 1 is natural in the sense that it does not suffer from fine tuning. In section 4, we first quantify how turning on can reconcile TeV range atmospheric neutrino data with the 3+1 solution to the LSND. We then study the bounds from MINOS and MINOS+. In section 5, we summarize our results.

Bounds on new gauge couplings of neutrinos and baryons
In this section, we outline the various bounds that constrain the couplings of Z to baryons (g B ) and to a sterile neutrino ν s (g s ). Let us start with the bounds on g B . These bounds are shown in Fig. 1.
Bound from fifth force: For m Z < 10 eV, the most stringent bound comes from the socalled 5th force search experiments. By introducing this new interaction, two test objects located at a distance of r from each other will exert a new force to one another with a potential given by where m i , A i and Z i (with i = 1, 2) are respectively the mass, mass number and atomic number of the test objects. Here, a e represents the ratio of the coupling to the electron to g B . For example, in case of B − L gauge symmetry, a e = −1. Since the test objects are electrically neutral, Z i gives the numbers of protons as well as that of electrons in atoms of test objects. For m Z ∼ 10 eV (i.e., range of 2 × 10 −6 cm), the measurement of the n e u t r o n -P b Figure 1. Here we display some of the most relevant bounds on the B − L gauge coupling in the mediator mass range of interest. These include stellar cooling in the Sun and horizontal branch stars in globular clusters [30,31], bounds from n−Pb scattering data [32], meson decays [33,34], and fifth force searches [35].
Casimir force (see fig. 27 of [35]) sets a strong bound on g B as shown in fig. 1 with a red line, marked with fifth force. Stellar cooling: Another strong bound on the model comes from the cooling rate of the Sun and Horizontal Branch (HB) stars. The Z boson can be produced in the stellar core via plasmon effect. For m Z ∼ 100 eV (for m Z ∼ keV) the production of both transverse and longitudinal Z in the sun (in HB stars) can be on-shell. Thus, for this mass range stellar cooling can set a very strong bound on couplings to matter fields. For smaller m Z , the on-shell production of the transverse Z is suppressed but longitudinal Z can still be produced on-shell contributing to the cooling of the star. Considering the upper bound on the cooling rate, Refs. [30,31] find strong upper bound on the Z coupling to electrons from Z production in stars. This bound is shown with a blue line in fig. 1 which corresponds to |a e | = 1 (e.g., for gauged B − L).
Bounds from neutron scattering and meson decay: Another bound comes from the n−Pb scattering experiment [32] which is shown in fig. 1 with dashed-dotted orange line. For m Z > 30 MeV, the strongest bound comes from K + → π + +missing energy [33,34] but this bound is irrelevant for m Z <keV. Notice that in our model, Z decays to ν sνs appearing as missing energy so the recent searches for dark photon from NA48/2 which looks for π 0 → Z γ, Z → e − e + are not relevant for our model. Moreover, the NOMAD bound on π 0 → Z γ [36,37] comes from search for subsequent production of π 0 from Z scattering off nuclei in detector. The NOMAD bound does not apply to our case because in our model, Z decays to ν sνs before reaching the detector. From the three-body decay of charged mesons (K + (or π + ) → νZ µ + (or e + )), strong bounds on the coupling of Z to active neutrinos (g αβνα γ µ ν β Z µ ) can be derived α |g eα | 2 , α |g µα | 2 2 × 10 −9 (m Z /10 eV) [38][39][40] which in the context of B − L models also applies to g B and is shown in fig. 1.
Let us now discuss the bound on g s . Since no sterile neutrino has so far been discovered, it is no surprise that the bounds on g s are not very strong. However, from cosmological observations, the following constraints can be already set on g s .
Neutrino decay: Because of the neutrino mixing, the coupling of g s to ν s can lead to decay of the heavier neutrino mass eigenstates to lighter ones. If Z is heavier than ν 4 , the two-body decay will be forbidden but the three-body decay ν i → ν jνk ν l can take place where ν i can be any neutrino mass eigenstate other than the lightest one. The decay rate of ν i with energy E ν can be estimated as where the mass of the final neutrinos are neglected. The factor m i /E ν takes care of the time dilation because neutrinos decay in flight. Taking U s4 ∼ 1, U si | i∈{1,2,3} ∼ 0.1, m Z ∼ 10 eV and m 4 ∼ 1 eV, we find cτ 4 = cΓ −1 4 ∼ 5 × 10 31 cm(3 × 10 −4 /g s ) 4 (E ν /300 MeV). This means that ν 4 produced in the oscillation of the atmospheric neutrinos, in long (and, of course, short) baseline neutrino beams or in the way from the Sun cannot decay before reaching the detector. However, ν 4 produced in the oscillation of cosmic neutrinos can decay. Moreover, as long as g s 3 × 10 −4 , ν 4 produced in the early universe can decay before recombination. The decay of ν i (where i stands for 2 and 3 (1) for normal (inverted) ordering) will be further suppressed by (m 6 i /m 6 4 )U 2 si 1 so the decay cannot be relevant even for cosmic neutrinos.
Bounds from cosmology: We first discuss the bounds from m ν and free streaming of neutrinos at recombination. We then address the production of ν s before neutrino decoupling and the contribution from ν s to extra-relativistic degrees of freedom, (δN ν ) ef f .
As discussed above for g s 3 × 10 −4 , ν 4 decays into lighter neutrino states before the onset of structure formation so the bound on the sum of neutrino masses does not restrict our scenario (see also Ref. [41]). We should however check the bounds on extra relativistic degrees of freedom and on free streaming of neutrinos.
The self-interactions of ν s induce a nonzero V s in early universe which suppresses the active sterile mixing and therefore the sterile neutrino production before neutrino decoupling and BBN era. The bounds from BBN can be therefore relaxed if the g s coupling is larger than O(10 −3 −10 −2 ) and the mass of the gauge boson coupled to ν s is smaller than MeV [24]. The second condition can be readily fulfilled in our model. Let us check whether g s ∼ O(10 −3 ) − O(10 −2 ) can prevent neutrinos from free streaming at structure formation. As discussed before, for g s 3 × 10 −4 , the ν 4 component of the active neutrinos decays away. However, the ν i component also has a coupling of g s U 2 siν i γ µ ν i Z µ which leads to a self-interaction cross section of g 4 s |U si | 8 T 2 /(4πm 4 Z ) at T 2 m 2 Z . For g s U 2 si > 5 × 10 −6 , this means that neutrinos will not be free steaming at the onset of structure formation so for U si ∼ 0.1, g s ∼ O(10 −3 ) − O(10 −2 ) is ruled out by the free streaming of ν at recombination. The suppression of effective mixing at T >MeV therefore requires another mechanism which we will introduce later. Thus, we take g s ∼ 3 × 10 −4 to guarantee the free streaming of all active neutrinos as well as fast decay of ν 4 . Notice however that for T > m 4 process ν iνi → ν 4ν4 and ν i ν i → ν 4 ν 4 can take place with mean free path shorter by a factor of U 4 si . This sudden increase in the mean free path of ν i in the threshold of structure formation may leave an observable effect that requires more thorough investigation beyond the scope of the present paper. For a fixed mixing angle, suppressing ν i + requires smaller values of g s (less than 5 × 10 −5 ) but then ν 4 cannot decay to ν i | i<4 fast enough. Introducing lighter degrees of freedom with relatively large coupling to Z can make fast enough decay of ν 4 possible but we shall not explore this addition.
Within the scheme of the 3+1 solution to the LSND, ν s can be produced through oscillation in the early universe when (∆m 2 41 The produced ν s contribute to extra relativistic degrees of freedom on which there are strong bounds from BBN and CMB. Let us see whether the new neutral current (NC) interaction between sterile neutrinos and quarks can induce an effective potential leading to suppression of effective neutrino mixing. As is well-known at high temperatures, because of charged current (CC) interactions between (−) ν e and e ± , an effective potential proportional to G 2 F T 4 E ν emerges which is non-zero even when the densities of e − and e + are equal. Since the new interaction that we turn on between sterile neutrinos and quarks is of neutral current type, we do not expect a similar effect here. The matter effects on ν s need asymmetry between baryons and anti-baryons at T < GeV. Since at these temperatures, baryons are non-relativistic, we can write V s = √ 2G F (3n B ) where 3 reflects the fact that there are three valence quarks in each baryon. In order to have suppression of mixing due to NSI at temperature T , NSI should satisfy the bound Inserting n B = η B n γ ∼ 10 −10 n γ and T ∼ MeV, we find 10 9 . Considering the strong bounds on the coupling of nuclei to new particles it seems impossible to have such large . Moreover, such large is ruled out by neutrino oscillation experiments [42]. However within asymmetric dark matter (χ) scenario [43][44][45][46][47][48][49][50][51][52][53], we may be able to achieve desired suppression of active sterile mixing [18] invoking NSI of form where b is an arbitrary real number. For definiteness, let us take the benchmark point n χ = n B (and m χ /m p 5) that is motivated by the scenario within which the dark matter and baryonic matter asymmetries are simultaneously created by the same mechanism. The suppression of active sterile mixing then requires Such effective couplings shown in Eqs. (1.2,2.2) can originate from integrating out the intermediate U (1) gauge boson, Z . Let us denote the coupling of Z to quarks with g B /3 and those to ν s and χ with respectively g s and g χ . In general, χ and ν s can have different U (1) charges. A long as g χ < 10 −4 , the condition m χ α χ m Z is satisfied so we are in the perturbative regime [54]. From Bullet clusters, the bound σ(χχ → χχ)/m χ < 1.25 cm 2 /gr has been derived [55,56]. For non-chiral interaction with b = 0, we can write [57] so the bound σ(χχ → χχ)/m χ < 1 cm 2 /gr for χ with velocity of v χ ∼ 10 −3 in the galaxy implies g χ < 0.05 which is readily satisfied in the perturbative range. However, for b = 0 the cross section σ(χχ → χχ) The enhancement comes from the q µ q ν /m 2 Z part of the t-and u-channel propagator and the fact while q µχ2 (q + p)γ µ χ 1 (p) = 0, the axial part does not vanish: q µχ2 (q + p)γ µ γ 5 χ 1 (p) = m χχ2 (q + p)γ 5 χ 1 (p). To avoid too strong self-interaction, we set b = 0. That is we take the coupling to dark matter to be non-chiral.
In principle, χχ → χχZ can lead to dissipation. To prevent this, we should require where the factor of 10 takes care of the IR logarithmic enhancement. The m 2 Z dependence comes from the longitudinal component of Z . The dissipation constraint then leads to g χ < 0.008 which is again satisfied within the perturbative range.
Finally let us discuss the possibility that at temperatures higher than the electroweak symmetry breaking new processes involving χ lead to thermalization of ν s and Z . The densities of ν s and Z will be diluted because of the entropy dump in plasma so their contribution to N ef f will be negligible, avoiding the bounds from BBN and CMB. At temperatures MeV< T < 100 MeV where ∆m 2 21 /T H, the active neutrinos are mainly in form of incoherent (effective) mass eigenstate, ν i which have a coupling of formν i γ µ ν s Z µ given by g s U si [∆m 2 41 /(2E ν V s )]. The matter potential due to ADM is where η χ is the DM asymmetry parameter; η χ ≡ (n χ − nχ)/n γ . Since Ω χ 5Ω B , one typically requires m χ η χ = 5m p η B , suggesting m χ 5 GeV when η B = η χ . Let us check if the scattering of ν a off relic ν s can populate ν s at T > 1 MeV via ν a +ν s → ν s +ν s . The rate of scattering of each relic ν s can be written as Γ = n ν σ s where for T m Z , the process has cross section of σ s ∼ (g 4 which for values of g s and g χ giving rise to χ > 10 9 can be readily satisfied.
In summary, combining the conditions of decay of ν 4 before recombination and free streaming of light neutrinos at that epoch implies Moreover, the suppression of ν s (and therefore of Z production) before neutrino decoupling implies χ 10 9 . Bounds from neutrino oscillation experiments: As mentioned in the introduction, the MiniBooNE collaboration has suggested resonance matter effect with E res ∼ 300 MeV as a solution to the MiniBooNE anomaly. This value of resonance for energies and densities of interest to MiniBooNE corresponds ∼ 10 4 . Such large within our model requires g s ∼ 10 −3 and g B ∼ 10 −13 − 10 −12 at m Z ∼ 10 eV which are allowed. However, with ∼ 10 4 , the effective mass splitting between active neutrinos for E ν > 1 GeV obtains a correction given by ∆m 2 41 (1 − |U s4 | 2 ) [58] which is much larger than ∆m 2 31 . Thus, the pattern of oscillation in long baseline and atmospheric neutrino data will be dramatically changed, constraining to values smaller than O(10) [42].
Impact on supernova evolution: The coupling of Z particles to standard model particles is too small to allow an efficient production of Z in a supernova. If ν s particles exist inside the supernova, they can produce Z particles. However, inside the supernova the matter effects suppress the effective mixing between active and sterile neutrinos. Moreover, since the rate of scattering of active neutrinos off nuclei via electroweak interaction is much larger that oscillation time, active neutrinos stay in coherent active form with no coupling to Z so processes such a ν aνa → Z Z or ν sνs cannot take place 2 . Thus, there is no possibility of ν s production in the core. When the active neutrinos or antineutrinos stream out of the core and reach densities of O[10 4 gr/cm 3 (10/ )], depending on the sign of V s , they may undergo resonant conversion to sterile neutrinos which will leave its imprint in the spectrum of neutrinos observed on the earth. Studying the whole effect is beyond the scope of the present paper.

The Model
In this section, we discuss how to build a U X (1) gauge model with gauge boson of mass ∼ 10 eV that can give rise to sizeable non-standard interactions of sterile neutrinos with matter fields. We shall call this model the 3+1+U (1) model. Let us take the gauge coupling equal to g B and assign the following U X (1) charges to the standard model fermions: B + a e L e + a µ L µ + a τ L τ .
(3.1) Notice that we have assigned equal U X (1) charges to quarks of all three generations. Had we assigned unequal charges to quarks of different flavors, Z would have obtained couplings of form Z µq i γ µ q j where q i and q j are quarks of different masses. Such a coupling could lead to q i → q j Z enhanced by (m q i /m Z ) 2 . As long as a e + a µ + a τ = −3, the chiral anomalies cancel out. For a e + a µ + a τ = −3, new chiral fermions charged under U X (1) should be added to the standard model to cancel the anomalies.
The effective coupling between quarks and active neutrinos will be then of form In the parameter range of our interest, the NSI couplings of active neutrinos can be estimated as αα ∼ 10 −12 a α (g B /10 −16 ) 2 (10 eV/m Z ) 2 . The non-standard interactions of active neutrinos will be therefore too small to be observable. Remember that in our scenario, active neutrinos are not supposed to have sizeable NSI.
The sterile neutrino is taken to be an electroweak singlet fermion with U X (1) charge equal to a s = g s /g B . However, mixing between ν s and ν α breaks both SU (2) × U (1) and U X (1) which can be achieved by the mechanism described below [59]: Let us introduce a new singlet Weyl fermion N R neutral under both the SM and U X (1) and a new scalar φ which is charged only under U X (1) with a charge equal to that of ν s . We can then write the following Yukawa couplings where H is the SM Higgs. VEVs of H and φ respectively induce masses of form m aN = λ aN H and m νs = λ sN φ . Taking m aN /m νs ∼ U a4 ∼ 0.1 and m νs ∼ 1 eV, we obtain λ aN ∼ 10 −12 so the branching ratio of invisible decay mode H → νN is quite suppressed. Moreover this interaction at early universe cannot bring N R to thermal equilibrium λ 2 aN T H = T 2 /M * P l for T 200 GeV. The mass term can be written as where m a is the mass of active neutrinos, which might be produced by any of the mechanims introduced in the literature. Taking m 2 a m 2 aN m 2 νs , the mass matrix can be diagonalized via with the mass eigenvalues equal to {m a , −m νs , m νs }. Notice that up to corrections of O (m aN /m νs ) 3 , m a will not receive corrections. In other words, the contribution from "seesaw-like" mechanism to active neutrino mass is m νs (m aN /m νs ) 3 which is much smaller than ∆m 2 atm so m a should originate from another mechanism. Moreover, for the first approximation ν a does not mix with N R . The interesting point is that even after turning on V s , ν a does not mix with N c . This can be understood as follows. The Weyl fermions ν s and N R together form a Dirac fermion, ψ where ψ L = ν s and ψ R = N R . Regardless of whether we turn on matter effects, only ψ L will be involved in oscillation. This is a well-known result which comes from the fact the oscillation is given by m † ν .m ν rather than by m ν . The VEV of φ gives a contribution to m 2 Z given by g 2 s φ 2 . Taking m νs = λ sN φ ∼ 1 eV, the condition φ ≤ m Z /g s = 33 keV(m Z /10 eV)(3 × 10 −4 /g s ) implies λ sN ≥ 3 × 10 −5 . Taking the quartic self coupling of φ to be λ φ , it is natural that the mass of φ to be m φ ∼ λ φ φ < 33 keV . Taking λ φ ∼ 0.1 − 1, φ particles can be produced in the early universe at temperatures around 100 keV when ν s come into equilibrium with active neutrinos. At these temperatures the φ mass obtains a correction of λ sN T ∼ (m νs / φ )T which is smaller than φ and m φ for φ > 300 eV and therefore negligible. When the temperature drops below m φ , φ decays to sterile neutrinos, N and ν s which in turn decay into lighter neutrino states. As mentioned above, N and ν s form a Dirac fermion whose left-handed component has a U X (1) charge. This induces a U X (1) − U X (1) − U X (1) chiral anomaly. As discussed in the previous section, the coupling of the dark matter field, χ, has to be non-chiral so it cannot help with anomaly cancellation. To cancel anomaly, we may add another chiral field, ν R with U X (1) charge equal to those of ν s and φ. If the mass of ν R comes from a VEV of a scalar charged under U X (1), it cannot be heavier than ∼ 30 keV. ν R can have a small coupling of form λ RN φ † ν T R cN R with λ RN λ sN , providing it with a decay mode so in case it is produced in the early universe, it can decay away immediately.
As mentioned before, with a e + a µ + a τ = −3, chiral anomalies involving the SM fermions cancel without a need for extra fermionic degrees of freedom which are chiral under the gauge group. The B −L symmetry (i.e., taking a e = a µ = a τ = −1 ) is compatible with mixing in the lepton sector. That is obtaining a mixing between different neutrino flavors does not require an extra scalar whose VEV breaks U X (1). Of course, anomaly cancellation requires adding the right-handed neutrinos. These right-handed neutrinos can help to give Dirac mass to neutrinos. If no Majorana mass (which breaks U X (1)) is induced, these righthanded neutrinos will be degenerate with light active neutrinos and can kinematically be produced in the early universe. Taking g B 10 −13 , they can never come to thermal equilibrium: It is also possible to implement seesaw mechanism by introducing an extra scalar whose U X (1) charge is equal to −2 times U X (1) charge of leptons S ν T R cν R . The upper bound on S from the mass of Z is 100 TeV(m Z /10 eV)(10 −13 /g B ).
Since in the electrically neutral mediums such as the Earth, the number densities of the electrons and the proton are equal, for the case of gauging the B − L symmetry, the contributions from the electron and proton to V s cancel out and V s turns out to be proportional to the neutron number density: where the factor of 3 reflects the fact that there are 3 quarks in a neutron 3 . This factor of 3 compensates the factor of 1/3 in Eq. (2.5). Lastly note that if the VEV is zero until very late times (i.e. T < MeV), then the cosmological bounds are also significantly weakened. This is because the presence of a nonzero VEV is required in order for the active and sterile neutrinos to mass mix, and the only mechanism for thermalization for the steriles is via mixing. These scenarios have been explored in [23,60]. The parameter range of interest comes out naturally without a need for fine-tuning. In Table 1 we present a fiducial summary of parameters needed in the model consistent with all bounds.

The Oscillation Picture
Despite the 6 σ evidence for a ∼ 1 eV sterile neutrino from LSND and MiniBooNE along with other weaker hints, there is compelling evidence from IceCube, MINOS, MINOS+ against the 3 + 1 scenario with a ∼ 1 eV sterile neutrino. In this section we address whether the addition of a new interaction in the sterile sector modifies these constraints.
In the 3+1 scenario, the unitary mixing matrix U is the PMNS matrix [61,62] extended to include a fourth generation via U = R 34 (θ 34 )R 24 (θ 24 )R 14 (θ 14 )U PMNS . We choose this definition to make the U e4 = sin θ 14 and U µ4 = cos θ 14 sin θ 24 terms (which are relevant for LSND and MiniBooNE) as simple as possible. We have taken all new CP violating phases to be zero for simplicity. Section 4.1 describes the predictions of the 3 + 1 and 3 + 1 + U (1) scenarios for high energy atmospheric neutrinos at IceCube. Section 4.2 summarizes IceCube results for the 3 + 1 + U (1) scenario. Section 4.3 summarizes the bounds from MINOS and MINOS+ on the 3 + 1 + U (1) scenario. Section 4.4 discusses the bounds from IceCube and MINOS/+ allowing all the mixing angles including θ 34 to be nonzero.  Figure 2. The ν µ → ν µ survival probabilities at cos θ Z = −0.8 for three different cases. For the SM the standard oscillation parameters are taken from [65,66]. For the other cases we take our benchmark sterile parameters: θ 14 from reactor data and θ 24 from a global fit including LSND and MiniBooNE. The NSI case is evaluated at our best fit point from IceCube (see section 4.2 below) = −1.3.

The Atmospheric Neutrino Constraint from IceCube
IceCube has provided one of the strongest constraints on ∼ 1 eV sterile neutrinos by measuring atmospheric ν µ disappearance probabilities at energies 1 TeV [63]. Beyond these energies, the oscillation length becomes larger than the Earth diameter so neutrinos are not expected to oscillate in the Earth. Thus, at E 1 TeV without any new physics we have P (ν µ → ν µ ) ≈ 1 as shown in black in fig. 2. The existence of both a ∼ 1 eV sterile neutrino and NSI's will alter this as was shown in [14].
We calculate the effect of a sterile neutrino with a new interaction as it would affect IceCube's measurement. For the matter density profile we use the Preliminary Reference Earth Model [64] and calculate the oscillation probabilities through the Earth numerically. Phenomenologically, our model contains five new parameters beyond the standard oscillation parameters: ∆m 2 41 , θ 14 , θ 24 , θ 34 , and . For the standard oscillation parameters we use the result of the global fit and assume the normal mass ordering (which is preferred at ∼ 3 σ) from [65,66].
LSND and MiniBooNE interpret their results in a two flavor picture where the probability is given by This mixing angle is related to the angles given above by sin 2 2θ µe = sin 2 θ 24 sin 2 2θ 14 .
We fix sin 2 θ 14 = 0.095 which is the best fit point of the global fit to reactor ν e disappearance without input from theoretical reactor fluxes [67]. While a recent analysis from Daya Bay suggests that the Reactor Antineutrino Anomaly (RAA) may be due (in part or in full) to nuclear effects [68], a more recent article indicates that more analysis is required [69]. Regardless, even in the event that the RAA is due to nuclear effects we take 0.095 as a reasonable upper limit on sin 2 θ 14 . We want θ 14 large to minimize the amount  Figure 3. The 90% exclusion region from IceCube atmospheric neutrino data with 500 GeV < E ν < 10 TeV for 2 degrees of freedom is shown as a function of ∆m 2 41 and sin 2 2θ 24 where θ 14 is fixed to the best Reactor Antineutrino Anomaly value. Also shown is the 3 σ allowed region from the global fit to ν µ → ν e appearance data including LSND's decay-in-flight data; the best fit point is denoted with a star [67]. On the left is the sterile only model (i.e., 3 + 1), and on the right an additional minimization is performed over | | ≤ 1.7 for the 3 + 1 + U (1) scenario. of sterile mixing in the ν µ sector where the IceCube (and MINOS and MINOS+ below) constraints are strong.
For the remaining two sterile parameters, we define our benchmark values based on a global fit to the ν µ → ν e data including LSND's decay-in-flight data, which gives sin 2 2θ 24 = 0.0664 and ∆m 2 41 = 0.559 eV 2 [67]. While this does not include the latest MiniBooNE results, the new results do not change the region of interest much, only the significance. Figure 2 shows the oscillation probabilities for the SM, our benchmark sterile neutrino only model (i.e., 3 + 1), and our benchmark sterile neutrino model with NSI (i.e., 3 + 1 + U (1)) fixed to what will be our best fit result from IceCube. The large resonant oscillation probability for anti-neutrinos in the 3 + 1 model is the result of the MSW resonance [70,71] in the Earth, and is the signal that IceCube is particularly sensitive to. The addition of large NSI significantly suppresses the resonant conversion and mostly returns the probability to that in the standard 3ν scheme.
We use the publicly available two year upward-going muon neutrino flux from IceCube with ∼ 35, 000 events [72,73] with 501 GeV < E ν,p < 10 TeV and cos θ Z < 0 binned into 13 energy bins and 10 angular bins where E ν,p is the energy proxy used by IceCube. The atmospheric flux is provided by [74]. It is extended to the TeV range as described in [75] using the cosmic ray flux provided in [76] to account for the knee, and there is an estimated uncertainty of ∼ 25% at E ν ∼ 1 TeV. IceCube's effective area as a function of energy, angle, and flavor along with the Digital Optical Module (DOM) efficiency is provided by IceCube [73] and then used to calculate the expected number of events in each bin: where E ν is the true neutrino energy, the E ν,p and cos θ Z integrals are taken over the bin size, A is IceCube's effective area which includes DOM efficiency, Φ is the atmospheric flux, and P is the ν µ disappearance probability. We then construct a χ 2 test statistic, where the sum is over the E ν,p and cos θ Z bins, and we have accounted for the atmospheric flux normalization uncertainty σ x = 0.25 with a pull term [77].

IceCube Results
We perform a scan over the two sterile parameters θ 24 and ∆m 2 41 for the case with no NSI (i.e., the 3+1 scheme) and for the 3 + 1 + U (1) scenario where we minimize over | | ≤ 1.7. The results are shown in fig. 3. While the sterile neutrino picture is disfavored by the IceCube data, we find that the addition of NSI with | | ≤ 1.7 describes the data with 500 GeV< E ν < 10 TeV well.
Next, we fix the sterile parameters to their best fit values from [67] and vary only the NSI parameter in fig. 4. The 3 + 1 hypothesis is disfavored compared to the standard 3ν scheme at ∆χ 2 = 15.1. The addition of NSI in the sterile sector not only relaxes this tension, but provides a slightly better (∼ 1 σ) fit to the data with 500 GeV< E ν < 10 TeV for = −1.3 +0.2 −0.8 where the error bars are the 1 σ uncertainty. Finally, we plot the sum of both experiments in red. We find that the best fit point for both experiments for 3 + 1 + U (1) is at = −0.07 with ∆χ 2 = 26 compared to the SM. The 3 + 1 + U (1) is preferred over the best fit 3 + 1 point (at ∆χ 2 = 40) by an improvement of ∆χ 2 = 14.
Two additional features are of note in fig. 4. First, at = −1/3 the χ 2 increases by more than 200 due to an MSW resonance appearing in the neutrino channel at ∼ 700 GeV. Next at = −1/12 we see a sharp feature wherein the χ 2 nearly returns to its value within the standard 3ν scheme. This is because = −1/12 means that the sterile neutrino feels the same NC interaction that the active neutrinos do, and the only difference is due to the effect from vacuum oscillations which are relatively small at these large energies. While we were preparing this paper, Ref. [78] appeared which shows that for E ν < 500 GeV, NSI worsens the agreement with IceCube DeepCore atmospheric neutrino data. We note that while it would be possible to evade even the DeepCore bounds with very large NSI 20 which pushes the resonance below 10 GeV, at that point tight constraints from Super-KamiokaNDE [79] apply.

The Long-Baseline Constraint from MINOS and MINOS+
MINOS and MINOS+ also set strong constraints on the LSND and MiniBooNE best fit 3+1 point. These bounds are dominated by the MINOS+ CC analysis [28]. Using publicly available data and covariance matrices, we construct a χ 2 including neutrinos and antineutrinos, CC and NC, far and near detectors, and appearance and disappearance channels for both MINOS and MINOS+.  There is a small region where sterile neutrinos and NSI are actually slightly favored over the SM by IceCube with a best fit value of = −1.3 +0.2 −0.8 . For MINOS and MINOS+ the picture is never better than the SM and the improvement over the 3+1 picture is only marginal. The best fit point for both experiments is at = −0.07 at ∆χ 2 = 26 which is preferred over the best fit 3 + 1 (at ∆χ 2 = 40) by an improvement of ∆χ 2 = 14.
We find that, to some extent, it is possible to relax the MINOS and MINOS+ constraint on the 3+1 by including non-zero . The best fit 3+1 point is disfavored relative to the SM at ∆χ 2 = 24.9. Within the 3 + 1 + U (1) scheme, this can be slightly improved relative to the 3 + 1 case by ∆χ 2 = 1.9, although it is still disfavored compared to the SM at ∆χ 2 = 22.9. The best fit value of NSI is = 0.7 +0.4 −0.5 . See the right panel of fig. 5.

Complete Oscillation Picture
Setting θ 34 = 0, in section 4.3, we found that while the IceCube 3+1 picture for 501 GeV < E ν < 10 TeV is completely ameliorated by the addition of NSI, the MINOS and MINOS+ picture is only slightly improved. Moreover, the values of for each turn out to be quite distinct, with only moderate gains compared to the 3 + 1 picture when both experiments are considered simultaneously. We now check how the pictures change by allowing θ 34 to vary. We again fix θ 14 , θ 24 , and ∆m 2 41 as above and scan over and θ 34 . We find that the MINOS and MINOS+ picture can be considerably improved, but the agreement is still worse than that in the standard 3ν scheme. The best fit point is disfavored relative to the standard 3ν scheme at ∆χ 2 = 10.2 which is still better than the result without θ 34 of ∆χ 2 = 22.9. We perform the same analysis with the IceCube data and find that varying  Figure 5. The 1, 2, and 3 σ best fit regions (blue, green, and red respectively) at 2 d.o.f. in the , θ 34 space with the best fit point denoted by the star for IceCube (left) and MINOS and MINOS+ (right). While the best fit point for IceCube provides a good fit to the data and is slightly preferred over the SM, the best fit point for MINOS and MINOS+ is still disfavored compared to the SM by ∆χ 2 = 10.2. θ 34 from zero worsens the fit considerably for most values. Importantly, we did not find any point in the − θ 34 parameter space which provides suitable fit both for the MINOS and MINOS+ data and for the IceCube data at the same time. The 1, 2, and 3 σ contours for each case are shown in fig. 5.
We note that we have still kept the new CP violating phases fixed to zero and these could provide some additional relaxation, although their effect is expected to be small. A complete global fit varying all the sterile and NSI parameters is beyond the scope of the present paper.

Summary and discussion
In the dawn of neutrino precision era, studying possible effects of non-standard neutrino interaction with matter fields on long baseline and atmospheric neutrino experiments have received significant attention. From a theory point of view, it is challenging to build viable electroweak symmetric models that give rise to NSI with effective coupling large enough to have discernible effects in neutrino experiments. One possibility is to invoke light mediators (see also [80][81][82] for related work). It has been shown in the literature that invoking a U (1) gauge symmetry with a gauge boson (Z ) at the MeV-scale, phenomenologically relevant NSI can be obtained. For 10 keV < m Z <few MeV, the cosmological bounds on the extra relativistic degrees of freedom from the BBN and CMB constrain the maximal effective NSI coupling. In this paper, we have explored another mass window: m Z ∼ 10 eV where the bounds from cosmology are relaxed and the most stringent bound on the matter field coupling to Z (< few × 10 −11 ) comes from the fifth force searches. In this mass range, the strongest bound on neutrino coupling g αβνα γ µ ν β Z µ to Z come from Kaon decay K + → e + ν α Z , µ + ν α Z which constrains ( α |g eα | 2 ) 1/2 , ( α |g µα | 2 ) 1/2 10 −9 . In other words for m Z ∼ 10 eV, all g αβ except g τ τ are constrained to 10 −9 . Even with these stringent bounds, active neutrinos can obtain effective NSI coupling of order of G F or even larger.
In this paper, we have focused on building a model that provides NSI for sterile neutrinos with matter fields. We have called this scenario the 3 + 1 + U (1) model. The motivation for the model is the fact that by introducing NSI for the fourth neutrino of the 3 + 1 solution to LSND and MiniBooNE anomalies, the bound from the high energy IceCube atmospheric neutrino data can be relaxed (for energies greater than 500 GeV). Our analysis of the atmospheric IceCube data shows that within the 3 + 1 + U (1) model, not only the IceCube constraint is relaxed, the quality of fit is actually improved over the fit of standard 3ν scheme, however we urge the IceCube collaboration to publish a single sterile neutrino analysis spanning from DeepCore energies up through the highest available energies.
Obviously, to obtain an effective interaction between ν s and matter fields, the Z gauge boson has to couple to both ν s and matter fields. In our model, the coupling of Z to gauge boson comes from gauging the B − L symmetry so the contributions from electrons and protons of a medium to the effective neutrino mass cancel each other out rendering the matter effects proportional to the neutron number density and independent of the electron and proton density.
As expected, the ν s coupling to a new Z gauge boson is not restricted by terrestrial experiments. However, the requirement of free streaming of active neutrinos mixed with ν s at recombination from one side and the condition of the decay of ν 4 before that era from another side restrict g s to be around 3 × 10 −4 . With this value of g s and saturating the bound on the coupling to quarks, the four-Fermi effective coupling of ν s to matter fields, , can be much larger than 1 but for U si ∼ 0.1 (i = 1, 2, 3), such large values of NSI are ruled out by neutrino oscillation experiments themselves. An effective coupling of ∼ 1 which is favored for reconciling the 3 + 1 scheme with IceCube atmospheric neutrino data can be easily obtained with a value of coupling to matter fields five orders of magnitudes below the most stringent upper bound. The mixing of active neutrinos with ν s breaks both the new U X (1) and electroweak symmetry so it requires a new scalar field, φ, with a nonzero VEV. The requirement of lightness of Z then implies φ ∼ 30 keV. The U X (1) gauge coupling of SM fermions are too small to populate the new light states in the early universe but ν s and consequently the rest of light new state can be produced via ν a → ν s oscillation when T >MeV. To prevent such extra relativistic degrees of freedom in our model, Z also couples to asymmetric dark matter, χ. The forward scattering of ν s off the dark matter background for T >MeV induces a large effective mass for ν s suppressing the effective active sterile mixing and therefore preventing a new contribution to relativistic degrees of freedom at the BBN era. At T ∼ 100 keV (well below active neutrino decoupling from plasma) the active sterile mixing resumes its vacuum value and ν a → ν s can take place. The produced ν s can in turn populate other light degrees of freedom like φ and Z but they decay back to active neutrinos at the onset of structure formation. We found that if Z coupling to dark matter is chiral, the self-interaction of dark matter particles with such light Z will be enhanced by m 4 χ /m 4 Z violating the present bound. We therefore take the coupling to χ to be non-chiral. We carried out a numerical study to test the ability of such a model in reconciling the conflicting evidence on sterile neutrinos. The addition of a matter effect for the sterile sector provides no change to the short baseline experiments that provide evidence for the sterile neutrino. As expected, introducing NSI with an effective coupling of the same order of magnitude as G F completely removes the constraint from IceCube atmospheric neutrino data with 500 GeV < E ν < 10 TeV by returning the oscillation probability to one as it is in the Standard Model. The next most important constraint on sterile neutrinos comes from MINOS and MINOS+. We found that including NSI and turning on all three new mixing angles does provide moderate improvement to the quality of fit from the 3+1 case, it is never as good a fit to the data as the standard 3ν scheme. We examined whether nonzero θ 34 can help to relax the bound and found that the value of θ 34 necessary for relaxing the MINOS and MINOS+ bounds is not compatible with that from IceCube. Ref. [78] which appeared when this paper approached conclusion shows that including atmospheric neutrino data with E ν < 500 GeV from IceCube DeepCore [83] disfavors the 3 + 1 + U (1) solution to LSND. We eagerly anticipate more data from IceCube as their analysis progresses to draw conclusive verdict on the 3 + 1 + U (1) solution to the LSND and MiniBooNE anomalies.
Recently [84] has shown that reconciling MINOS+ bounds with the 3 + 1 solution of MiniBooNE and LSND requires non-standard charged current interaction of ν µ and µ with quarks with an effective coupling equal to 0.03G F . Such an effective coupling requires a charged mediator on which LHC has set strong bound. For example, if the mediator is a sequential W (W ), the bound from the LHC implies that the non-standard charged current effective coupling cannot be larger than 2 × 10 −4 . One should also bear in mind that the bounds from MINOS+ have been put into question by [29].
The next generation short baseline experiment, SBN, is being designed to test Mini-BooNE results [85]. At such short baselines SBN will not be sensitive to matter effects so the signature of our 3 + 1 + U (1) scenario will be similar to the 3 + 1 scheme. If SBN finds null results and the MiniBooNE anomaly is proved to be a systematic error (for example, due to underestimate of π 0 events), an upper bound on |U e4 | 2 |U µ4 | 2 will be set but the existence of ν s mixed with the active neutrinos will still remain an intriguing possibility which may be the source of other observed phenomena. Thus, the possibility of activating this fourth neutrinos via the 3 + 1 + U (1) scenario introduced in the present paper can be noteworthy independent of MiniBooNE and LSND data. For example, the ANITA experiment has observed two upward-pointing events with energies in the EeV∼ 10 18 eV range [86,87]. According to Ref. [88], these events can be interpreted as a sign for oscillation of active neutrinos to sterile neutrinos. The additional matter effects felt by a sterile neutrino in our 3 + 1 + U (1) scenario could have interesting implications for these ANITA events, a question we defer to future work.