Non-unitary evolution of neutrinos in matter and the leptonic unitarity test

We present a comprehensive study of the three-active plus N sterile neutrino model as a framework for constraining leptonic unitarity violation induced at energy scales much lower than the electroweak scale. We formulate a perturbation theory with expansion in small unitarity violating matrix element W while keeping (non-W suppressed) matter effect to all orders. We show that under the same condition of sterile state masses 0.1 eV2 ≲ mJ2 ≲ (1–10) GeV2 as in vacuum, assuming typical accelerator based long-baseline neutrino oscillation experiment, one can derive a very simple form of the oscillation probability which consists only of zeroth-order terms with the unique exception of probability leaking term C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{C} $$\end{document}αβ of O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(W4). We argue, based on our explicit computation to fourth-order in W, that all the other terms are negligibly small after taking into account the suppression due to the mass condition for sterile states, rendering the oscillation probability sterile-sector model independent. Then, we identify a limited energy region in which this suppression is evaded and the effects of order W2 corrections may be observable. Its detection would provide another way, in addition to detecting C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{C} $$\end{document}αβ, to distinguish between low-scale and high-scale unitarity violation. We also solve analytically the zeroth-order system in matter with uniform density to provide a basis for numerical evaluation of non-unitary neutrino evolution.


JHEP02(2019)015
The scenario of high-scale unitarity violation, based on the orthodox view of new physics at high-energy scales, has been studied extensively in the literature [21][22][23][24][25][26][27][28][29]. 3 On the other hand, there exist a good amount of activities in recent years hinting the possibilities of new physics at low energies. For scenarios which involve light sterile neutrino(s), see, e.g., [40][41][42], and the references therein. In the former case, due to the preserved SU(2) × U(1) symmetry at high scales, it is conceivable that the constraints from measurements using probes in the charged lepton sector play a dominant role. On the other hand, in the case of low-scale unitarity violation, neutrino oscillation experiments will play key role in constraining unitarity violation.
In a previous paper [20], we have proposed a model-independent framework for testing low-scale unitarity violation. It is based on the three active plus N sterile lepton (called neutrino) system, which is unitary in the whole (3 + N ) dimensional state space but restriction to observables in the active neutrino subspace renders the theory non-unitary in that subspace. It is referred to as the "(3 + N ) space unitary model". We have shown in the context of accelerator and reactor neutrino measurement that the restriction on the masses of sterile states 4 to 0.1 eV 2 m 2 J 1 MeV 2 (with J being sterile state index) is sufficient to make the observables sterile-sector model independent. That is, the neutrino oscillation probability can be written in such a way that it is independent of details of the sterile neutrino mass spectrum and mixing with active neutrinos. The model-independent nature of the framework will be translated into that of the constraints obtained, thereby making leptonic unitarity test more powerful.
As an outcome of our formulation we have pointed out a new way of distinguishing lowscale unitarity violation from high-scale one by observing the probability leaking term in the oscillation probability. The term signals existence of energetically accessible sterile states, which is characteristic to low-scale unitarity violation, and it has been included for the first time in the analysis of unitarity violation in [20] which uses a JUNO [43]-like setting. See refs. [28,44] for a comprehensive analysis of the currently available neutrino data with the active plus sterile framework, and [28,45] for analyses of the future experiments.
In this paper, we give a comprehensive treatment of the (3 + N ) space unitary model. We formulate a novel perturbative framework with small unitarity violating matrix element W as the unique perturbing parameter, which we call "small unitarity-violation perturbation theory". It allows us to calculate the oscillation probability in the presence of matter effect comparable in size to the vacuum mixing effect. It must be remarked that the sterile sector model-independent nature of the (3 + N ) space unitary model is demonstrated in ref. [20] only in vacuum and in matter to first order in matter perturbation theory. Hence, the first goal of this paper is to show that the model-independence holds after inclusion of sizeable matter effect. In fact, we observe that the same condition on the sterile neutrino masses guarantees this property, and the extremely simple expressions of the oscillation probabilities result even with our computation to fourth order. 3 Works have also been done on unitarity violation by sterile sector from somewhat different point of view, e.g., if it exists, how it could disturb measurement of lepton Kobayashi-Maskawa phase δCP, or mass ordering. See for example, [30][31][32][33][34][35][36][37][38][39]. 4 To be more precise, "masses of sterile states" implies hereafter masses of neutrinos which are mostly sterile. In this paper, for brevity, we use this simplified terminology in most of the places.

