Neutrino mixing, interval matrices and singular values

We study the properties of singular values of mixing matrices embedded within an experimentally determined interval matrix. We argue that any physically admissible mixing matrix needs to have the property of being a contraction. This condition constrains the interval matrix, by imposing correlations on its elements and leaving behind only physical mixings that may unveil signs of new physics in terms of extra neutrino species. We propose a description of the admissible three-dimensional mixing space as a convex hull over experimentally determined unitary mixing matrices parametrized by Euler angles which allows us to select either unitary or nonunitary mixing matrices. The unitarity-breaking cases are found through singular values and we construct unitary extensions yielding a complete theory of minimal dimensionality larger than three through the theory of unitary matrix dilations. We discuss further applications to the quark sector.


INTRODUCTION
Studies of neutrinos are at the frontier of contemporary research in particle physics. These fundamental particles influence crucially processes occurring inside the Sun [1], stars, and supernovae [2,3]. In order to learn about their properties, there are dozens of short-and long-baseline neutrino oscillation experiments studying, e.g., their appearence or disappearence [4]. Thanks to them, we know that at least two out of three neutrinos are massive, though their masses are extremely tiny, at most at the electronvolt level, m ν ∼ O(1) eV [5]. Gathering this information was a highly non-trivial task as neutrino experiments involve the challenge of low event statistics. Among unsolved and important problems in neutrino physics remains the issue of total number of neutrino species in nature. Do we really have only the three electron, muon, and tau neutrino flavors as prescribed by the neutrino theory of the Standard Model (SM) [2,3]? This knowledge is of paramount importance for progress in understanding particle physics and theory beyond the Standard Model (BSM), but also in astrophysics and cosmology of Big Bang, leptogenesis and baryogenesis, and Dark Matter [2,3,6,7]. The point is that additional neutrino species are likely massive, affecting dynamics of many processes and systems, including the Universe as whole. Their existence is theoretically appealing as they could provide an explanation of the smallness of masses of known neutrinos, e.g. by the celebrated see-saw mechanism [8][9][10][11]. There is currently no compelling experimental evidence for extra neutrino states, despite direct collider [12][13][14][15] and indirect electroweak precision studies [16][17][18][19][20][21] providing bounds on their masses and couplings. As a dim clue for their presence one may consider an outcome of the Large Electron-Positron collider (LEP) studies where the central value for the effective number of light neutrinos N ν was determined by analyzing around 20 million Z-boson decays, yielding N ν = 2.9840±0.0082 [22,23]. In fact, a natural extension of the SM by righthanded, sterile neutrinos leads to N ν value less than three [24]. There are also intensive studies concerning sterile eV-scale neutrinos, connected with oscillation phenomena. In the Liquid Scintillator Neutrino Detector (LSND) experiment an excess of electron antineutrinos appearing in a mostly muon antineutrino beam at 3.8σ-level was observed while the SM would predict no significant effect [25]. To solve this puzzle conclusively new experiments are underway [26]. For recent MiniBoone results, see [27]. The question whether sterile neutrinos exist is being researched by ongoing studies performing global analyses of neutrino oscillation data [26,28].
In the description of phenomena like neutrino oscillations, mixing matrices are the central objects. In the SM scenario with three neutrino species, the mixing matrix is known as the Pontecorvo-Maki-Nakagawa-Sakata matrix (PMNS) [29,30]. It is three dimensional, unitary, and can be parametrized by Euler angles. When evaluation of experiments is performed, the hope is from the BSM perspective that an inconsistency in data analysis -in particular violation of unitarity of the mixing matrixwould give a hint for the existence of new neutrino states. In this work, assisted by concepts and theorems taken from matrix theory and convex analysis [31][32][33][34][35][36] we describe an elegant approach to mixing phenomena capable of capturing SM and BSM within the same framework.
At the foundation of our studies lies the analysis of singular values of mixing matrices in form of an interval matrix which gathers knowledge of experimental errors. Firstly, we characterize physical mixing matrices by looking at the largest singular value of a given mixing matrix (which equals the operator norm) and derive on physical grounds that it must be less or equal to unity, a matrix property known as contraction. Using the notion of contractions we consistently stay within the region of physi-cal states with properly correlated mixing elements. Secondly, we study unitarity violation as witnessed by any of the singular values being strictly less than one which has a direct physical consequence and means that the three SM neutrinos mix with unknown species. Therefore identifying such a situation is a smoking gun signal for the existence of additional neutrinos. Finally, we employ the theory of unitary matrix dilations in order to find a unitary extension of any three dimensional mixing matrix which is physically admissible yet not unitary. We apply this method to an example from experimental data and discuss how this approach can be used to find a minimal number of necessary extra neutrino states in a BSM scenario, leading to a complete theory based on a higher dimensional unitary mixing matrix. SETTING We begin our discussion with the SM scenario of three weak flavor -electron, muon, tau -neutrinos. In this framework, mixing of neutrinos is modeled by singleparticle asymptotically free scattering states with a given momentum and spin which are emitted in a fixed flavor state |ν defined by [2] |ν (f ) The PMNS mixing matrix U PMNS is unitary and can be parametrized by [30,37,38] where we denote c ij ≡ cos(θ ij ), s ij ≡ sin(θ ij ) and the Euler rotation angles θ ij can be taken without loss of generality from the first quadrant, θ ij ∈ [0, π/2], and the CP phase δ ∈ [0, 2π]. The current global 3ν oscillation analysis [39,40] gives at 3σ CL These results are independent of the normal or inverse mass hierarchies [41,42] which is not of first concern in this work. The exact ranges can differ also slightly in other analyses [43,44]. In the above, it was assumed that mixing among light and active neutrino states is complete -hence the neutrino mixing matrix is unitary. However, the situation can be more complicated. In a BSM scenario other neutrino mass and flavor states can be present that we denote by | ν (m) j and | ν (f ) j for j = 1, . . . , n R , respectively. In this scenario mixing between an extended set of neutrino mass states {|ν Such block structures of the unitary U are present in many neutrino mass theories. Note, that (4) effectively implements an assumption of unitary mixing restricted to the level of single particle states only, e.g., neglecting interaction effects which are expected to be weak. Indices "l" and "h" in (4) stand for "light" and "heavy" as usually we expect extra neutrino species to be much heavier than known neutrinos, cf. the see-saw mechanism [8][9][10][11]. However, it does not have to be the case: They can also include light sterile neutrinos, which effectively decouple in weak interactions, but are light enough to be in quantum superposition with three SM active neutrino states and to take part in the neutrino oscillation phenomenon [45]. The observable part of the above is the transformation from mass |ν If V is not unitary then there necessarily is a lightheavy neutrino "coupling" and the mixing between sectors is non-trivial V lh = 0 = V hl . Without extra states | ν (m) , we end up with situation described in (1)-(3), V → U PMNS , i.e., there are either no BSM neutrinos or they are decoupled on the level of the joint mixing matrix.

PHYSICALLY ADMISSIBLE MIXING MATRICES ARE CONTRACTIONS.
In this section we will make precise the notion of physically admissible mixing matrices. To this end, we will study the singular values σ i (V ) of a given matrix V which are equal to the positive square roots of the eigenvalues λ i of the matrix V V † , i.e. σ i (V ) ≡ λ i (V V † ) for i = 1, 2, 3 [31]. Singular values generalize eigenvalues to all kinds of matrices, e.g. those not diagonalizable by a similarity transformation or even rectangular ones, and have useful properties that in particular can be related to the operator norm V ≡ max i σ i (V ). In the SM scenario one would only consider unitary matrices, hence V = 1 and all singular values are equal (see the Appendix A). In this work, we are also interested in constraints on V as a principal submatrix of a unitary U realizing some BSM scenario (4). For any such matrix V the operator norm is bounded by unity a matrix property known as contraction. In other words, if U is unitary then U = 1 and for any submatrix V of U it holds that V ≤ U = 1, see the Appendix A for a simple proof. Observe that V = 1 is not sufficient for V to be unitary and any significant deviation of any singular value from unity σ i (V ) < 1 signals BSM physics. Physically, measuring a mixing matrix with nonunit singular values means that a given neutrino mixes with other ones that are not being observed and hence the unitarity loss. Note that any observable mixing matrix must be a contraction both in the SM and BSM scenario. Moreover, singular values are suitable quantities while working with experimental data since they are stable under addition of perturbing error matrices and the resulting errors of the operator norm can be upper bounded, while the stability of eigenvalues can be in general very weak, e.g., violating Lipschitz continuity [31].
It can be achieved only if matrices after the perturbation remain normal [31], a condition that obviously cannot be fulfilled generally when considering experimental data.
In this work, we show how the contraction property allows to distinguish physically admissible mixing matrices. Namely, recall that if V is a submatrix of some larger unitary mixing matrix U then it must be a contraction. Conversely, as we will show presently any V which is a contraction can be completed into a unitary mixing matrix U whose minimal dimension can be read-off from the singular values of V . Hence we establish the following characterization useful for data analysis allowing to decide whether a candidate mixing matrix V is physically admissible.
Definition 1 (Physically admissible mixing matrix). A matrix V is a physically admissible mixing matrix if and only if it is a contraction i.e. V ≤ 1.

INTERVAL MATRICES, UNITARITY VIOLATION AND CONTRACTIONS.
Though the matrix U PMNS is unitary, information on BSM physics can be hidden there. To see this one should ask what would be the result of a fit assuming unitarity in the case that the mixing was actually non-unitary? In a BSM scenario (5) a unitary fit to (2) would hide the BSM physics in the error bars and hence the experimental Euler angle ranges may reflect not only measurement inaccuracies but also the hypothetical non-unitarity of the underlying mixing. For similar reasons, the search for BSM via unitary triangle analysis [46,47] is based on PMNS data. So far experimental analysis is not precise enough to confirm or exclude definitively BSM [5].
In order to find the non-unitary cases, we discretize the experimentally allowed ranges in (3), calculate the corresponding U PMNS matrices using (2) and collect the extreme values of each matrix element that occurred into an interval matrix V osc [39] UPMNS → Vosc = (7)  We will write V ∈ V osc whenever all entries of V lie in the intervals of (7). This interval matrix is real as we have fixed for simplicity in (3) the CP-phase δ to be zero but our analysis can also be applied to complex mixing matrices. The exact values in this interval matrix can differ slightly depending on global fits and considered approaches [39,43,44,48]. Our construction of the interval matrix is based on [39,40] where the interval matrix is obtained in the same way, i.e. by looking at the extreme values of the entries of the mixing matrices V ij for all possible Euler angles consistent with the oscillation data. As an alternative, we will also refine this procedure by looking at convex combinations of U PMNS matrices which should be even closer to the data by retaining correlations between matrix elements. In particular this allows to construct candidate BSM matrices as toy examples to study various methods related on mixing matrices close to the data. It is important to observe that it is not necessary to construct (7) from U PMNS . In principle, such an interval matrix could be derived directly from experimental data. In the neutrino sector direct experimental access to each of the entries of the 3 × 3 matrix individually is presently not possible and experimental analyses based on U PMNS are a natural choice. If this was possible then the interval matrix would become a useful way of bringing various experimental findings together. Indeed approaches to oscillation analysis independent of PMNS are possible [48].
Here α is a lower triangular matrix and U is unitary. Such a triangular structure of α is especially convenient for singular values analysis [60]. This parametrization is often used in oscillation analysis, e.g. [56,61,62]. We discuss both α and η parametrizations in the wider context of matrix analysis in the Appendix E.
Observe that decompositions given only by (8) or (9) need some extra conditions to produce contractions exclusively as in general it can happen that although a given matrix lies within experimental limits and is of the form given by Eq. (8) or (9), it will not be a contraction (for a proof see Proposition 1 below). In particular, such a condition can be provided by embedding a three-dimensional mixing matrix into a larger unitary one. Accordingly, it is standard in the neutrino unitarity-breaking literature to take one of the approaches (8), (9) together with that embedding as the precise definition of the so-called α or η parametrization (cf. [50,51,54,56,58,59,61,63,64]). Therefore, by combination with such additional conditions the contraction property of the mixing matrix is secured, see Appendix E for further discussion. However, it is common to present the data of such analyses in form of an interval matrix, where the correlations between elements are lost. If one would like to consider a mixing matrix η taken from such interval matrix as a point of departure, it is profitable to find a condition on that particular η matrix to be physical, i.e. to give rise to physically admissible mixing matrix. For this, Proposition 2 characterizes a particular sufficient condition securing any candidate mixing matrix taken from an interval matrix to be physically admissible. A similar argument can be proven for the α parametrization (9) and in fact, it has been shown in [56] that (8) and (9) are equivalent. Therefore only η case is regarded in the following.
(10) Then Θ contains matrices which are not contractions and hence are not valid mixing matrices i.e are unphysical.
Proof. Let V = (1 + η)U ∈ Θ satisfy η = 0. If η has a positive eigenvalue λ + > 0 then we use that I + η is diagonalizable and obtain V = I + η ≥ 1 + λ + > 1 using unitary invariance of the operator norm and the lower bound comes from the fact that there may be other eigenvalues that are still larger than λ + . If η 0, i.e., it has no positive eigenvalues, then we observe thatη = −η has at least one positive eigenvalue and the constraints of Θ are satisfied. Thus we find forṼ = (1 +η)U ∈ Θ that Ṽ > 1.
As an example, let us consider an interval matrix for η (see [50]) and its particular elements The subscript "max" indicates that we have chosen elements of η to have largest absolute values given the constraints of the respective interval matrix. This matrix is hermitian, as given in (10). As follows from Prop. 1, more stringent limits, e.g. [64,65] do not change the situation. Due to the fact that we bound only absolute values we consider the following two cases Performing a singular value decomposition [31] we obtain the following singular values We can see that both spectra, which correspond to the singular values of the matrices V = (I ± η max )U , contain eigenvalues larger than one which means that mixing matrices V constructed using this particular matrices are not contractions. We consider a general form of the η parametrization, where U ≡ Q in (8) is an arbitrary unitary matrix. Observe that, regarding our analysis of η max taken from Eq. (11), there is a subtle detail. To check whether a matrix V is a contraction, we do not use the unitary matrix U at all. This follows from the nature of the η parametrization which is in fact a polar decomposition. The contraction property is based on the operator norm which is unitarily invariant, which means that only the polar matrix contributes and the unitary part (by definition) does not change the norm. Thus for the analysis of singular values that we have done, the unitary part is irrelevant. It should be no surprise that such particular element η max could be unphysical in spite of its hermiticity, since the very construction of interval matrices destroys correlations between elements as discussed above. Nevertheless, we can restrict ourselves to physical matrices which are contractions by the following Proposition 2. If all ij ≤ are sufficiently small then restricting to negative semi-definite perturbations η 0 yields exclusively physically admissible mixing matrices (14) For any V = (1 + η)U ∈ Θ we have that the norm of V can be obtained from the largest eigenvalue of the diagonalizable matrix 1 + η.
Proof. It suffices that < 1 n where n is the dimension of the matrix η we will find η < 1. As the identity I and η are simultaneously diagonalizable and all eigenvalues λ i (η) of η are non-positive we find that V = I + η = 1 + max i λ i (η) ≤ 1 so all V ∈ Θ are contractions and thus are admissible mixing matrices.
We sum up this section in the following way. As mentioned already, there are parametrizations which allow to generate U PMNS -like 3 × 3 matrices which by construction are contractions [50,51,54,56,58,59,61,63,64], respecting the present experimental bounds. If not secured directly, the condition of negative semi-definite perturbation (14) can be used to ensure that the considered mixing matrices are physically admissible when working particularly with (8). In general, it is numerically efficient to check directly the contraction property (6) of examined mixing matrices for any parametrization.
We proceed by characterizing physical mixing matrices consistent with the experimental data. Firstly, let us note that the set of all (unrestricted) contractions is a unit ball in operator norm and hence is convex. This abstract property allows to describe B in terms of its extreme points which in the case of contractions are unitary matrices U 3×3 [66]. In fact, we can easily find that a convex com- For such combinations, restricted to experimentally determined U PMNS unitary matrices, we have V ∈ V osc because the interval matrix is constructed from extreme values of U i and convex combinations cannot change these bounds. Conversely, when V ∈ V osc is a contraction but cannot be written as a convex combination of unitary matrices within allowed angle ranges (3) then it means that the construction of V osc through extreme matrix elements simply introduced discrepancies with the data by disregarding correlations between matrix elements. Therefore, the set of all finite convex combinations of PMNS matrices given by comprises all contractions spanned by the experimental data, see Fig. 1. This definition takes into account possible non-zero values of the CP-phase δ. Currently, it is not possible to measure experimentally values of all elements of the neutrino mixing matrix in the three dimensional flavor space [48]. To determine missing elements one uses Euler angles obtained from available data and calculates unreachable matrix elements of the neutrino mixing matrix by (2). The set Ω could be explored in future in a broader context for data analysis and independent cross checks with experiments that measure entries of the 3 × 3 mixing matrix directly rather than through Euler angles. The matrices in Ω with M = 1 yield admissible PMNS matrices, while taking M ≥ 2 a a a U 1 Illustration of the neutrino mixing space. Eq. (6) states that physical mixing matrices V ∈ Vosc lie within an abstract operator norm unit ball represented by the ellipse. On the left are cases that are physically admissible, but are excluded by the experimental data (3). The middle region Ω represents relevant mixing matrices consistent with the experiment which are convex combinations of unitary PMNS matrices. The cord slicing Ω consists of convex combinations of two PMNS matrices U1 and U2, e.g., V = 1 2 U1 + 1 2 U2 which is further discussed in section 6. The rectangle on the right depicts the interval matrix form of the data Vosc that is largely unphysical and may include contractions spanned by unitaries outside of (3).
allows to obtain non-unitary contractions. Although the upper limit M in (15) is not unique, in principle it can be bounded from above by Carathéodory's theorem, which states that if a point x ∈ R n lies in the convex hull of some set X, then x can be written as a convex combination of s-many points from X such that s ≤ n + 1 [67].
Since matrices under study (in the CP-invariant case) are points in R 9 , elements of (15) are narrowed down to combinations of at most M = 10 unitary U P M N S matrices. Thus one obtains an upper bound for the number of free parameters under study in this approach. From the point of view of particle physics phenomenology (but also optimization theory), it would be interesting to refine M even further and look for the smallest possible M , called the Carathéodory number, that would allow to span Ω (see e.g. [68]). While certainly important, this issue goes beyond the present study.
In the next section we give an example of two unitary PMNS matrices within the accuracy of the interval matrix. This exemplifies how to find non-unitary BSM cases within admissible set Ω, through the analysis of singular values. For BSM mixing matrices it is possible to find minimal model extensions consistent with the data again using singular values. A unitary dilation is an operation that extends a matrix which is a contraction to a uni-tary matrix of an appropriate dimension. Our approach to find a unitary dilation of possible smallest dimension employs the special case of cosine-sine (CS) [32] decomposition of unitary matrices as follows. It can be proven that any unitary matrix U ∈ M (n+m)×(n+m) (C) can be brought to a canonical form We use this result to extend any contraction V ∈ Ω to a unitary matrix. First, we find a singular value decomposition of V , i.e. V = W 1 ΣQ † 1 where W 1 , Q 1 are unitary, Σ comprises the singular values σ i (V ) and is diagonal. Next, we determine the number r of unit singular values defining I r and collect the rest into a diagonal matrix C. This yields Σ = I r ⊕C. Finally, we define S = √ I m − C 2 and choose W 2 , Q 2 to be arbitrary unitaries of appropriate dimension. Conjugating the CS matrix constructed in that way by W and Q yields the unitary dilation U of V . Below an example of a non-unitary contraction V ∈ Ω with m = 2 will be discussed, extended into a unitary matrix U of dimension 5. Any larger unitary dilation of V can by obtained by the general form of CS decomposition, see Theorem 3 in the Appendix D. There, we also prove that m, also known as the dimension of the defect space, is the minimal number of new neutrino species necessary to ensure unitarity. To obtain this number, one thus has to take experimental errors into account. Assuming that the data V includes an error matrix E and is of the form V + E, we can establish stability of the defect space. We use Weyl inequalities [69,70] for decreasingly ordered pairs of singular values of V and V + E which read In our case E should be taken as E ij ∼ 0.001 by (7), and the uncertainty in the precise value of singular values is bounded by E = 0.003. Note that this criterion applies both to the selection of contractions from the full interval matrix (7) and to determination of minimal dimension of matrix dilation.
To show dilation procedure in action we restrict all matrix elements to real numbers, hence the complex phase δ is equal to zero and thus we work with orthogonal matrices. The first step is to pick a contraction from the convex hull Ω (15). As an example let us consider two unitary matrices obtained from the experimental ranges The chosen convex combination will be constructed as a sum with an equal contribution of above matrices In order to make use of CS decomposition and parametrize the unitary dilation U of the matrix U 11 ≡ V , first we have to find its singular value decomposition where We will parametrize only the most interesting case of unitary dilation of a minimal dimensionality, hence of a minimal number of additional neutrinos i.e. the number of singular values strictly less than 1. Since the matrix Σ determines singular values of V , this number equals 2. Hence it is possible to construct unitary dilation U of the minimal dimension 5 × 5.
To complete the construction, we are left only with two free unitary 2 × 2 "parameters" W 2 and Q 2 , for the sake of this example we choose them randomly Having all ingredients and making all necessary calculations we find the following form of the unitary dilation of V given by U = ( The matrices U, U differ by both off-diagonal block and by the bottom-diagonal block. However, the scale of the off-diagonal block is comparable in both cases. The reason of this lies in the fact that to construct each of these blocks we use C and S fixed by the singular value decomposition of V matrices. The biggest difference can be observed in the bottom diagonal block since only the matrix S is fixed in both cases. However the global scale of each block (global in a sense of the Frobenius norm which is a entrywise norm defined in [31]) is conserved in each of these cases. Since this norm is unitarily invariant, the choice of W 1,2 and Q 1,2 does not change its value.
The dilation procedure described above is based exclusively on mixing matrices. In contrast, there are constructions in literature which refer in addition to the mass spectrum, see for instance [64,71,72]. In the approach taken in this work, the information on the number of additional neutrinos, i.e. the dimension of complete unitary mixing matrix, is nicely seen through the number of nonunit singular values. As discussed in the Setting section, our approach based on singular values and the dilation procedure is general, no matter if extra neutrino states are heavy (e.g. seesaw mechanisms) or light and sterile. As far as the present situation in neutrino physics is concerned, the minimal 3+1 neutrino scenario is still not excluded, though LSND and MiniBoone results make it a less probable scenario. For global analysis see [26,28,73]. Here we considered an example of extension to 3 + 2 dimension. However, it is still possible to find one of singular values strictly less than 1 while remaining two are The interval mixing matrix (7) contains unitary (SM), non-unitary contraction (BSM) and unphysical matrices (the latter have to be discarded by the contraction property). We have found that statistically about 4% of matrices V ∈ V osc are contractions, while some unphysical ones have norms as large as V = 1.178. This result was obtained within 0.003 accuracy, by uniformly sampling elements of the intervals of V osc with sufficiently high statistics. All calculations presented in this work has been made in Mathematica [74]. The statistical analysis of distributions of contractions in V osc has been performed under the assumption of mixing parameters errors having a uniform distribution. This implies that also values in V osc were treated uniformly. A discretization of intervals in V osc was made with a step 0.001 to match the precision of extreme values. Up to 10 9 randomly generated matrices have been produced within V osc ranges for which singular values have been found. Next, the largest singular value for each random matrix was compared to the number 1.003, to be consistent with precision ensured by the stability of singular values, splitting this way matrices into two sets of contractions and non-contractions.
Likewise, we analyze distributions of contractions for a given element in V osc . We fix a value of one of the elements of V osc and then randomly generate matrices and make the same analysis as above. As an illustration, Fig. 2 presents contraction distributions for two exemplary matrix elements taken from the full interval matrix. While one may argue that these diagrams show only statistical density, Proposition 4 in the Appendix C shows that there is a sharp matrix boundary (surface in C 9 ) with interior composed solely of contractions.
If we shrink errors in (3) to 1σ C.L., we get 11% of contractions, instead of 4% discussed at the beginning of this section. Narrowing angle ranges (3) usually increases the amount of contractions in V osc , however for arbitrary angle ranges this does not always occur. Concerning new physics, it is interesting to ask how strong a contraction can be found in V osc . The minimal value of the norm for V ∈ V osc is V min = 0.961 and can be obtained by sufficiently fine discretization of V osc . Alternatively, V min can readily be obtained by semi-definite programming which is a very useful numerical tool when analyzing properties of interval matrices [75].
It should be stressed that we apply our methods to data in order to illustrate our matrix machinery in applications to neutrino mixing matrices but do not attempt to make a definite analysis. We have made rough estimations based on a construction where experimental data and PMNS formalism are used, though as mentioned already, the interval matrix can be obtained even directly without restriction to PMNS parametrizations when nonunitarity is assumed from the very beginning [48]. This interesting and universal option is left for separate and detailed future studies.

QUARK DATA ANALYSIS.
Our scheme in Fig. 1 is general enough to be used in the quark sector as well. For quarks the unitary CKM mixing matrix [37,76] can be parametrized in the same way as the PMNS mixing matrix for neutrinos Experimental results have established the following hierarchy of mixing parameters s 13 s 23 s 12 1.
Due to this order it is convenient to present V CKM in an approximate parametrization proposed by Wolfenstein [77] which reflects the above hierarchy. The mixing parameters (27) are connected with Wolfenstein parameters in the following way This results in the following structure of the quark mixing matrix We are interested how contractions are distributed within V CKM with respect to experimental values of the mixing parameters [42] λ = 0.22506 ± 0.00050, where,ρ = ρ(1 − λ 2 /2) andη = η(1 − λ 2 /2). The application of above results to (29) gives us the following experimental intervals for elements of the mixing matrix Let us recall that in the neutrino case, minimal and maximal deviations from unity are 0.961 and 1.178, respectively. It shows how much as far as precision of analysis in the neutrino sector is concerned still must be done there.
It is interesting that a vanishing fraction of matrices within V CKM has norm strictly less than one. This can be a sign that the only contractions in the quark sector are unitary matrices. However, since we have used only leading order of the mixing matrix expressed by the Wolfenstein parameters, additionally more refined analysis of this sector is necessary. In principle, we do not have to rely on the Wolfenstein parametrization and the analysis can be done directly on quark data in form of an interval matrix. At the LHC there are already direct measurements of V tq (q = d, s, b) by studying top production as well as its decays and charge asymmetry [78][79][80]. Our approach based on interval matrix will become very interesting in the context of future collider experiments, like FCC-hh, with center of mass energies few times larger than LHC [14], where all elements of the interval mixing matrix can be probed directly with much better precision.

SUMMARY AND OUTLOOK.
We have shown how to recover physically admissible mixing matrices from the interval matrix representation of neutrino or quark data, namely any contraction matrix within the interval matrix is physical and has properly correlated matrix elements. This characterization is complete, as any contraction can be completed into a unitary matrix via a unitary dilation procedure which yields an extension of minimal dimension. The approach is universal in the sense that it does not invoke any specific parametrization and is based on general features of the interval matrices. Physical mixing matrices consistent with the experiment are shown to have the structure of a convex hull over admissible PMNS matrices.
Singular values play a special role in our analysis. The general observation is that whenever we find singular values smaller than one, it is a signature of BSM. This observable seems to be an interesting alternative to other quantifiers of unitarity breaking so far employed in literature. We are commenting on possible analysis in the quark sector and our estimations based on Wolfenstein parametrization point out very little space for nonunitary effects there.
Finally, assuming a BSM scenario we show how to construct a unitary mixing matrix of minimal dimension larger than 3 consistent with data. It allows us in particular to construct a dilation procedure to determine the minimal number of extra neutrino species, compatible with experimental data in a BSM scenario. This is potentially a very fertile area of studies. Should a BSM signal be found, dilation theory will be a promising point of departure for further analysis. Of course our studies are not complete with this commencing paper. The estimation of errors to judge unambiguously on deviations of singular values from unity will be crucial in future. In this work we estimate errors on singular values through Weyl inequalities.
Our methods are based on advanced matrix analysis, studying the singular values of mixing matrices. We apply a model-independent analysis based on the interval matrix to the present data, in a way that may become significant in future experiments that will measure entries of this matrix directly. It can be also useful through Propositions 1 and 2 to cross-check with other analyses based on specific parametrizations, since the contraction condition is easy to apply.
We shall go further in this direction and merge our studies on mixings (eigenfuncions problems) with masses (eigenvalues). For instance we could study the angle between subspaces of the mass matrices to connect neutrino masses with their mixing. This approach is closely related to the methodology presented in our work. Moreover, a separate analysis of the properties of the neutrino mass spectrum could be done exclusively. For this we might adopt many advanced methods of matrix analysis, e.g. Gershgorin circles. Clearly further potential for practical applications of our procedures is there.

ACKNOWLEDGEMENTS.
This publication was made possible through the support of a grant from the John Templeton Foundation (Grant No. 60671) and the support of the Polish National Science Centre (NCN), Grant No. DEC-2013/11/B/ST2/04023.

APPENDIX.
In the following we give the technical details supplementing the results of the main text. We begin in the Appendix A by providing more details on contractions as principal submatrices of unitary matrices. Then in the Appendix B we give a very simple example of how contractions allow to restrict parametrizations of mixing matrices. We then provide a section in the Appendix C describing the interval matrices within convex geometry. In the Appendix D we provide a description of the theory of matrix dilations. In the Appendix E various non-unitary parametrizations are classified. Their relation to contractions is discussed.

A. CONTRACTIONS.
A matrix norm is a function · from the set of all complex matrices into R that satisfies for any A, B ∈ M n×n the following properties In other words, a matrix norm is a vector norm (first three conditions in (33)) with additional condition of submultiplicativity. The most important norm in our work is the operator norm A = max x =1 Ax for which one can prove that it is equal to the largest singular value i.e. singular values are the positive square roots of the eigenvalues of AA † denoted by λ i (AA † ). We note that there exist other matrix norms that bring different properties into focus [31] but are less important for mixing matrices.
We now consider any principal submatrix V of a unitary matrix U and show that it is a contraction, i.e., V ≤ 1 in operator norm.
Proof. It is straightforward to see that for any unit x ∈ C m there is a unit embedding y ∈ C n of x such that (namely, by inserting zeros at entries of y corresponding to columns of A deleted to obtain B). Furthermore, the range of this embedding is a subspace of C n , hence which gives the result. The next observation is almost trivial, yet is crucial in the analysis of neutrino mixing matrices in the main text.
Corollary 1. Let U ∈ M n be unitary. Then U = 1 and any submatrix V of U is a contraction.
Proof. The equality λ i (U U † ) = λ i (I) implies that U = 1. By Proposition 3, for any submatrix V of U it holds that V ≤ U = 1, hence V is a contraction.

B. UNITARITY AND CONTRACTIONS: TOY EXAMPLE.
Here we provide more details on problems occurring when studying non-unitary U PMNS through a particular parametrization. For U PMNS it holds that sum of probability of neutrino oscillations equals 1 α P iα = 1, e.g. P ee + P eµ + P eτ = 1.
However, for a non-unitary U analogous relation is not fulfilled. Let us see it in a simple case of two flavors (the same can be done for dimension 3 modified U PMNS matrix), when U is defined as (Θ 2 = Θ 1 + ) In this case we get, We can see that the sum can be either larger or smaller than 1. This example was given in [81], however, no clue in that time was how to interpret possible results where sum of probabilities does not equal 1. Here we show that matrix (38) is not the right way to parametrize BSM effects. Let us find the norm which helps to interpret the matrix (38).
First, we calculate U U T and U T U for (38), s(c) a ≡ sin(cos)Θ a As for the real A we have A T A = AA T = A 2 , we can focus only on one of these products. We write U U T in the following form This can be simplified into where s 3 ≡ sin Θ 3 = sin(Θ 1 − Θ 2 ). B is symmetric and its eigenvalues are equal ±s 3 . Let V be a unitary matrix such that V T BV = D = diag(s 3 , −s 3 ). Since operator norm is unitarily invariant [31], we write Since I + D equals its operator norm, i.e., the largest singular value is equal So we can see that by adding B to identity matrix we cannot decrease operator norm Thus U ≥ 1.
As discussed in the main text, physically meaningful theory should include only fields for which contraction relation U ≤ 1 is fulfilled, and U > 1, being a part of some more complex complete theory based on unitarity cannot describe BSM effects at all. The result (49) implies that not all parametrizations which violate unitarity are a proper choice, and a toy mixing matrix (38) is superfluous from physical point of view. It fulfills U = 1 for = 0, but then trivially unitary matrix is recovered.

C. CONVEX GEOMETRY AND INTERVAL MATRIX ANALYSIS.
Here we gather necessary facts and comments that refer to convex geometry, which plays a crucial role in the paper in two-fold way: it gives a very convenient parametrization of contraction matrices (see Theorem 1) and provides some decisive conditions on distributions of (non-)contractions in interval matrices (see Proposition 4).

Definition 2. [35] A nonempty set
of all convex combinations of finitely many points of A.
The following theorem states that there is an analogue of a linear span in convex geometry, such that the span is over all extreme points of the set A, i.e points that are not interior points of any line segment lying entirely in A.

Theorem 1. (Krein-Milman) [36]
Let X be a topological vector space in which the dual space X * separates points. If A is compact, convex set in X, then A is closed, convex hull of its extreme points.
Proposition 4. Once a set of matrix contractions is given, the convex hull with vertices at this set contains only contractions.
Proof. Let n be fixed and consider the non-empty polytope P = × n 2 k=1 [a k , b k ]. To every p ∈ P we associate a matrix A (p) with entries A (p) i,j = p ζ(i,j) where ζ : [n] 2 → [n 2 ] is the bijective map defined by ζ(j, k) = (j − 1)n + k. We will show that the subset of matrices based in P is convex i.e. for p, q ∈ P if A (p) ≤ 1 and A (q) ≤ 1 then A (λp+(1−λ)q) ≤ 1 for 0 ≤ λ ≤ 1.
We now explicitly calculate = λA From the triangle inequality we obtain This means that if one verifies that for a set of points p 1 , . . . , p N the matrices are contractions then for all matrices in the convex hull p ∈ conv{p 1 , . . . , p N } the matrix A (p) will be a contraction.
In particular this means that if one checks that contractions are vertices of some Q = × n 2 k=1 [a k , b k ] ⊆ P then no matrix inside Q will have norm larger than 1.

D. UNITARY DILATIONS.
To find a complete theory for BSM mixing matrices we need to find a matrix that has a non-unitary V as a principal submatrix and is unitary. In 1950 Halmos [33] noticed that any contraction A acting on a Hilbert space H can be dilated to a unitary operator which acts on H ⊕ H space by Few years later Sz.-Nagy [34] generalized this idea. In the Halmos construction we see that for a n × n matrix A its unitary dilation U will have dimension 2n × 2n.
There exists a further theorem [83] which allows to dilate a contraction to a unitary matrix of possibly lower dimension than 2n, yet some additional conditions must be satisfied.
Recall that rank of a matrix can be defined as the number of its non-zero singular values. We use this theorem to show that m is optimal.
Corollary 2. Let A, U be as above and m = rank(I − A † A). Then the minimal dimension of U is n = k + m.
Proof. Suppose n < k + m. From Theorem 2 we have m ≤ min{k, n − k}, hence m ≤ n − k in particular. Thus n ≥ k + m which contradicts the assumption.
In the main text we construct the minimal extension and use the fact that rank(I − A † A) is equal to the number of singular values of A strictly less than one, which is a direct consequence of a rank definition given above. The construction is achieved through the cosine-sine (CS) decomposition of unitary matrices. In [32] it has been shown how that Halmos construction (55) is a particular example of the CS decomposition. This construction in its generality allows for dilations of dimension determined by Corollary 2. Again, singular values play a crucial role here. If m ≥ n, then there are unitary matrices W 1 , Q 1 ∈ M n×n and unitary matrices W 2 , Q 2 ∈ M m×m such that where C ≥ 0 and S ≥ 0 are diagonal matrices satisfying C 2 + S 2 = I n .
If n ≥ m then it is possible to parametrize a unitary dilation of the smallest size.
Corollary 3. The parametrization of the unitary dilation of smallest size is given by where r = n − m is the number of singular values equal to 1 and C = diag(cos θ 1 , ..., cos θ m ) with | cos θ i | < 1 for i = 1, ..., m.

E. BSM PARAMETRIZATIONS OF NEUTRINO MIXINGS AND CONTRACTIONS
There exist three different matrix factorizations that decompose matrix into product of two matrices of which one is unitary, namely [31,84] 1. Polar decomposition, 2. QR decomposition, 3. Mostow decomposition.
Two first are used frequently in neutrino physics in a context of parametrization of non-unitarity effects in the neutrino mixing matrix. These are the polar decomposition and a modified version of the QR decomposition. Thus, let us take a closer look at these two parametrizations. The polar decomposition factorizes a given square matrix A into the following product where matrix P is a positive semidefinite Hermitian matrix and U is a unitary matrix. The polar factor P is uniquely determined and is given by √ AA † while unitary part is also uniquely determined if the initial matrix is non-singular.
To our knowledge an application of the polar decomposition to parametrize a deviation from unitarity in the neutrino sector appears for the first time in [51]. There, the polar factor is further decomposed in the following way where a matrix η describes deviation from unitarity of the neutrino mixing matrix. As we recall from the main text, physical mixing matrices must be contractions, i.e. matrices with spectral norm less or equal to one or equivalently with the largest singular value less or equal to one. Let us notice that in general the polar decomposition do not provide this property. To see this let us look at a simple example, where we take matrix η in simple diagonal form 1 where 0 < ≤ 1.
Observe that this results in a positive semidefinite matrix P = I − η which is necessary for a polar factor. However, such P is not a contraction since one singular value will be always larger than one, independently how small is.
Recently, the polar factor in a form of (60) was identified with a matrix I − ΘΘ † 2 which arises in the context of the complete unitary mixing matrix [54] (for similar construction see also [63,64]). Thus in a scenario such that the complete unitary mixing matrix is considered, the polar factor is from definition a contraction. In this approach, to ensure polar factor I − η be a contraction, a necessary condition for the matrix η follows in a form of a positive semidefinite matrix. Using a fact that the operator norm is unitarily invariant, it can be shown that for sufficiently small entries of the matrix η also the inverse is true, i.e, if the matrix η is positive semidefinite then P = I − η must be a contraction. Scenarios that employ such unitarity-breaking constructions are usually called top-down approaches.
The second of currently used factorizations in neutrino physics is the QR decomposition. It factorizes a given matrix into product of a unitary matrix Q and an upper triangular matrix R and was proposed as a parametrization of the neutrino mixing matrix in [58,59]. For this purpose a modified version of the QR factorization is used, namely the LQ decomposition, where L corresponds to a lower triangular matrix and Q is a unitary matrix. Moreover, in the context of the neutrino mixing this lower triangular matrix is further split into the following form where the matrix α is lower triangular and describes a deviation from unitarity of the U P M N S . Recently, a correspondence between the polar and QR parametrizations in the case of neutrino mixing was found [56].
In the end let us look briefly at the last factorization, i.e. Mostow decomposition. It decomposes any nonsingular complex matrix A in the following way where U is a unitary matrix, K is a real skew symmetric matrix and S corresponds to a real symmetric matrix.
To this point we discussed matrix decompositions commonly used to parametrize a possible deviation from unitarity of the mixing matrix. Currently they are mostly used in a top-down analyses [50,51,54,56,58,59,61,63,64] which means that they are considered as a part of a complete unitary matrix each time. As we shown, such an approach trivially ensures contraction property for these matrices. Let us note that top-down parametrizations are based on general treatment of unitarity breaking effects described by matrix factorization and there is lack of exact description based on entrywise parametrization of mixing matrix which would fulfill automatically contraction property. Such a construction would be very useful. So far, parametrizations which are constructed fulfill a condition of contractions involving a general I − ΘΘ † 2 representation of the matrix η, parametrizing a matrix Θ in such a way that ΘΘ † 2 will fit into currently known limits on η.
Actually in our strategy we come back to the bottomup scenario, as our analysis starts from the present state of knowledge on U PMNS mixing data in a form of an interval matrix and we examine directly whether matrices within are physically meaningful (i.e. are contractions). An extension of this idea allows us to define complete region of physical mixing matrices as a convex hull of U PMNS matrices which ensures that any physical mixing matrix can be constructed as a convex combination of U PMNS matrices. Now, let us emphasize the relation of our approach to the polar decomposition. In our analysis we use singular values as an indicator whether a given matrix is a contraction. However, it is known that eigenvalues of the polar factor, which follows from the definition, are equal to singular values of an initial matrix. Thus from that perspective a polar decomposition can be treated as a compact version of singular value decomposition. Nevertheless, from the numerical analysis perspective singular value decomposition algorithms are more natural, since they arise from eigenvalues decomposition of matrices AA † and A † A. Thus in most cases, in order to obtain an algorithm for a polar decomposition we have to translate algorithms for the singular value decomposition.
iciently small entries of the matrix η also the inverse is true, i.e, if the matrix η is positive semidefinite then P = I − η must be a contraction. Scenarios that employ such unitarity-breaking constructions are usually called top-down approaches.
The second of currently used factorizations in neutrino physics is the QR decomposition. It factorizes a given matrix into product of a unitary matrix Q and an upper triangular matrix R and was proposed as a parametrization of the neutrino mixing matrix in [58,59]. For this purpose a modified version of the QR factorization is used, namely the LQ decomposition, where L corresponds to a lower triangular matrix and Q is a unitary matrix. Moreover, in the context of the neutrino mixing this lower triangular matrix is further split into the following form where the matrix α is lower triangular and describes a deviation from unitarity of the U P M N S . Recently, a correspondence between the polar and QR parametrizations in the case of neutrino mixing was found [56].
In the end let us look briefly at the last factorization, i.e. Mostow decomposition. It decomposes any nonsingular complex matrix A in the following way where U is a unitary matrix, K is a real skew symmetric matrix and S corresponds to a real symmetric matrix.
To this point we discussed matrix decompositions commonly used to parametrize a possible deviation from unitarity of the mixing matrix. Currently they are mostly used in a top-down analyses [50,51,54,56,58,59,61,63,64] which means that they are considered as a part of a complete unitary matrix each time. As we shown, such an approach trivially ensures contraction property for these matrices. Let us note that top-down parametrizations are based on general treatment of unitarity breaking effects described by matrix factorization and there is lack of exact description based on entrywise parametrization of mixing matrix which would fulfill automatically contraction property. Such a construction would be very useful. So far, parametrizations which are constructed fulfill a condition of contractions involving a general I − ΘΘ † 2 representation of the matrix η, parametrizing a matrix Θ in such a way that ΘΘ † 2 will fit into currently known limits on η.
Actually in our strategy we come back to the bottomup scenario, as our analysis starts from the present state of knowledge on U PMNS mixing data in a form of an interval matrix and we examine directly whether matrices within are physically meaningful (i.e. are contractions). An extension of this idea allows us to define complete region of physical mixing matrices as a convex hull of U PMNS matrices which ensures that any physical mixing matrix can be constructed as a convex combination of U PMNS matrices. Now, let us emphasize the relation of our approach to the polar decomposition. In our analysis we use singular values as an indicator whether a given matrix is a contraction. However, it is known that eigenvalues of the polar factor, which follows from the definition, are equal to singular values of an initial matrix. Thus from that perspective a polar decomposition can be treated as a compact version of singular value decomposition. Nevertheless, from the numerical analysis perspective singular value decomposition algorithms are more natural, since they arise from eigenvalues decomposition of matrices AA † and A † A. Thus in most cases, in order to obtain an algorithm for a polar decomposition we have to translate algorithms for the singular value decomposition.