Cornering (3+1) sterile neutrino schemes

Using the most recent atmospheric neutrino data, as well as short-baseline, long-baseline and tritium $\beta$-decay data we show that the joint interpretation of the LSND, solar and atmospheric neutrino anomalies in (3+1) sterile neutrino schemes is severely disfavored, in contrast to the theoretically favored (2+2) schemes.


Introduction
Reconciling the existing data on solar [1] and atmospheric [2,3] neutrinos with a possible hint at the LSND experiment [4,5] (indicating the existence of (−) ν µ → (−) ν e transitions) is a challenge to the simplest standard model picture. In the absence of exotic mechanisms and/or new neutrino interactions, such as neutrino transition magnetic moment [6], one requires the existence of neutrino oscillations involving three different scales. As a result a joint explanation for all the data (including the LSND anomaly) requires a fourth light neutrino which, in view of the LEP results, must be sterile [7,8] 1 .
There have been several theoretical models and phenomenological studies of 4neutrino models [9][10][11]. Two very different classes of 4-neutrino mass spectra can be identified: the first class contains four types and consists of spectra where three neutrino masses are clustered together, whereas the fourth mass is separated from the cluster by the mass gap needed to reproduce the LSND result; the second class has two types where one pair of nearly degenerate masses is separated by the LSND gap from the two lightest neutrinos. These two classes will be referred to as (3+1) and (2+2) neutrino mass spectra, respectively [12]. All possible 4-neutrino mass spectra are shown in Fig. 1.
Theoretically the existence of a light sterile neutrino sets a challenge. One possibility is to postulate a protecting symmetry [7] such as lepton number. Alternatively, in models based on extra dimensions, one may appeal to a volume suppression factor [13] in order to account for the light sterile neutrino. The theoretical origin of the splittings depends on the model. In Ref. [11] R-parity violating interactions are used, while in the original proposals the splittings were due to calculable two-loop effects. Such theories lead to a (2+2) symmetric scheme where two of the neutrinos combine to form a quasi-Dirac [14] or pseudo-Dirac [15] neutrino, whose splitting accounts for atmospheric oscillations, while the oscillations among the two low-lying states explain the solar neutrino data, with the overall scale accounting for the LSND result.
Although less motivated theoretically, it has been argued recently that (3+1) schemes are not strictly ruled out phenomenologically [12,16,17]. However, using the most recent atmospheric neutrino data, as well as short-baseline and tritium data we show that such schemes are severely disfavored. In order to accomplish this we extend the analysis of neutrino oscillation data in the framework of (3+1) neutrino mass spectra performed in Ref. [18]. In addition to the full data from the short-baseline (SBL) experiments Bugey [19], CDHS [20], KARMEN [21] and the result of the long-baseline reactor experiment CHOOZ [22], we also include the full and updated data set of atmospheric neutrino experiments and the data from the ν µ → ν e oscillation search in NOMAD [23]. With this information we derive a bound on the LSND amplitude A µ;e within a Bayesian statistical framework. We find that the inclusion of the full atmospheric zenith-angle distribution data considerably strengthens the bound on the LSND amplitude A µ;e for low ∆m 2 , whereas the NOMAD data strengthens the bound for high ∆m 2 values. In contrast, the details of the solar data do not play a key role in the analysis, so that there is no need to perform a full-fledged global fit of solar data. Likewise the most recent SNO results [24] will hardly affect our results. Finally, we also perform a different statistical analysis of the data in order to also include information from the tritium β-decay experiments [25,26]. This sets additional strong bounds on 4-neutrino spectra of the type (3+1)B (see Fig. 1).
In contrast to the (3+1) schemes, the intrepretation of the LSND anomaly in terms of (2+2) spectra is in good agreement with SBL data [18,[27][28][29][30]. It is a general prediction of (2+2) spectra that the sterile neutrino must take part either in solar or atmospheric neutrino oscillations (or both [11]). Atmospheric neutrino data prefer ν µ → ν τ oscillations over oscillations into sterile neutrinos (see e.g. Refs. [31,32]). Also fits to solar neutrino data in terms of active neutrino oscillations are better than sterile neutrino oscillation fits (see e.g. Refs. [33,34]), especially after the first SNO results [24]. However, a joint analysis of both solar and atmospheric neutrino data in a 4-neutrino framework shows that an acceptable fit can be obtained for the (2+2) schemes [35]. This paper is organized as follows. In Sec. 2 we fix our notations. In Sec. 3 the fit of atmospheric neutrino data in the framework of (3+1) mass spectra and its implications for parameters relevant in short-baseline oscillation experiments are discussed. In Sec. 4 we present the bound on the LSND amplitude A µ;e whereas in Sec. 5 implications of tritium β-decay experiments are discussed. In Sec. 6 we draw our conclusions.

Notation
The Standard Model can be extended with an arbitrary number of singlet Majorana leptons, as they carry no gauge anomalies [36]. The minimal case is to have simply one single neutrino [37] which we assume to remain light (due, for example, to some symmetry) and therefore able to take part in the oscillations 2 . In any such 4-neutrino gauge scheme the charged current weak interaction is characterized by the lepton mixing matrix K αj . This is a rectangular 3 × 4 matrix arising from the unitary 4 × 4 neutrino mixing matrix diagonalizing the neutrino mass matrix and the corresponding unitary 3 × 3 matrix diagonalizing the left-handed charged leptons. This matrix K αj contains in general six mixing angles and six CP phases [36].
All neutrino oscillation probabilities in vacuo are determined by the structure of the matrix K αj . For the case of solar and atmospheric neutrinos matter effects in the solar and/or Earth interiors must also be taken into account.
The probability of SBL ν e transitions relevant for the accelerator experiments LSND, KARMEN and NOMAD is given by a very simple twoneutrino-like formula [27] where L is the distance between source and detector and E is the neutrino energy and the parameters d α are defined as Note that for SBL oscillations we can safely neglect solar and atmospheric splittings relative to the LSND gap. With our labeling of the neutrino masses indicated in Fig. 1 the mass separated by the LSND gap is denoted by m 4 . This is the heaviest mass in spectra of type (3+1)A and the lightest in (3+1)B spectra. As a result ∆m 2 SBL ≡ ∆m 2 ≈ |m 2 4 − m 2 1 | in all cases. The LSND experiment gives an allowed region in the (∆m 2 , A µ;e ) plane.
The survival probabilities relevant in the SBL disappearance experiments Bugey and CDHS are given by where α = e refers to the Bugey and α = µ to the CDHS experiment. The result of the Bugey experiment requires d e to be very small or very close to 1. One can show that for the (3+1) spectra the survival probability of solar neutrinos is bounded by P ⊙ νe→νe ≥ d 2 e [27], so that d e must be small, and we can include this information from solar neutrinos in the analysis through the approximation d e (1 − d e ) ≈ d e [18].

Atmospheric data and short-baseline oscillations
We now discuss the implications of atmospheric neutrino data on SBL oscillation parameters. In Ref. [29] it was shown that the up-down asymmetry of atmospheric multi-GeV µ-like events measured in the Super-Kamiokande experiment [2] can be used to constrain the parameter d µ to values smaller than around 0.5. In the following we will see that a detailed fit to the full atmospheric neutrino data gives a much stronger constraint on d µ . For the analysis of atmospheric neutrinos we use the latest experimental data used in Ref. [35]: e-like and µ-like data samples of sub-and multi-GeV and upgoing muon data including the stopping and through-going muon fluxes from Super-Kamiokande [32] and the latest MACRO [40] up-going muon samples. For further details see Refs. [34,35,41].
The analysis of atmospheric neutrino data presented in Ref. [35] was performed for the (2+2) spectra adopting the approximations ∆m 2 ⊙ = 0 and ∆m 2 LSND → ∞. Moreover, in that analysis the electron neutrino was considered as completely decoupled from the atmospheric oscillations. This approximation is well justified for (2+2) spectra, because in this case the projection of ν e over the atmospheric states is severely restricted by the very strong Bugey bound. In contrast, in (3+1) spectra the contribution of electron neutrinos to atmospheric oscillations is limited only by the somewhat weaker CHOOZ bound. However, in Ref. [41] it was found that a ν e contamination small enough not to spoil the results of the CHOOZ experiment has only a very small effect on the quality of the fit of atmospheric neutrino data. Therefore, even in the context of (3+1) schemes it is reasonable to assume that electron neutrinos are decoupled, so that atmospheric oscillations actually involve only three neutrino flavours (ν µ , ν τ and ν s ). Under these approximations, the (2+2) and (3+1) spectra become identical (except for irrelevant signs), and therefore it is possible to use the analysis given in Ref. [35] also in the context of (3+1) schemes. In particular, the parameter d µ used in the present work corresponds to the parameter s 2 23 = |U µ1 | 2 + |U µ2 | 2 of Ref. [35].
In the left panel of Fig. 2 we show ∆χ 2 atm = χ 2 atm − (χ 2 atm ) min from the fit to atmospheric neutrino data as a function of the parameter d µ . For each value of d µ the χ 2 is minimized with respect to all other undisplayed parameters necessary to fit the atmospheric neutrino data. In the right panel of Fig. 2 we show the 90% and 99% CL bounds on d µ obtained from a combination of all the atmospheric neutrino data with the ν µ disappearance experiment CDHS in a Bayesian framework 3 .
In the lower part of this plot the constraint on d µ comes from atmospheric neutrino data alone, as the CDHS bound disappears for ∆m 2 0.3 eV 2 . Hence, atmospheric neutrino data leads to the bound d µ ≤ 0.090 (0.13) at 90% (99%) CL .
Let us note that Fig. 2 and the bound given in Eq. (4) are valid also for the (2+2) spectra if the parameter d µ is identified with |K µ1 | 2 + |K µ2 | 2 and neutrino masses are labeled according to Fig. 1 (right).

An upper bound on A µ;e
In this section we discuss the upper bound on the LSND amplitude A µ;e obtained by combining the data of the SBL experiments Bugey, CDHS, KAR-MEN and NOMAD, with those of the atmospheric neutrino experiments and the CHOOZ experiment. Here we focus mainly on the results of the extendend analysis, technical details can be found in Ref. [18].
To combine all the oscillation data we use the likelihood function The likelihood functions L Bugey , L CDHS and L KARMEN are described in Ref. [18] and L atm ∝ exp − 1 2 χ 2 atm [42]. To calculate L NOMAD we perform a re-analysis of the ν µ → ν e oscillation search at NOMAD by using the data given in Ref. [23]. The result of the CHOOZ experiment is included via L CHOOZ , which is obtained with the maximum likelihood method described also in Ref. [18].
For a fixed value of ∆m 2 the likelihood function Eq. (5) is transformed into a probability distribution p(d e , d µ ) by applying Bayes' theorem (see for example Refs. [42,43]) and assuming a flat prior in the physically allowed region d e , d µ ≥ 0 and d e + d µ ≤ 1. Choosing a CL β, we find the corresponding upper bound 3 An analysis using a χ 2 -cut method gives very similar results.
The bounds at 95% and 99% CL are shown in Fig. 3 together with the regions allowed by LSND at 90% and 99% CL. We find that there is no overlap of the region allowed by our bound at 95% CL with the region allowed by LSND at 99% CL [18]. If we take our bound at 99% CL there are marginal overlaps with the 99% CL LSND allowed region at ∆m 2 ∼ 0.9 and 2 eV 2 , and a very marginal overlap region still exists around 6 eV 2 . The overlap region found in Ref. [18] between 0.25 and 0.4 eV 2 is now excluded by our bound at 99% CL because of the inclusion of the full atmospheric neutrino data set 4 . 4 We note that for ∆m 2 10 eV 2 there are additional constraints on the amplitude A µ;e from the experiments BNL E776 [44] and CCFR [45], which are not included in our analysis. Therefore, the (in any case very marginal) overlap region at ∆m 2 ∼ 10 eV 2 in Fig. 3 is irrelevant.
For the alternative case of symmetric (2+2) spectra the bound on the LSND amplitude A µ;e is dominated by the Bugey bound on d e and by the KARMEN bound on A µ;e . Hence, an improvement of the bound on d µ by the inclusion of the full atmospheric neutrino data does not lead to a stronger bound on A µ;e in this case so that the results of Ref. [18] apply. The conclusion is that (2+2) spectra are in good agreement with SBL data (see Fig. 2 of Ref. [18]). As shown in [35] they are also in agreement with the totality of solar and atmospheric neutrino data.

Implications of tritium β-decay
Experiments studying the electron spectrum dN/dE e from tritium β-decay can obtain information on the quantity m 2 β determined by the relation where E e is the energy of the electron and E 0 is the total decay energy. The latest result obtained by the Troitsk experiment is [25] The Mainz collaboration recently presented two values, obtained from different analyses of their data [26]: In Eqs. (8)-(10) the upper bounds are at 95% CL.
Let us now consider the implications of these measurements for the (3+1)A and (3+1)B neutrino mass schemes (see Fig. 1). In the presence of neutrino mixing relation (7) has to be modified to (see e.g. [46]) In view of the very strong constraint on d e from Bugey we can safely neglect the contribution from |K e4 | 2 to the sum in Eq. (11). Further we take into account that mass splittings implied by solar and atmospheric neutrino data cannot be resolved in tritium decay experiments [46]. Hence, in spectra (3+1)A the value m β is given by the lowest neutrino mass and is independent of ∆m 2 to a good approximation. Therefore, we do not include any information from tritium β-decay in this case. However, for spectra of the type (3+1)B one obtains [47] m 2 β ≈ m 2 4 + ∆m 2 . 5 We include this result in our statistical analysis by using the likelihood function where the sum is over the three experimental values of m 2 β given in Eqs. (8)-(10) and σ i is the corresponding error (statistical and systematic errors added in quadrature). For fixed values of m 4 we perform now an analysis with the two parameters A µ;e and ∆m 2 , in contrast to Sec. 4, where the analysis is performed only with the parameter A µ;e for each value of ∆m 2 .
As a first step the total likelihood function obtained from Eqs. (5) and (12) is transformed into a probability distribution p(d e , d µ , log ∆m 2 ) by using Bayes' theorem. We assume a flat prior distribution for d e and d µ in the physical region and a flat prior distribution for log ∆m 2 . This ensures that we introduce no bias concerning the order of magnitude of ∆m 2 , a priori all scales are equally probable. Then we perform a transformation of the variables d e and d µ to 6 and integrate over the variable t to obtain the probability distribution for the variables we are interested in: p(A µ;e , log ∆m 2 ). We calculate an allowed region at the 100β% CL by demanding dA µ;e d(log ∆m 2 ) p(A µ;e , log ∆m 2 ) = β .
The boundary in the (A µ;e , ∆m 2 ) plane is determined such that the value of p(A µ;e , log ∆m 2 ) along this line is constant.
In Fig. 4 we show the allowed regions at 99% CL for spectra of the types (3+1). In the case of (3+1)B we show the regions for m 4 = 0, 1 and 2 eV. As we do not know the true value of m 4 the curve corresponding to vanishing m 4 is the most conservative one. In this case tritium β-decay rules out values of ∆m   The following remarks are in order: (1) To normalize the probability distribution p(A µ;e , log ∆m 2 ) one has to choose a lower integration bound for log ∆m 2 as p(A µ;e , log ∆m 2 ) does not vanish for log ∆m 2 → −∞. The allowed regions in Fig. 4 depend somewhat on this lower bound. However, if one chooses the lower bound sufficiently small (log ∆m 2 −3) the allowed regions become independent of it. In the case of the spectra (3+1)A the allowed region depends also on the upper integration bound for log ∆m 2 . However, again the dependence disappears, if this bound is chosen large enough (log ∆m 2 2).
(2) Note that the statistical meaning of the bounds in Figs. 3 and 4 is different. The method applied in Sec. 4 to produce Fig. 3 allows to place an upper bound on A µ;e for a given value of ∆m 2 at a certain CL, whereas the meaning of the 99% CL bounds shown in Fig. 4 is the following: the true values of A µ;e and of ∆m 2 are expected to lie at the left of the curves shown in the figure with probability 0.99. This explains the small differ-ences between the curve for the (3+1)A spectra in Fig. 4 and the 99% CL curve in Fig. 3. However, the general agreement of both methods renders confidence to our analysis. Finally, we note that the non-observation of neutrinoless double β-decay may also place important restrictions on (3+1)B spectra [47], where the electron neutrino has a substantial component along the heaviest neutrinos. However, this will be very model dependent as the resulting bounds are subject to possible destructive interference due to cancellations among different neutrinos [14,48].

Conclusions
We have extended the analysis of neutrino oscillation data in the framework of (3+1) neutrino mass spectra performed in Ref. [18]. In addition to the full data from the short-baseline experiments Bugey [19], CDHS [20], KARMEN [21] and the result of the long-baseline reactor experiment CHOOZ [22], we have included also the full and updated data set of atmospheric neutrino experiments, the data from the ν µ → ν e oscillation search in NOMAD [23] and tritium β-decay data [25,26]. We have shown that the interpretation of the LSND anomaly within such schemes is severely disfavored by the combined data, in contrast to the case of the theoretically preferred (2+2) schemes. Since the details of the solar data do not play an important role in our (3+1) analysis, the most recent SNO results [24] will not affect the conclusions derived in our paper, nor contribute to enhance the bounds we have derived. Our present results put an additional challenge on some recent attempts to revive (3+1) schemes [12,16,17].