JHEP02(2019)015
The second goal of this paper is to utilize the oscillation probability formulas to uncover in which region of energies and baselines unitarity violating effect is large, and to examine the possibility of sizeable W 2 corrections which distinguishes between high-and low-scale unitarity violation. These exercises may be useful in the application of our framework to some of the ongoing and next generation neutrino oscillation experiments [16,17,[46][47][48][49][50][51][52]. To carry it out, we derive an exact expression of the oscillation probability in leading order in perturbation theory for uniform matter density. In summary, the framework can be used in dual modes: it serves (1) as a suitable framework for leptonic unitarity test in neutrino oscillation experiments, and (2) as a hunting tool for unitarity violation effects, which could serve for another way of distinguishing low-scale unitarity violation from high-scale one.

Essence of the present and the previous papers
In this section, we present essence of the present and the previous [20] papers, in which an adequate formulation is given to describe neutrino oscillations with unitarity violation caused by new physics at energies much lower than m W . In section 2.1 we define the system, section 2.2 serves for reviewing the content of ref. [20], and section 2.3 is to summarize the key points of this paper.

Unitary 3 active + N sterile neutrino system with partial decoherence
The system we are considering consists of 3 active + N sterile neutrinos which is unitary in the whole state space, but serves for a model of non-unitarity when restricted to observables in the active neutrino subspace. The sterile-sector model independence is realized due to decoherence between active-sterile and sterile-sterile states, which essentially wipes out detailed informations of sterile sector such as mass spectrum and mixing structure with active neutrinos. Generically, the decoherence condition associated with energy resolution reads (see [20]) 1 MeV 2 , or (1 − 10) GeV 2 , but otherwise the readers must assume that m J obeys the general conditions above.

Non-unitary evolution of neutrinos in vacuum or with small matter effect
Here we summarize the main findings of ref. [20]. Thanks to partial decoherence, fast oscillations in active-sterile and sterile-sterile channels are averaged out, which leads to a very simple form of the active neutrino oscillation probability in vacuum where x denotes baseline and in the appearance (α = β) as well as in the disappearance (α = β) channels with α, β = e, µ, τ . In eq. (2.2), the indices i, j, k = 1, 2, 3 and J = 4, 5, · · ·, N + 3 are, respectively, for (mostly) active and (mostly) sterile neutrino mass eigenstates. The active neutrino flavour states ν α are connected to mass eigenstates (ν i , ν J ) through where m j and m J denote the active and sterile neutrino masses, respectively, and E denotes neutrino energy.
The characteristic features of the oscillation probability in (2.2) are: 1. The non-unitary matrix U replaces the standard unitary three-flavour mixing matrix often parametrized with Particle Data Group convention U PDG [19].
2. Probability leakage term C αβ > 0 appears reflecting the nature of low-energy unitarity violation in which the probability can flow out from active neutrino space to the sterile state space, and vice versa.
3. Due to non-unitarity of the U matrix, δ αβ term in the unitary case is modified to

JHEP02(2019)015
Notice that each term of C αβ in (2.3) allows interpretation that "probability leaking from active to sterile state spaces" and coming back. The simple terminology of probability leaking assumes that the latter process must also exist which is ensured by generalized T invariance. Another aspect of the probability leaking term, which has a form of incoherent sum of the products of probabilities of active-to-sterile and sterile-to-active transitions clearly illustrates the decoherence associated to the sterile states. For instance, near the upper end of the sterile state mass region quoted in section 2.1, it describes effect of decoherence caused by separation of wave packets between active and sterile neutrinos. The points 2 and 3 above are important ones and the clarifying remarks about them are in order: • Presence or absence of the probability leakage term C αβ distinguishes between lowenergy and high-energy unitarity violation [20]. Nevertheless, C αβ may be small because it is of fourth order in W .
• Difference in normalization factor, the second term in (2.2), between unitary and nonunitary cases is of order ∼ W 4 (∼ W 2 ) in the appearance (disappearance) channels.
To understand the latter point, we notice that unitarity in the (3 + N ) space unitary model can be written as in the appearance channel (α = β), and in the disappearance channel (α = β), which justifies the above statement.
We emphasize, therefore, that the probability leaking term C αβ and the another con- in the oscillation probabilities are the same order, O(W 4 ), in the appearance channels. Hence, we do not see any good reasons why the former can be ignored, as was done in the existing literatures. It is also worth to note that O(W 2 ) difference in normalization in the disappearance channel would make detection of unitarity violation more feasible. It is one of the reasons for high sensitivity to unitarity violation that could be reached in disappearance measurement in the JUNO-like setting [20].
In the same work, by including small matter effect up to first order, we have found that as far as we remain in the region of unitarity violating element |W | 0.1, 6 or somewhat larger, the matter effect does not alter the above features of the oscillation probability in (2.2) under the same restriction on sterile neutrino masses. Notice that |W | 0.1 implies that the unitarity violating effect in the probability is of the order of |W | 4 ∼ 10 −4 , except for the O(W 2 ) difference in normalization constant in the disappearance probability. It is practically the limit of order of magnitude that can be explored by the next generation neutrino oscillation experiments.

Non-unitary evolution of neutrinos in matter to all orders
Given the fact that setup of some of the next generation accelerator LBL experiments require consideration of the matter effect comparable with the vacuum mixing one, it is clear that a better treatment is necessary to understand the influence of the matter effect in the (3 + N ) model. Then, we formulate in this paper the small unitarity-violation perturbation theory, a systematic and controlled way of treating small unitarity violation effect while including all order matter effect. We derive a simple expression of the oscillation probability in matter which retains the favourable feature of the vacuum formula (2.2), the sterile sector model independence under the same sterile neutrino mass condition as in vacuum. That is, the model-dependent terms are either averaged out, or made small due to large sterile state mass denominator suppression. We must note here that our treatment of the matter effect in this and the previous papers is restricted to the case of uniform matter density.
The resulting oscillation probability in matter between active flavour neutrinos in the (3 + N ) space unitary model to fourth order in W can be written as where h i (i = 1, 2, 3) denote the energy eigenvalues of zeroth-order states of active neutrinos in matter, and X is the unitary matrix which diagonalizes the zeroth-order Hamiltonian used to formulate our perturbation theory. C αβ is the same as we have in the vacuum case in (2.3). The expression is valid under the same restriction on sterile neutrino masses we have in vacuum, 0.1 eV 2 m 2 J (1 − 10) GeV 2 for |W | 4 10 −4 assuming neutrino energy and baseline (and the associated matter density) which correspond to accelerator LBL experiments. For more precise conditions we require and for the restriction needed on the sterile state masses for smaller W , see section 3.5.
The expression (2.6) is a very transparent result in the sense that (1) the vacuum non-unitary mixing matrix U is "dressed" in a simple way by the matter effect represented by X, and (2) the probability leaking term C αβ and the normalization term stay as they are in vacuum. The latter feature is perfectly natural, given the nature of these terms as probability leaking and (mis-) normalization at zero distance. 7 The detailed derivation of eq. (2.6) is carried out in section 3. While in section 4, we derive an exact analytic expression for the matter dependent part of the oscillation probability (2.6). The combinations 7 A comment is ready for the normalization term, the second term in (2.6). Its original form is 3 j=1 (U X)αj(U X) * βj , which is natural because it comes from the contribution of zeroth-order Hamiltonian with all orders effect of the matter potential. It is easily reduced to the vacuum form in (2.6) (or in (2.2)) by using unitarity relation 3 j=1 X kj X * lj = δ kl .

JHEP02(2019)015
of X matrix elements that used in the derivation can also be utilized to calculate higher order corrections in W . In section 5.1, the regions of visible effect of unitarity violation is illuminated by plotting the probabilities with/without unitarity violation in wide ranges of E and baseline L.
After understanding the general feature of perturbative series based on explicit calculation to order W 4 , we postulate the "Uniqueness theorem" which states that the oscillation probability formula eq. (2.6) is valid to all orders in W expansion under the same conditions on the sterile state mass and the kinematical region as used in the discussion of fourth-order formulas. See section 6.3. The reasons for this interesting feature, the same mass conditions as in vacuum to guarantee the sterile sector model independence prevail in matter, will be partially explained at the end of section 3.5.
Finally, but probably most importantly, we point out that outside the region of validity of our above theorem, there are regions of neutrino energy and baseline that condition for suppression due to the large sterile state mass denominators is not fully effective. We show that in such region, second order correction terms in W , together with the leaking term C αβ , may not be totally negligible, and it could be detectable. It would offer yet another way of distinguishing low-scale unitarity violation from high-scale one. These new terms are derived in section 3.4 and their effects are quantified in section 5.2.

Small unitarity-violation perturbation theory of neutrino oscillation in matter
We formulate a perturbation theory of the (3 + N ) state unitary model using an expansion parameter of matrix elements of W signifying unitarity violation effect, assuming it small. It will be done aiming at constructing a model-independent framework for leptonic unitarity test. It necessitates the conditions on the sterile neutrino mass as discussed in section 2.1. In most of the discussions in this section we presume, as an appropriate setting for unitarity test, terrestrial neutrino experiments, i.e., accelerator LBL experiments, and/or atmospheric neutrino measurement. Use of reactor and accelerator neutrinos at short baselines offers an alternative way for testing leptonic unitarity but with only minor matter effect.
In the main text we mostly confine ourselves to the formulas to second order in W , but include fourth order terms whenever it is necessary. We take for simplicity the uniform number density approximation for electrons and neutrons in matter. However, extension to the varying density case is, in principle, straightforward as far as adiabaticity holds. Usage of the same probability formula as a hunting tool of unitarity violation and discriminator between low-scale and high-scale unitarity violation will be discussed in section 5.

JHEP02(2019)015
and the oscillation probability is given by The neutrino evolution in flavour basis in the (3 + N ) space unitary model is governed by the Schrödinger equation Given the flavour basis Hamiltonian H, the S matrix is given by where T symbol indicates the "time ordering" (in fact "space ordering" here). The righthand side of (3.4) may be written as e −iHx for the case of constant matter density. The flavour basis Hamiltonian H is (3 + N ) × (3 + N ) matrix: Here, m i (m J ) denote the mass of mostly active (sterile) neutrinos and E is the neutrino energy. ∆ A and ∆ B are related to Wolfenstein's matter potential [4] due to charged current (CC) and neutral current (NC) reactions, a and b, as eigenstate basis as ν ζ = U ζzνz , where ζ runs over active flavour α = e, µ, τ and sterile flavour s = s 1 , · · ·, s N indices, z runs over mostly active i = 1, 2, 3 and mostly sterile mass eigenstate J = 4, 5, · · ·, N + 3 indices. For simplicity, we introduce a compact notation which writes (3 + N ) × (3 + N ) matrix in a form of 2 × 2 matrix. By defining the active 3 × 3 matter potential matrix the flavour basis Hamiltonian is written as where ∆ a = diag(∆ 1 , ∆ 2 , ∆ 3 ) and ∆ s = diag(∆ 4 , ∆ 5 , · · ·, ∆ N +3 ).
As an application of our framework, we anticipate leptonic unitarity test in the LBL accelerator neutrino experiments which utilize atmospheric-scale neutrino oscillations. Therefore, we assume that the system satisfies the following conditions in formulating our perturbation theory where L denotes baseline, ∆m 2 ji ≡ m 2 j − m 2 i , and They probably ensure that our oscillation probability formulas have applicability to the terrestrial LBL and atmospheric neutrino experiments with baseline up to ∼ 10 4 km and energies from low to high, up to E ∼ 100 GeV. More precise discussions on where our formulas are valid will be given in sections 3.5 and 6.2.

Vacuum mass eigenstate basis, or tilde basis
To formulate perturbative treatment it is convenient to consider the vacuum mass eigenstate basis, the tilde basis, introduced in the previous sectioñ The tilde basis Hamiltonian is related to the flavour basis one as (3.14) The explicit form ofH is given bỹ

JHEP02(2019)015
We parameterize the (3 + N ) × (3 + N ) dimensional flavour mixing matrix U as The matrix U and V are 3 × 3 and N × N matrices, respectively, and W and Z have sizes that just fill in the space. In our (3 + N ) model, unitarity is obeyed in the whole (3 + N ) state space: Then, the HamiltonianH in vacuum mass eigenstate basis is given bỹ As in vacuum, the neutrino oscillation is governed only by the U and W matrices, and is independent of Z and V matrices. It is natural that V matrix does not show up in physical Hamiltonian matrix because the rotations inside sterile basis does not have any physical meaning, if we observe the system only by the Standard Model interactions. However, the flavour basis Hamiltonian H in (3.5) obviously depends on Z and V . The apparent puzzle will be resolved in appendix A.

Formulating small unitarity-violation perturbation theory
We now construct the small unitarity-violation perturbation theory. It is natural to consider the framework in which the tilde-basis HamiltonianH is decomposed into the un-perturbed and perturbed parts,H 0 +H 1 , as follows: Therefore, what we mean by "expansion by unitarity violation effect" is an expansion by the W matrix elements. 8 We assume, for simplicity, that all the W matrix elements are small and have the same order s . Then, 3 × N (N × 3) sub-matrix elements inH matt are of order s , while the pure sterile space N × N sub-matrix elements are of order 2 s . For simplicity, we often use the expression "expanding to order W n " which means to order n s in this paper. 8 Through unitarity (3.17), U matrix elements have some dependence on W matrix elements. We choose not to expand U matrix elements by this W dependence. In this sense, we use a "renormalized basis" (in the same sense as in ref. [53]) in which some higher order effects are absorbed into the zeroth-order state.

Hat basis
To formulate perturbation theory withH 0 andH 1 given above we transform to a basis in which the un-perturbed part of the Hamiltonian is diagonal, which we call the "hat basis". Since the 3 × 3 sub-matrix ∆ a + U † AU inH 0 is Hermitian, it can be diagonalized by the unitary transformation with X being the 3 × 3 unitary matrix. Then,H 0 can be diagonalized by using the zeroth-order Hamiltonian in the hat basis. SinceĤ 0 is diagonal it is easy to compute e ±iĤ 0 x : (3.23) Then, the perturbed Hamiltonian is given bŷ The eigenvalues ofH 0 is therefore h 1 , h 2 , h 3 , and ∆ J (J = 4, · · ·, 3 + N ). Therefore, the sterile neutrino masses are affected neither by the active states nor the matter potential in our zeroth-order unperturbed basis. It must be a good approximation because we have assumed that the sterile neutrino masses are much heavier than the active ones, and we are interested in the energy region implied by a ∼ ∆m 2 31 . To do real calculations of the S matrix elements we must solve the zeroth order Hamil-tonianH 0 . This task will be carried out in section 4.2, in which we derive explicit expressions of the eigenvalues h i and the unitary matrix X. Now, we formulate perturbation theory with the hat basis Hamiltonian,Ĥ 0 in (3.22) andĤ 1 in (3.24) after a clarifying note in the next subsection.

The relationship between quantities in various bases
So far we have introduced the tilde-and the hat-basis:

JHEP02(2019)015
where X is given by eq. (3.21). Therefore, Notice that both U and X are unitary, and hence UX is unitary too. The relationship between wave functions of various basis are given bŷ where y denote the hat-basis indices. Using the explicit parametrization of the U matrix we have It may be helpful for our discussions later to understand the relationship between S andŜ matrix elements. For this purpose, we denote them in the block form where the subscripts a and S indicate that they act (for the right index) to the active or the sterile subspaces. Notice that S aa and S aS , for example, are 3 × 3 and 3 × N matrices, respectively. Then, the relationship between S andŜ matrix elements can be written explicitly as

Computation ofŜ matrix elements
To as Ω(x) obeys the evolution equation

JHEP02(2019)015
where Then, Ω(x) can be computed perturbatively as where the "space-ordered" form in (3.35) is essential because of the highly nontrivial spatial dependence in H 1 . Upon obtaining Ω(x),Ŝ matrix can be obtained aŝ By knowingŜ matrix elements, the S matrix is obtained by using (3.27), or (3.31).
The perturbing Hamiltonian H 1 defined in (3.34) has a structure (3.37) That is, (H 1 ) ij = 0 in the whole active neutrino subspace. The non-vanishing elements of H 1 are as follows: Inserting eq. (3.38) into (3.35), we can compute all the Ω matrix elements. The simplest ones in first order in H 1 , the second term in (3.35), are given by

JHEP02(2019)015
which serve as a building block of the perturbation series because of the structure in (3.35). The notation " [1]" implies that the terms come from first order perturbation with H 1 . For more about notations, see appendix B. We need to compute up to fourth order in H 1 because we want to keep all the order W 4 terms. The requirement arises because the probability leaking term, whose observation is crucial to distinguish between low-energy and high-energy unitarity violation, is of order W 4 . The other normalization term, the second term in (2.2), also deviates from the one in unitary case by a quantity of order W 4 in the appearance channels, but in an implicit way. The resulting expressions ofŜ matrix elements to order W 4 are summarized in appendix B.
There exists important consistency check in the calculation. That is, the identity relation betweenŜ matrix elements that follows from generalized T invariance: 9 whereŜ Ji is obtained by performing the exchange h i ↔ ∆ J inŜ iJ . The generalized T invariance relation is explicitly verified by the computed results ofŜ matrix elements to fourth order in W given in appendix B. 10

Computation of S matrix elements
Given the results ofŜ matrix elements it is straightforward to calculate S matrix elements by using the formulas in eq. (3.31). The active neutrino space S matrix elements can be written in perturbative forms, Using (3.41) the explicit expressions of S matrix elements can be easily obtained with use ofŜ matrix elements given in appendix B. For example, S αβ in zeroth and second orders in W are given, respectively, by As in the Standard Model in particle physics T invariance is broken in our system only by complex numbers in the mixing matrix. 10 SinceĤ system is a consistent dynamical system it is legitimate and easier to verify generalized T invariance in theŜ level, though it can be done in the S matrix level as well. A pedagogical treatment for proving generalized T invariance is given in version 1 of this work, arXiv ePrint: 1712.02798.

The oscillation probability to second order in W
In this section, we discuss the oscillation probability to second order in W . It is to illuminate the principle of calculation, how averaging over the fast oscillation works, and to show which constraints are obtained on the sterile state masses by the requirement of suppression by the large sterile state mass denominators to make these sterile-sector model dependent terms negligible. Of course, we will calculate in this paper all the oscillation probabilities P (ν β → ν α ) in matter to fourth order in W to keep the necessary term, the probability leaking term C αβ , as mentioned earlier. The key features of the fourth-order terms will be described in the next section 3.6.
The oscillation probability P (ν β → ν α ) is given to second order in W as
Notice that there is no matter dependent terms without suppression either by highfrequency oscillations ∝ cos(∆ K −h m )x (or sin), or by large sterile state mass denominators ∝ We take averaging over fast oscillations due to active-sterile and sterile-sterile mass squared differences which leads to where . . . stands for averaging over neutrino energy within the uncertainty of energy resolution, as well as averaging over uncertainty of distance between production and detection points of neutrinos. 11 The second approximate equalities in (3.45) assume that there is no accidental degeneracy among the sterile state masses. That is, we assume that the relation |∆m 2 JK | |∆m 2 31 | always holds. After averaging out the fast oscillations, P (ν β → ν α ) is given to second order in W as The zeroth-order term P (ν β → ν α ) (0) is nothing but the one in eq. (2.6) except for dropping the probability leaking term while the W 2 correction terms are given by

JHEP02(2019)015
The W 2 correction terms in P (ν β → ν α ) (2) , together with the probability leaking term C αβ in eq. (2.6), will be utilized in section 5.2 to explore the possibility of distinguishing between low-and high-scale unitarity violation. If such terms are detected, the sterile sector model-dependence in P (ν β → ν α ) (2) would serve for identifying the structure of the sterile sector.

Suppression by the large sterile state mass denominator
In this section, we study the conditions under which P (ν β → ν α ) (2) in eq. (3.48) can become negligibly small. It would allow us to use P (ν β → ν α ) (0) + C αβ (= eq. (2.6)) for leptonic unitarity test in a sterile sector model-independent manner. We start by examining the effect of suppression by the large sterile state mass denominator which characterizes transition between active-sterile states, 1/(∆ K − h k ). We demand that the matter dependent terms in (3.48) be smaller than the probability leaking and the normalization terms of order ∼ W 4 . It leads to where L is the baseline distance and i and J denote, respectively, generic indices for active and sterile states. For notational convenience, we define λ i (i = 1, 2, 3) to be the eigenvalues of 3 × 3 submatrix 2EH 0 in (3.19) corresponding to the active neutrino mass squared in matter and hence λ i = 2Eh i . In region λ i ∼ |∆m 2 31 | and near the atmospheric oscillation maximum, L ∼ 2E GeV , which further suppresses the first and the second items in (3.49) unless ρE 10 (g/cm 3 )GeV. Therefore, in this region the last one in (3.49) gives the severest constraint (taking the matter potential due to CC in A and removing the factor 1 2E ) (3.50) Notice that, in order for the first inequality in (3.50) to be valid, we have restricted the energy region for a given matter density such that λ i remain in the order of active neutrino masses. Roughly speaking, it corresponds to −50 (g/cm 3 )GeV Y e ρE 50 (g/cm 3 )GeV where the negative sign is relevant for antineutrinos. See e.g., figure 3 of ref. [53]. Clearly, it excludes the interesting region of "IceCube resonance" due to sterile neutrino mass of eV scales [54], for which an entirely different theoretical framework would be necessary. Then, we notice that in a regime |W | 2 ∼ 10 −2 , the condition in (3.50) is valid given the estimation (assuming Y e = 0.5) unless ρE 10 (g/cm 3 )GeV. That is, the second-order matter dependent correction terms can be ignored in comparison with O(W 4 ) terms if ∆m 2 Jk 0.1 eV 2 , which is already JHEP02(2019)015 required in vacuum. If we want to treat the regime |W | 2 10 −n , we need to limit the sterile masses to ∆m 2 Jk m 2 J 10 (n−3) eV 2 to keep our (3 + N ) space unitary model insensitive to details of the sterile sector [20]. We note, however, that terms of order |W | 4 ∼ 10 −4 may be the limit of exploration for near future neutrino oscillation experiments.
The condition (3.50) is identical with the one obtained using the first order matter perturbation theory [20], which may look strange to the readers. Let us understand the reason why taking care of all order matter effect does not alter the condition obtained by first-order treatment in matter perturbation theory. The matter-dependent term in the zeroth-order HamiltonianH 0 only involves U matrix, but no W matrix. Since we treat H 0 in an unperturbed fashion it produces all-order effect of the matter potential which is however independent of W matrix elements. On the other hand, perturbative effects that come from single or double powers in W inĤ 1 are always accompanied by the matter potential in the form of W A or W † A, as in eq. (3.37). That is, perturbative effect of W is always accompanied by matter potential, and hence can always be dealt with matter perturbation theory. 12 It is the reason why the matter perturbation theory is able to yield the same condition on sterile masses as obtained in a fuller treatment of matter effect done in this paper.

The oscillation probability in fourth order in W
The oscillation probability in fourth order in W contains the two terms (3.52) We will show in appendix C.1 that the first term in (3.52), after averaging over the fast oscillations and using the suppression by large sterile state mass denominator as discussed in the previous section, leaves the unique term, the probability leaking term C αβ in eq. (2.6), which can be seen in eq. (C. 6). An interesting feature of C αβ in matter is that it is identical to the one in vacuum, eq. (2.3) without any matter effect dressing. In our computation the term comes from the hat basis S matrix in zeroth order, the first term in the last line of eq. (B.4), and hence it is free from the matter potential. 13 We will also show in appendix C.2 that the second term in (3.52), under the same treatment for the first term, gives vanishing contribution. Therefore, no matter-dependent fourth order term survives after large sterile state mass denominator suppression is used and averaging over the fast oscillations is performed.
In conclusion, the oscillation probability in matter between active flavour neutrinos in the (3 + N ) space unitary model to fourth order in W in our small unitarity-violation 12 An example of this feature can be observed in eq. (7.13) in ref. [20]. We must remark, however, that this reasoning does not prove that the first order in matter perturbation theory is sufficient to obtain all the necessary conditions on the sterile state masses. 13 One may suspect that including higher order corrections could alter the feature of matter potential independence of C αβ . However, one can show (see section 6.3) that higher order W corrections to the piece of S matrix elements relevant to C αβ organize themselves as a phase factor, so that C αβ has no matter effect dressing. The rest of the correction terms are suppressed due to the dimensional reason, an extra matter potential must be accompanied by an energy denominator.

JHEP02(2019)015
perturbation theory can be written as in eq. (2.6) in section 2. We hope that it serves as a useful tool to test leptonic unitarity in various ongoing and future neutrino oscillation experiments.

Analytical and numerical methods for solving non-unitary evolution in matter
In this section, we describe the numerical and analytical methods for calculating the neutrino oscillation probability by solving non-unitary evolution in matter.

Numerical method for calculating neutrino oscillation probability
We describe a numerical method for computing the oscillation probability in matter. This method can be used, assuming adiabaticity, in cases with varying matter density. We show that in zeroth order in W the system simplifies to an evolution equation in the 3 × 3 active subspace.
We solve the Schrödinger equation in the vacuum mass eigenstate basis ("tilde basis"), ν z = (U † ) zζ ν ζ with HamiltonianH in (3.18): where i = 1, 2, 3 and J = 4, 5, · · ·, 3 + N denote mostly active and mostly sterile neutrino mass eigenstate labels, respectively. The initial condition with only active component impliesν Therefore, in the mass-basis formulation only U and W are involved, which is consistent with our experience in W perturbation theory. An apparent contradiction to this property that one faces in the evolution equation in the flavour basis is resolved in appendix A. A drawback of this method is that we have to solve explicitly the evolution of the sterile states which are coupled to the active states. Then, we need to specify the sterile sector model, and have to know how to deal with averaging over the fast modes.
We notice, however, that in the zeroth-order in W the system simplifies. Since the HamiltonianH is block-diagonal it suffices to solve the equation only in the 3 × 3 active neutrino subspace:

JHEP02(2019)015
The initial condition (4.2) and final reverse-back formula (4.3) involve only U matrix elements. Therefore, the oscillation probability in the zeroth-order in W can be calculable in a manner independent of sterile sector models. 14

An exact solution of zeroth-order oscillation probability
Here, we describe a method for obtaining the analytical solution of the zeroth-order Hamiltonian. The exact solution, as well as the numerical one described in the previous section, provides the basis for computing the higher order corrections in W . We calculate an exact form of the oscillation probability P (ν β → ν α ) in leading order in our perturbative framework, the one in (2.6) except for C αβ , in the case of uniform matter density.
The zeroth-order S matrix element S (0) αβ in (3.42) can be written as 5) and the factor in parenthesis can be calculated by the KTY technique [55]. We want to diagonalize the Hamiltonian For our notational convenience we call this form of H 0 as H d . Note that h i = λ i 2E where where with (4.10) 14 As we remarked in footnote 8 the non-unitary mixing matrix U has some W dependence through unitarity of the U matrix in the whole (3 + N ) space. Therefore, the nature of the eq. (4.4) as the zeroth-order in W is ambiguous. However, following [20], we remain in the treatment with this "W effect renormalized basis" in this paper.

JHEP02(2019)015
The adjugate of H 0 is defined as AdjH 0 ≡ (H 0 ) −1 detH 0 . Notice that T , A and D are invariant under unitary transformation of H 0 → KH 0 K † with K any unitary matrix and so are λ i . Following the notation in [55] we define p ij and q ij as (i, j = 1, 2, 3) (4.11) Notice that p ij and q ij are written only by the known (or given) quantities. Then, the equations together with unitarity of X, become the equations to determine XX † : They lead to the solution (k = 1, 2, 3) where k, l, m is cyclic, and sum over k is not implied in (4.14). Therefore, to zeroth-order in W expansion, the S matrix elements are given by , (4.15) and the oscillation probability by P (ν β → ν α ) = |S (0) αβ | 2 . Finally, armed with the solution (4.14), we can also calculate all higher order terms in oscillation probability for e.g. those in eq. (3.48) since only such combination X ik X * jk (no sum over k implied) can appear.
5 Where are the unitarity violation and W 2 corrections?
Having formulated the small unitarity violation perturbation theory, we now utilize it to answer the following questions: (1) Where is the regions of energy E and baseline L in which the effect of unitarity violation is significant?, and (2) how large can the W 2 corrections be? We address the questions (1) and (2) in sections 5.1 and 5.2, respectively.

Comparison between the oscillation probabilities with and without unitarity violation
To know where the effect of unitarity violation is large, and how large it is, we calculate

JHEP02(2019)015
as a function of E and L, where P (ν β → ν α ) standard and P (ν β → ν α ) non-unitary imply, respectively, the oscillation probabilities calculated with the standard unitary mixing matrix and the leading order (i.e., W 0 ) one with non-unitarity. The probability leaking term C αβ in eq. (2.6) as well as the W 2 correction terms in eq. (3.48) are not included in the analysis here. Therefore, the results given in section 5.1 apply to both high-scale unitarity violation as well as low-scale one in its leading order in W . 15 On intuitive ground, at zeroth order in W our system describes high-scale unitarity violation. There is no "W corrections" in high-scale unitarity violation because the energy scale is so high that the high-mass sector is truncated. We examine the three channels ν µ → ν e , ν µ → ν τ , and ν µ → ν µ . However, we do not enter into any quantitative analyses, nor attempt to cover the whole parameter space.
Here is a brief note on how the standard mixing and the unitarity violating parameters are chosen: we take the (3 + 1) model in which the constraints on the parameters are best understood [27,[56][57][58]. In consistent with the current constraints we have chosen: sin 2 θ 14 = 0.02, sin 2 θ 24 = 0.01, and sin 2 θ 34 = 0.1 for ∆m 2 41 = 0.1 eV 2 , and set all the CP phases to zero. Then, we cut out the 3 × 3 active neutrino mixing matrix, which is nonunitary. 16 For the standard leptonic mixing parameters in U PDG , we take sin 2 θ 12 = 0.3, sin 2 θ 23 = 0.5, sin 2 (2θ 13 ) = 0.09, and the mass squared differences ∆m 2 21 = 7.4 × 10 −5 eV 2 and ∆m 2 31 = 2.4×10 −3 eV 2 , and set the CP phase δ CP to zero. The uniform matter density is taken as ρ = 3.2 g cm −3 over the entire baseline, which may not be realistic.
non-unitary is small. However, we identify the two regions where P (ν µ → ν e ) (0) non-unitary is relatively large, 0.3. One of them is at low energy, E a few hundred MeV, and baseline L 1000 km. The other one is a region E ∼ 10 GeV and L ∼ 10000 km. The former may be understood as due to the solar MSW enhancement, and the latter as the atmospheric MSW enhancement [3,4]. Roughly speaking, the regions with relatively large |∆P (ν µ → ν e )| overlap with these regions. 15 In fact, it is in agreement with the formulations in ref. [28] with which we share the same evolution equation (4.4) in the vacuum mass eigenstate basis. See also [21]. However, it appears that the flavour basis formulation of neutrino evolution in matter in high-scale unitarity violation poses some nontrivial features such as non-Hermitian Hamiltonian [21], or the evolution equation i d dx να = j U ∆a + U † AU U † αβ ν β [29]. The latter is not equivalent to (4.4) in the vacuum mass eigenstate basis due to non-unitarity of the U matrix. 16 It can be re-parameterized in terms of the "α matrix parameterization" defined in ref. [26].
non-unitary , roughly speaking, the relation P (ν µ → ν µ ) non-unitary is small. It must be the case in the unitary case, but even in non-unitary case the relation holds approximately because unitarity violation is small in our choice of the parameters. Therefore, P (ν µ → ν µ ) non-unitary is small, and vice versa, as seen in figure 3. It appears that the anticorrelation is inherited to the relationship between ∆P (ν µ → ν µ ) and ∆P (ν µ → ν τ ). Relatively large ∆P (ν µ → ν τ ) in first a few oscillation maxima, or similar large depletion of ∆P (ν µ → ν µ ), would allow detection of non-unitarity if the detector has a good τ (in the former channel), or µ (in the latter channel) detection capabilities. If the detector can detect the both, anticorrelation between µ and τ yields must help.
Some comments on observational aspects: in the two regions where |∆P (ν µ → ν e )| is large, and ∆P (ν µ → ν µ ) in energy region E 10 GeV may be explored by high-statistics atmospheric neutrino observation by Super-K, Hyper-K/HKK, or DUNE [46,[48][49][50]. The atmospheric MSW enhanced region of P (ν µ → ν e ) would be a good target for PINGU extensions of IceCube and KM3NeT-ORCA [51,52]. P (ν µ → ν τ ) and P (ν µ → ν µ ) would be explored by them, with possibility of seeing anticorrelation between µ and τ yields. Although it is very interesting to investigate these experimental prospects, a detailed examination of these questions is beyond the scope of this paper.

5.2
The probability leaking and W 2 correction terms

Low-scale versus high-scale unitarity violation
In leptonic unitarity test, a clear understanding of the relationship between low-scale and high-scale unitarity violation may be one of the key issues. We have stressed in our previous paper [20] that observing the probability leaking term C αβ in eq. (2.3) would testify for lowscale unitarity violation. As mentioned in section 3.6 the leaking term is not dressed by the matter effect, which is perfectly natural for the effect of probability leakage. In this paper, we propose yet another way of distinguishing low-scale unitarity violation from high-scale one. That is, detection of the W 2 correction terms in eq. (3.48). In this section 5.2, we give a brief sketch of how and where we might see visible effects of the probability leaking and the W 2 correction terms.

JHEP02(2019)015
non-unitary is presented. The parameters used are the same as in figure 1.  Let us go back to the expression of the oscillation probability to second order in W , eq. (3.48), in section 3.4 to know where we might see visible effects. If we enter into the region ρE 10 (g/cm 3 )GeV at around the first oscillation maximum, the first two terms in eq. (3.48) can become large apart from W 2 suppression,

JHEP02(2019)015
Comparing to the conditions (3.49), they can be larger than W 4 terms and hence cannot be neglected. After taking account of W 2 suppression of ∼ 0.01 (assuming W 0.1), | AAL (∆ J −h i ) W 2 | ∼ 3 × 10 −2 at E ∼ 100 GeV, assuming ∆m 2 Ji = 0.1 eV 2 . To know more quantitatively the sizes of W 2 corrections and their E or L dependences, we have to fix the W matrix elements which have large arbitrariness. We defer this technical discussion to appendix D, which describes the recipe we took to fix them with a common m 2 J = 0.1 eV 2 . 18 We plot in figure 4, δP (ν µ → ν α ) ≡ P (ν µ → ν α ) (2) + C µα , that is, the order W 2 correction terms in P (ν µ → ν α ), eq. (3.48), plus the probability leaking term C µα , α = e (top panel), α = τ (middle panel), and α = µ (bottom panel). In other words, δP (ν µ → ν α ) is equal to the total probability minus P (ν µ → ν α ) (0) , if the fourth and the higher-order in W correction terms with matter are neglected. In each panel the three cases are examined. N = 1 case with maximal C µα (solid line), the universal scaling model 19 with N = 3 (dotted line), and the order W 2 correction only (dashed line). The last case corresponds to the universal scaling model with N = ∞. The blue lines are for E = 10 GeV, and the red for E = 100 GeV.
At longer distance and in appearance channels, we see enhancement. At E = 10 GeV, we observe a factor of several enhancement in |δP (ν µ → ν α )| for both α = e and α = τ 18 We are aware that the assumption of equal sterile neutrino masses is contradictory to the assumption of no accidental degeneracy in the sterile mass spectrum we made in section 3.4. It was done not to complicate term by term evaluation of the perturbative series, and to avoid using degenerate perturbation theory. Fortunately, we can remove this assumption to second order in W in which no purely sterile state mass splitting denominator is involved. 19 The universal scaling model is defined in appendix E. It prescribes a way of distributing Wα4 matrix element in 3 + 1 model to the W matrix elements in 3 + N model in such a way that the size of order W 2 correction terms in (3.48) remains unchanged when all the sterile masses are equal. 20 Of course, there is an issues of how to separate effects of W 2 correction terms from unitarity violation through U matrix in leading order.

3) for definition) in
is plotted assuming a common m 2 J = 0.1 eV 2 . The top, middle and bottom panels are for α = e, τ , and µ, respectively. In each panel the three cases are shown: N = 1 case with maximal C µα (solid line), the universal scaling model with N = 3 (dotted line), and the order W 2 correction terms only (dashed line). The last case corresponds to the universal scaling model with N = ∞. The blue lines are for E = 10 GeV, and the red for E = 100 GeV. The leaking terms in the N = 1 model (shown without superscript (N = 1) in the legend) have values C eµ = 2 × 10 −4 , C τ µ = 9.5 × 10 −4 , and C µµ = 9.6 × 10 −5 .
A final remark on C αβ vs. W 2 corrections. Since C αβ is a constant term in the oscillation probability, it can in principle be distinguished from the other normalization term which shares U matrix element dependences with the oscillation terms. In particular, they can dominate for large m 2 J since the W 2 correction terms are suppressed by at least ∼ 1/m 2 J . In this case, they will be the sole indicator of low-scale unitarity violation. In general (though not in the N = 1 model), the order W 2 terms depend upon details of the sterile sector, e.g., matrix structure of W . Therefore, once the effect is seen it would give us useful information on the structure of low-scale leptonic unitarity violation.

Some remaining theoretical issues and extending
In this section, we will give some remarks on the theoretical basis in our framework, basic one as well as on its perturbative aspects. They include our treatment of decoherence, generic structure of higher-order corrections and its relation to the "Uniqueness theorem" (see section 2.3), absence of enhancement due to small solar mass splitting denominator, and its relation to the other non-standard physics.

Decoherence imposed onto coherent evolution system
We have started with the Schrödinger equation (3.3) with Hamiltonian (3.5) assuming that all the neutrino states remain coherent. We have shown in this and the previous papers that the coherence between active and sterile, and sterile and sterile states are not maintained for sterile mass differences larger than 0.1 eV 2 . The effect of decoherence is taken into account by making average over the fast oscillations. We feel it desirable for the current treatment be replaced by the real quantum mechanical one using wave packets, in which the effect of decoherence would automatically come in. Yet, we do believe that our present framework is able to describe effectively the right physics derived from such improved treatment.

Smallness of expansion parameters and higher order corrections
Here, we discuss general structure of the perturbation series without recourse to averaging out the fast oscillations. The effective expansion parameters in our perturbative framework are the following four,

JHEP02(2019)015
We already saw them, except for the last one, in the discussion in section 3.5, and it can be seen by inspecting the expressions of the oscillation probabilities up to the fourth orders given in section 3.4 and appendix C. Formally, the expansion parameter is the first one in (6.1) in view of (3.35) with Ω [1], the kernel, in (3.39). But, the spatial integration in (3.35) produces different effective expansion parameters, the second and the third ones in (6.1). The extra factor of W 's without the kinematical factors is provided when transforming from theŜ to S matrices, as seen in section 3.3.4. For simplicity of the discussion in this section, we limit ourselves to the case of |W | 0.1. Under the same conditions we have imposed in section 3.5, the first one in (6.1) is 7.6 × 10 −4 for ∆m 2 Ji = 0.1 eV 2 and ρE = 10 (g/cm 3 )GeV while the second and the third, which are comparable to each other at around the first oscillation maximum, are estimated to be 2.3 × 10 −2 . Therefore, the smallness of the expansion parameter is ensured unless ρE 10 (g/cm 3 )GeV. In fact, a close examination of the order W 4 terms in the oscillation probability (see appendix C) shows that all the formally W 4 terms are actually further suppressed. The largest term in the fourth-order oscillation probabilities is of the one suppressed by a factor AW ∆ J −h i (ALW ) W 2 1.7 × 10 −7 , which is as small as ∼ 10 −4 even in the case |W | = 0.5. Therefore, we expect that the formula for the oscillation probability in (2.6) works under much relaxed conditions than the one in (3.50).

On Uniqueness theorem and matter-dependent dynamical phase
We have shown in sections 3.4 and 3.6 that there is no surviving matter dependent correction term in the oscillation probability up to order W 4 after averaging out fast oscillations and using the suppression by large sterile state mass denominators. Should we expect that this feature is stable against higher order corrections beyond order W 4 ? We argue that the answer is Yes. Based on the feature of perturbative series we have learned, we postulate the following theorem: Uniqueness theorem. All the matter dependent perturbative corrections in W in the oscillation probability either vanish or can be ignored after averaging over the fast oscillations and using the suppression due to the large sterile state mass denominators, leaving only the probability leakage term C αβ , the first term in eq. (2.6) with (2.3).
It must be remarked here that unitarity violation effects which are hidden in non-unitary active space mixing matrix U produces zeroth-to higher order effects of W . The above theorem is only about the terms generated by explicit perturbative corrections in W . We first note that higher-order corrections in terms of W are computed by using Ω [1] as the kernel, as indicated in eq. (3.35). Notice also that all the elements of Ω [1], except for Ω [1] JJ , carry the sterile state mass denominator, as shown in (3.39). Then, higher order correction terms are always accompanied by the sterile state mass denominators which are composed of some of the first three in (6.1), and therefore they are suppressed. The unique exception for it is the terms generated only by Ω [1] JJ which lacks the sterile state mass denominator. Therefore, apart from this special case, we have shown that higher-order corrections in W does not produce the surviving terms after averaging over fast oscillation and using the sterile state mass denominator suppression. It is consistent with what JHEP02(2019)015 we saw in our explicit computation to order W 4 . This concludes our justification of the Uniqueness theorem.
We need to clear up the issue of special type of perturbative correction terms which involve only Ω[1] JJ as the kernel in (3.35). It produces the unique form ofŜ JJ aŝ a collection of terms of matter-dependent higher order renormalization to J W αJ W * βJ , the probability leaking term at the amplitude level. However, it exponentiates and has contribution to the S matrix element as 21 3) The unique form of S matrix, in principle, raises an interesting issue of dynamically generated phase produced jointly by unitarity violation and the matter effect. 22 In our setting, however, it either disappears from the amplitude squared, or has vanishing effect when the high frequency oscillation is averaged out. Finally, we should remark that our discussion to justify Uniqueness theorem in this section assumes the same kinematical region as in the treatment of order W 2 and W 4 correction terms in sections 3.4 and 3.6, in particular, ρE 10 (g/cm 3 )GeV. However, as mentioned at the end of section 6.2, it is likely that the region of validity of the probability formula (2.6) with vanishingly small higher order corrections is wider. At present, the precise boundary of kinematical region for its validity is not known to us.

Absence of enhancement due to small solar mass splitting denominator
In perturbation theory one has to sum up intermediate states including off mass shell states. Therefore, even though we sit in the kinematic region where atmospheric-scale oscillations are large, the denominator can become small, to the order of solar ∆m 2 mass splitting. Then, one might question whether the correction terms blow up at the small denominator, which would invalidate our perturbative treatment.
Fortunately, one can show that the "singularity" which could be produced in the limit of small solar mass splitting always cancels against the small numerator of the similar size. This problem exists already in the second-order expression of the oscillation probability (3.44). See the second term in second order (in W ) term. If we denote h l − h k ≡ the term would have 1/ singularity in the limit of → 0. However, one can see by inspection by eye that the expression inside the square parenthesis is antisymmetric under l ↔ k, and hence it is of order or higher. Therefore, the singularity cancels. Notice that the antisymmetry under l ↔ k is not required for the whole expression including the matrix element factor. 21 It might be easier to obtain the phase factor if we use a different decomposition ofH from (3.19) by absorbing W † AW intoH0. 22 The phase itself needs not be small. Taking the matter potential of CC reaction and the earth diameter, . Therefore, ALW 2 can be order unity for |W | 0.4.

JHEP02(2019)015
The situation is a little bit more complicated in the fourth-order expression of the oscillation probability given in appendix C. In addition to 1/ singularity similar to the one we already saw, there exist apparent singularity of 1/ 2 type. See, for example, the second term in (C.11) and the last term in (C.12). But, an explicit calculation shows that the 1/ 2 singularity always cancels against order 2 numerator in the limit of small solar splitting.
This phenomenon is reminiscent of the finiteness of the oscillation probability at the small solar mass splitting limit in helio-perturbation theory with the unique expansion parameter (or a renormalized one), see e.g., [53] and the references therein. Possible interpretation of applicability of the perturbative framework to the region of solar level crossing has been discussed [59,60]. Another example for the similar phenomena is the one at the small atmospheric mass splitting limit with additional expansion parameter sin θ 13 . In this case it is observed that near the atmospheric resonance region not only the oscillation probability is finite but also its accuracy improves when the higher order terms to fourth order in sin θ 13 is added [61].
Then, one might ask if our small unitarity violation perturbation theory gives quantitatively accurate result at around the denominator with small solar mass splitting. However, we note that this problem is not relevant in our case because all these terms with apparent singularities vanish after averaging over the high-frequency oscillations and using the suppression by the large sterile mass denominators. Yet, we must remark that if we investigate possible enhancement of the correction terms outside the condition (3.50), as done in section 5.2, the quantitative accuracy of the expression may become an issue.

Leptonic non-unitarity and the other non-standard physics
This final subsection is to mention the related but different approaches, and to make some clarifying remarks. Our 3 + N model has obvious relation with the various versions of active plus sterile neutrino models proposed in the context of LSND-MiniBooNE anomaly, as reviewed in [62], see also the references therein. The clear difference exists in the attitude of the treatment of the model, in our case seeking the conditions to make the predictions as model-independent as possible, while in the others pursuing the particular model which provide the best fit to the data. Unless we use our model-independent simplified formula (2.6), we would have to marginalise over the huge parameter space of the (3 + N ) model to obtain the bound on non-unitarity 3 × 3 mixing matrix U .
It is pointed out that one can establish a mapping between parameters in the mass eigenstate basis which describe non-unitary leptonic mixing and the ones for non-standard neutrino interactions (NSI) under certain conditions between neutron and electron number densities [28]. However, when we rotate back to the flavour basis, a non-unitary mixing matrix is involved in unitarity violating case, but not in the NSI case, as far as propagation in matter is concerned.

Concluding remarks
In this paper, we have presented a comprehensive treatment of the three active plus N sterile neutrino model in the context of leptonic unitarity test. We have formulated an JHEP02(2019)015 appropriate perturbative framework with expansion in small unitarity violating W matrix elements, while keeping (non-W suppressed) matter effect to all orders.
What we have done in this paper is mainly threefold: • We have shown that the oscillation probability in matter between active states can be made sterile-sector model independent, apart from N dependence in the lower bound on probability leaking term C αβ [see eq. (D. 2)]. The property holds under the environment of active and sterile neutrino evolution with decoherence in activesterile and sterile-sterile channels, which requires 0.1 eV 2 m 2 J (1-10) GeV 2 for the typical kinematical setting of LBL experiments. It leads to a very simple expression of the oscillation probability in matter, eq. (2.6).
The model-independent nature of the observable is demonstrated by showing that perturbative corrections to eq. (2.6) either vanish or are negligible after averaging over fast oscillations and using large sterile state mass denominator suppression. It is done by an explicit computation to fourth order in W which lefts the unique non-vanishing vacuum term, the probability leaking term C αβ . We have argued by postulating the "Uniqueness theorem" that this feature prevails to all orders in W perturbation theory.
• We have used the oscillation probability formula, eq. (2.6), to analyze ν µ → ν α channels (α = e, µ, τ ), to know in which region of energy and baseline the effect of unitarity violation is large. As a general tendency the effect is sizeable in regions where standard oscillation probability is large, with notable amplification in the two regions corresponding to the solar and the atmospheric MSW enhancement. We have observed relatively large effect in ν µ → ν µ and ν µ → ν τ channels, and pointed out, though qualitatively, that anticorrelation of signals between them would enhance the sensitivity to unitarity violating effects.
• We have discussed the question of how to distinguish low-scale unitarity violation from high-scale one. We have pointed out that outside the region of validity of our Uniqueness theorem, ρE 10 (g/cm 3 ) GeV, the second order W correction [eq. (3.48)] to the leading order could become large and, if detected, it would signal low-scale unitarity violation, offering new way of discriminating between low-and high-scale unitarity violation. Then, it would allow us to probe structure of W matrix elements which bridges between active and sterile sectors. This is to add to the method of detecting the probability leaking term C αβ discussed in [20], which may have a broader applicability by relying on the existence of "sterile", or undetectable but communicable, sector at low energy scales, a generic feature beyond the (3 + N ) model.
Notice that in "constraining mode" of unitarity violation, the model-independence of the framework translates into a universal nature of the bounds, thereby making them more powerful. Whereas in "discovery mode" of unitarity violation, the model dependence, in particular through the W dependent correction terms, is welcome because it serves for identifying the structure of the sterile sector.

JHEP02(2019)015
During the course of this work, we have obtained the new results and had some interesting observations including: • We have obtained an exact solution, eq. (4.15), of S matrix for the Hamiltonian (4. 6) with uniform matter density. It describes neutrino evolution in low-scale unitarity violation in zeroth-order in W , which applies also to the case of high-scale unitarity violation. It has been utilized in section 5 to calculate the oscillation probability in the leading-order as well as its higher order corrections in W . When applied to each shell inside the earth, it could provide a semi-quantitative way of simulating non-unitary neutrino evolution for the terrestrial experiments.
• The value of C αβ , if detected, could reveal structure of the hidden sterile sector. In this paper, this point is illustrated only in a toy model of equally distributed W matrix elements within each flavour, as defined in appendix E. In this model, the probability leaking term scales as 1/N depending upon number of sterile states.
We emphasize that neutrino experiment is the most powerful way to execute leptonic unitarity test in scenarios of low-scale unitarity violation, though it is unlikely in the case of high-scale unitarity violation. Nonetheless, we have to admit that our observations on what we could do for experimental detection of possible non-unitarity effects are rather qualitative to make any definitive claim for possible detection in the future. Clearly, more detailed analyses are called for.
While we worked exclusively on the (3 + N ) state unitary model as a model of lowscale unitarity violation, we do not know if it is the unique choice, or it merely reflects our ignorance. Even in the case there exist more generic class of models for low-scale unitarity violation, the phenomenon of probability leaking is likely to survive. It is because the probability leaking must take place whenever the extra light sector exists and communicates with the three active neutrinos.

A Neutrino evolution equation in flavour basis
The Schrödinger equation takes the form with flavour basis Hamiltonian H in (3.10) where ν a (ν s ) denotes 3 (N ) component vector in active (sterile) space. Apparently, the system depends not only on U and W , but also on Z and V matrix elements, which is not the case in our treatment using the mass eigenstate basis.
Here, we show that the dependence on Z and V is superficial. Since there is no physical meaning of the particular basis for the sterile sector fermions we can redefine it by doing the transformation where Y is a N × N unitary matrix. In the primed basis the Hamiltonian becomes We can arbitrarily choose Y = V † . Then, one can show by using unitarity (3.17) that .
(A.4) Therefore, our system depends only on U and W , and the dependence on Z and V is superficial.

BŜ matrix elements
The method of computingŜ matrix elements is outlined in section 3.3.3 and here, we will collect the results. We denote computed results of the matrix elements ofŜ asŜ[n] to indicate that it is the one that comes from n-th order contribution in H 1 . Since the elements of H 1 are of order either W or W 2 ,Ŝ[n] generally has order W n or higher. To show that a particular contribution is of order W m we use the superscript "(m)". That is, S (m) [n] denotes contribution toŜ that arizes from n-th order perturbative contribution in H 1 and is of order W m . In the following, we will denote

B.1 Contribution toŜ matrix elements from zeroth and first order in H 1
The zeroth and first orderŜ matrix elements can be calculated as follows: The terms above are invariant under generalized T transformationŜ pq (U, W, X, A) → S qp (U * , W * , X * , A * ) [eqs. (3.40)].

B.2 Contribution toŜ matrix elements from second order in H 1
Likewise,Ŝ matrix elements can be calculated in second order inĤ 1 by using the formula for Ω in (3.35) andŜ-Ω relation in (3.36) aŝ iJ,K = f

JHEP02(2019)015
In the expressions above, the combinations of the couplings remain invariant taking the complex conjugate together with p ↔ q while one can verify directly that f (2) pq,r = f (2) qp,r . Hence the expressions are T invariance [eqs. (3.40)].
What we should do in the rest of appendix is to computeŜ matrix elements perturbatively to fourth order in H 1 . In the rest of the appendix, we present only the terms which are required to compute S matrix elements to order W 4 . In view of the relations betweenŜ and S matrix elements given in eq. (3.41),Ŝ IJ ,Ŝ (4) IJ , andŜ (4) iJ (andŜ (4) Ji ) are all unnecessary. We only give the results of manifestly generalized T invariant form ofŜ matrix elements with which it must be straightforward to prove generalized T invariance.

B.3 Contribution toŜ matrix elements from third order in H 1
For the third order terms in H 1 , we havê iJ,kL = f IJ,kL =

JHEP02(2019)015
Notice that forŜ (4) ij [3],Ŝ (3) iJ [3] andŜ (3) Ij [3], the combinations of the couplings remain invariant under complex conjugation together with (p ↔ q) and hence we need that f (3) ab,cd = f (3) ba,cd as can be verified in the expressions above. On the other hand, for S For the fourth order terms in H 1 , we havê Notice that forŜ (4) ij [4], the combinations of couplings remain invariant under complex conjugation with (i ↔ j) and (L ↔ M ) and hence we need that f (4) ij,kLM = f (4) ji,kM L as can be verified from the expressions above. As forŜ (4) IJ [4], the combinations of the couplings remain invariant under complex conjugation with (I ↔ J) and (k ↔ l) and hence we need that f (4) IJ,klM = f We do not display explicitly the expression of each term in (C.1). But, the notation of S (4) αβ [n] diag and S (4) αβ [n] offdiag will be transported to the notation for the oscillation probability such that 2Re S Lk , whose notations are also transported to the oscillation probability.
The oscillation probability to second order in W is given in eq. (3.44) in section 3.4. What is left is, therefore, the expressions of the oscillation probability in fourth order in W , the explicit form of the two terms in (3.52), P (ν β → ν α ) = S Apart from the last line in (C.6) all the terms are suppressed by the two sterile state mass denominators with ∆m 2 Jk which doubly suppress the active-sterile state transition. The first term in the last line is the probability leaking term mentioned in section 2.2.
The nature of each term is explicitly indicated as follows: where S is a 3 × 3 matrix which diagonalizes 1 3×3 − U U † , w is diagonal matrix which consists of eigenvalues of 1 3×3 − U U † , and R is an arbitrary 3 × N complex matrix obeying RR † = 1 3×3 . The construction makes sense for N ≥ 3. Therefore, for a given N there is a large arbitrariness on the choice of the W matrix, and hence on the sizes of the W corrections and C αβ . Lacking a guiding principle of how to choose the R matrix in (D.1), we examine the cases with largest and smallest possible values of C αβ for given values of unitarity violation 1 − 3 j=1 |U αj | 2 (α = e, µ, τ ). It is shown that in the (3 + N ) model C αβ is bounded from above and below as [20] 1 In the (3 + 1) model, the W matrix elements are unique, with the upper and lower bound being equal. For the numbers given in section 5.1, we have W e4 = 0.141, W µ4 = 0.099, and W τ 4 = 0.141 assuming that they are real. Then, the leaking terms have the unique values, C (N =1) eµ = 2 × 10 −4 , C (N =1) µµ = 9.6 × 10 −5 , and C (N =1) τ µ = 9.5 × 10 −4 . The lower bound is realized in the "universal scaling" model described in appendix E, which predicts (J = 4, 5, · · ·, 3 + N ). 23 It is shown in appendix E that under the assumption of equal sterile state masses the universal scaling model predicts the same W 2 correction terms as those of the (3 + 1) model.

E Universal scaling model of N sterile sector
Suppose that we obtain a particular parametrization of U matrix by taking N = 1 sterile sector, as we did in section 5.1. In this (3+1) model, the W matrix elements are completely JHEP02(2019)015 determined, up to phase, by unitarity for a given U matrix (E.1) Now, we attempt to create a toy model of N sterile sector by "universal scaling". We postulate that all the W matrix elements are real and equal: which is consistent with (3 + N ) space unitarity. In this universal scaling model, the order W 2 correction terms in (3.48) remains unchanged provided that we further assume that all the sterile masses are equal. 24 It is because the W matrix elements enter into the W 2 terms in the form where n = 1 or 2. However, the leaking term C αβ becomes smaller by a factor of N in the universal scaling model. In the (3 + 1) model, C αβ takes the largest value, the upper limit in eq. (D.2). Because C αβ is fourth order in W it is evident that in the universal scaling model, which is the lower limit of (D.2).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.