Light Majorana neutrinos in (semi)invisible meson decays

We reconsider decays of pseudoscalar mesons (P) to neutrino pairs and possibly additional photons in presence of (light) Majorana neutrinos. For this purpose we derive a model-independent general parametrization of neutrino mass matrices with physically interpretable and irreducible set of parameters. The parametrization is valid for any number of neutrinos and interpolates smoothly between the heavy Majorana and the (pseudo)Dirac neutrino limits. We apply the new parametrization to the study of P→νν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P \rightarrow \nu \nu $$\end{document} and P→ννγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P \rightarrow \nu \nu \gamma $$\end{document} decays within the SM extended by additional singlet fermions. We update the SM predictions for the branching ratios of Bs,d→ννγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{s,d} \rightarrow \nu \nu \gamma $$\end{document} and discuss the sensitivity of the Bs,d→Emiss(γ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_{s,d} \rightarrow E_{\mathrm{miss}} (\gamma )$$\end{document} decays to neutrino mass and mixing parameters.


Introduction
The discovery of neutrino oscillations [1] implies the existence of at least two massive neutrino species. On the other hand, many theoretical models of neutrino mass generation, including the simplest canonical see-saw mechanism [2][3][4][5][6], predict the existence of additional electromagnetically neutral massive fermions. In general, the neutrino spectrum consists of 3 + n N fermions. Three of them are the so-far observed standard model (SM) like (ν M j ) neutrinos. Possible additional n N massive neutrinos (N M ) have not yet been observed. In the following and without loss of generality we assume them to be of Majorana type. 1 In the last few decades many different mechanisms have been proposed to explain the smallness of the observed neutrino masses. In a e-mail: blaz.bortolato@ijs.si b e-mail: jernej.kamenik@cern.ch (corresponding author) 1 The model of Dirac SM-like neutrinos is a special case with n N = 3 and with all Majorana mass terms set to zero. Its spectrum consists of 3 Dirac neutrino fields, which can be written as a superposition of ν M j and N M j fields.
the canonical see-saw mechanism for example, heavy N M k neutrinos induce small Majorana masses for the observed ν M j neutrinos via their Yukawa interactions. In this scenario N M k neutrinos are typically too heavy to be directly produced in terrestrial experiments [7]. On the other hand, in models with approximate lepton number conservation, these new degrees of freedom could also naturally appear at low energies, see e.g. [8]. In fact there are several circumstantial motivations for considering additional light N M k neutrinos. Massive neutral fermions which are long lived enough compared to the age of the universe and have mass in the range 2 keV m N k 5 keV [9,10] are good warm dark matter candidates [11][12][13][14]. Additional N M k neutrinos with masses in the range 1GeV ≤ m N k 100GeV are also predicted in models of cosmological Baryon asymmetry generation through neutrino oscillations [15,16]. Finally, persistent tensions in the interpretation of some neutrino oscillation experiments and cosmological observations might imply the existence of additional N M k neutrinos with masses at the eV scale, see e.g. [17].
An important aspect of neutrino mass model building involves consistently taking into account existing experimental information on low energy neutrino masses and mixings. In principle these inputs can be used to reduce the number of free model parameters. In practice however, this requires a detailed a priori knowledge of how elements of neutrino mass and mixing matrices are connected with each other. In the limit of heavy N M k neutrinos such connections are given explicitly by the Casas-Ibara parametrization [18]. In the last few years, more general parameterizations have been proposed, which do not rely on expansions in small mass ratios and are thus valid away from the limit of heavy N M k neutrinos. To date such parametrizations have been found for the case of two [19,20] or three [21][22][23] additional N M k neutrinos. Most recently, a master parametrization applicable for the most general case including neutrino mass generation beyond see-saw models has also been proposed [24]. However, the generality comes with several drawbacks: its param-eters lack intuitive physical interpretability, the connections between the N M k neutrino masses and the SM-like neutrino Yukawa couplings are somewhat obscured. In this paper we build upon and extend previous work [19,23] and derive a model-independent general parametrization of the neutrino mass matrices that covers and interpolates between all seesaw like scenarios, including the heavy Majorana mass limit and the pseudo-Dirac case for any number of additional N M k neutrinos (see e.g. Ref. [25] for an explicit model realization of such a scenario). Its main purpose is to better map out the possible neutrino mass parameter space and to help create a consistent picture of N M k neutrinos at low energies, which is a starting point for developing UV neutrino mass models. 2 We demonstrate the usefulness of the parametrization using the example of P → νν and P → ννγ decays, previously studied in Ref. [27] in the context of light dark matter searches, where P is a neutral pseudoscalar meson and ν includes both ν M j as well as possibly N M k if kinematically allowed. We estimate the contribution of these two decay topologies to the effective invisible decay widths of neutral mesons (assuming N M k are long lived enough to escape the detectors) and show how they are sensitive to neutrino mass and mixing parameters. Along the way we also update the theoretical predictions for B s,d → ννγ in the SM using state-of-the-art inputs [28] for the relevant hadronic parameters and the associated uncertainties. The P → νν decays are helicity suppressed and therefore negligible in the SM with Dirac neutrinos [27] as well as in the limit of heavy N M k . However in models with light N M k , such that they can appear in the final state, their branching fractions can become significant. On the other hand P → νννν decays are not helicity suppressed [29], however their contributions to the invisible P decay widths turn out (within our assumptions) to be completely negligible. Experimentally, the best sensitivity is expected from B s,d meson decays into unobserved decay products (registered as missing energy E miss in the detector) which have already been searched for by the Belle [30] and BaBar [31] collaborations. At present the tightest upper limit of Br (B d → E miss ) < 2.4 × 10 −5 at 90% confidence level is provided by BaBar [31]. While searches for invisible B s decays have not yet been attempted, they are planned at the Belle II experiment, which is also expected to improve significantly the upper bound on Br (B d → E miss ) [32].
The paper is organized as follows. In Sect. 2 we derive a general parametrization of neutrino mass matrices for an arbitrary number of additional massive fermionic singlets and explore the heavy Majorana neutrino limit and the pseudo-Dirac limit. We also present the basic properties of the parametrization and how these can be used to extract neu- 2 We note that in UV models where B − L is gauged (like in U (1) B−L [26] or left-right symmetric models [3,4]), anomaly cancellation requires exactly three right-handed neutrinos (n N = 3). trino parameters from experiments. In Sect. 3 we study the P → νν and P → ννγ decays separately, estimate their contributions to the invisible decay widths of B d,s mesons, and discuss their dependence on the neutrino parameters. We summarize our findings in Sect. 4. Analytical expressions for the P → νν and P → ννγ decays as well as the details on the perturbative and non-perturbative QCD inputs used in this work are given in Appendix A, while Appendix B contains the details on the derivation of a lower bound on the Frobenius norm of the neutrino mixing matrix.

Setup and notation
We consider a family of neutrino models at low energies which are described by the Lagrangian 3 L = L SM + L N , where L SM is the Standard Model (SM) Lagrangian and L N is given by: The first term is the Dirac mass term where the Yukawa coupling matrix y is implicit in where v is the Higgs VEV. The second term is the Majorana mass term of chiral right-handed neutrinos N bR . The mass matrices of models which preserve SM local symmetries and are renormalizable (chiral left handed neutrino mass terms are forbidden) form a symmetric block matrix M of the form 4 The matrix M can be diagonalized by the unitary matrix L in the following way Here where (ν j L ) m and (N k R ) m are the mass eigenstates which are connected to the gauge interaction eigenstates via Without loss of generality, all diagonal elements of D ν and D N can be taken as real (via suitable unphysical phase rotations of the neutrino fields). With this notation one can write down the interaction terms in the Lagrangian expressed by ν M j and N M k neutrino mass eigenstates and by l m ∈ {e, μ, τ } charged lepton mass eigenstates,

Derivation
In this section we consider a model with n ν neutrinos ν M j and n N neutrinos N M k . Matrices U , V , X and Y clearly depend on the neutrino masses, thus it is appropriate to have a simple parametrization of these matrices in terms of physical neutrino parameters. Parametrization of this type can be derived with a few simple steps. We start by decomposing the V matrix into V = gS : n N > n ν , g P : n N ≤ n ν , where g is a n ν ×n ν complex matrix, S is given by S n ν ×n N = [I n ν ×n ν ,Ŝ n ν ×(n N −n ν ) ] and P n ν ×n N is a projection matrix which in the case of a normal hierarchy of ν M j neutrino masses is given by TheŜ matrix is a general complex n ν × (n N − n ν ) matrix. Its elements can be chosen freely. By using the definition of the transition matrix, L M diag L T = M one obtains the relation which, depending on n N , can be rewritten as Equation (9a) is already in the desired form. With a few assumptions it can be written as R R T = I , where R is an invertible orthogonal complex squared matrix, which connects U D 1/2 ν with g(−S D N S T ) 1/2 . On the other hand Eq. (9b) is not yet of the desired form. In the special case n N = n ν the projection matrix P becomes the identity matrix, and we can decompose the left-hand side of Eq. (9b) into U D ν U T = U D 1/2 ν D 1/2 ν U T . Notice that on both sides all matrices have the same shape. In the case n N < n ν the diagonal mass matrix D ν is not invertible. This can be seen directly by looking at the block mass matrix M, Of the first n ν rows, maximally n N are independent. Similarly, of the last n N rows, also maximally n N are independent. Therefore one can construct at least n ν +n N −2n N = n ν −n N independent eigenvectors for the matrix M with zero eigenvalues. From here, there are maximally n N non zero masses m ν j in D ν . This property can be written as where P T D ν P is a n N × n N diagonal matrix of all non zero eigenvalues of D ν , therefore an invertible positive definite matrix. By specifying D ν one can use this property to find P. Using Eq. (11) one can finally rewrite the left-hand side of Eq. (9b) in terms of n N × n N matrices, Here we have multiplied the equation by a n ν × n N arbitrary matrixP from the right-hand side and byP T from the lefthand side. At this point we assume that all matrices, which are products of matrices inside brackets in Eqs. (12) and (9a), are invertible, D N matrix is invertible for all pairs (n ν , n N ) and that matrices U and g are also invertible in the case n ν < n N . Cases in which these assumptions do not hold can still be handled by the parametrization we are deriving, by taking appropriate limits. Within our assumptions, both Eqs. (9a) and (12) can be expressed in the form R T R = I or R R T = I equivalently, where R is a min(n ν × n ν , n N × n N ) general complex orthogonal matrix. It links U and V matrices through where the n ν × n N matrix Q is given by In case n ν > n N one obtains the relationP T (V − U Q) = 0 which leads to V = U Q, sinceP T is arbitrary up to the above assumptions. Eq. (13) together with the unitarity condition L L † = I finally leads to the desired parametrization of U and V matrices Here A is a n ν × n ν unitary matrix. The I + Q Q † hermitian matrix is positive definite, which means that it has precisely one positive definite square root which is also invertible.
Parametrizations of X and Y matrices are similarly obtained from the unitary condition L L † = I using the above derived results.

Main formulae
Below we give the full set of equations which define the parametrization of the neutrino mixing matrices: where A n ν ×n ν and B n N ×n N are unitary matrices and Q is defined by The M D = v/ √ 2y and M M matrices are therefore parametrized via where D ν and D N are diagonal neutrino mass matrices defined up to arbitrary phases (for each diagonal element). The derived parametrization expresses neutrino mixing matrices appearing in the Lagrangian in terms of unitary matrices A and B, complex orthogonal matrix R, neutrino masses contained in D ν and D N and in case n N > n ν also a general complex matrixŜ which is hidden inside S = [I n ν ×n ν ,Ŝ n ν ×(n N −n ν ) ]. The projector P in case n ν ≥ n N must be chosen such that P T D ν P is a diagonal matrix of all non zero eigenvalues of D ν and that the relation P P T D ν P P T = D ν holds. If n ν > n N , the choice of P depends on the hierarchy of ν M j neutrino masses.

Physical case n ν = 3
In the following we discuss the explicit form of our parametrization, it's limits and parameter counting, for the realistic case n ν = 3 and various possible choices of n N . In the case n ν = n N = 3 previously studied in Ref. [23] the R matrix can be parametrized with 3 complex angles: where c j = cos(φ j + iθ j ) and s j = sin(φ j + iθ j ) with φ j ∈ [0, 2π) and θ j ∈ IR for j ∈ {1, 2, 3}. Free signs in each of the three matrices must be equal (both + or both − in each matrix). The projector P in this case reduces to the identity matrix.
In the case n ν = 3 with n N = 2 previously studied in Ref. [19] one of the ν M j neutrinos is massless, thus the D ν matrix can be parametrized by in case of inverted hierarchy (IH). In both scenarios the R matrix is described by one complex angle: where φ ∈ [0, 2π) and θ ∈ IR. As before the free signs in the R matrix must be equal (both + or both −). The projector P for NH (IH) is given by: Also, in the case n ν > n N one can rename the P R matrix to R n ν ×n N , since the projector P is always multiplied by R from the right-hand side.
On the other hand, matrix B is a general n N × n N unitary matrix for all pairs (n ν , n N ), while matrix A is a general n ν × n ν unitary matrix only if n ν ≤ n N , since if n ν > n N not all generators of U(n ν ) are required to form A so that the (M D ) n ν ×n N and (M M ) n N ×n N matrices are fully parametrized with all n 2 N + n N (2n ν + 1) parameters. The number of real parameters for each matrix which appears in the derived parametrization for the case n ν = 3 is given in Table 1.

Heavy N k limit
In the limit where ||D ν ||/||D N || 1 and ||Q|| 1 the parametrization simplifies as and In case n N ≤ n ν we immediately recognize the original Casas-Ibara parametrization [18]. By combining the expressions of M D and M M matrices one can also construct the well known effective ν M j neutrino mass matrix This formula is valid for all pairs (n ν , n N ). In this scenario the PMNS matrix which is given by (5)) is approximately unitary. To clarify, the PMNS matrix maps fields of observed neutrinos from the mass basis into the flavor (gauge) basis. In the heavy neutrino limit low energy processes can only involve ν M j neutrinos.

Dirac neutrino limit
In the scenario with n ≡ n ν = n N and M M = 0 neutrinos can be described by Dirac fields (linear combinations of Majorana fields which are not Majorana fields). The condi- The phases in front of neutrino masses are arbitrary and do not affect any measurable quantities. Therefore we can fix the phases by choosing Signs are arbitrary and independent of each other. The transition matrix L then takes the following form At this point one can diagonalize M D through a biunitary transformation Therefore the mass term in the Lagrangian can be written as where are Dirac neutrino fields (not Majorana fields). One can use this equation to define ν M j and N M j Majorana fields starting from the Dirac field. Such definition of Majorana fields can then be applied outside the Dirac limit. In this way, the obtained model is fully consistent with both the general Majorana neutrino model (outside the Dirac limit) as well as with the Dirac neutrino model (in the Dirac limit). We use this model in Sect. 3. Finally, the Lagrangian (for the case n = 3) can be written as From the last equation one can recognize the (unitary) PMNS matrix as Its matrix elements are precisely the same as in the heavy neutrino limit if we keep A matrix unchanged. However, now the observed neutrinos are 3 Dirac fermions which are specific linear combinations of ν M j and N M j fields. Note that the measured values of the PMNS matrix can be directly related to the underlying neutrino parameters only in the heavy neutrino and (pseudo-)Dirac limits. Outside these two limits the precise relations are non-trivial, since the PMNS matrix is no longer unitary. However, results from experiments which measure the PMNS matrix elements can still be used to constrain parameters of the general low energy neutrino models.

Scaling relations
Matrices U † U , U † V and V † V satisfy the relation: where ||U † V || 2 = j j |(U † V ) j j | 2 is the Frobenius norm and n ν is the number of ν M j neutrinos. The property can be derived by a few simple steps. First define matrix U as Its relevant properties are U † = U and U 2 = U. The last one holds due to the unitarity condition UU † + V V † = I n ν ×n ν . From here, one gets n ν +n N E=1 U C E U E D = U C D . Therefore where the identity Tr{V † V } = Tr{V V † } was used. By the definition of the Frobenius norm, relations like |(U † V ) jk | < ||U † V || hold. These can be used to constrain matrix elements especially in scenarios with nearly degenerated N M k neutrino masses. In the heavy Majorana mass limit m ν /m N 1 Frobenius norms can be used as direct measures of nonunitarity in the 3 × 3 light neutrino sector, complementary to the PMNS matrix (U † O L ), since in this limit exact 3 × 3 uni- In general ||UU † ||, ||U V † || and ||V V † || norms can be expressed as a function of eigenvalues μ j of the Q Q † matrix: where the first identity can be related to the determinant of the PMNS matrix defined in the m ν /m N 1 limit since det(U † U ) = n ν j=1 (1 + μ j ) −1 . Number of ν M j neutrinos is also number of required quantities to completely describe norms.
In a special case in which n ν = n N = 3 and neutrino masses are degenerated, D ν = m ν I and D N = m N I , equations above can be simplified, where λ is the largest eigenvalue of the R † R matrix. Properties used to derive these equations are described in Appendix B. In this special case norms depend only on λ ∈ [1, ∞) and m ν /m N ratio. However if m ν /m N 1 equations above can be reduced to: where ξ = √ m ν /m N √ λ. Expressions in Eq. (37) can be regarded as scale relations for norms with respect to a dimensionless scale parameter ξ . Again the first identity is closely related to det(U † U ) ≈ 1/(1 + ξ 2 ) measuring PMNS matrix unitarity. If R matrix depends only on one θ angle, we find λ = exp(2|θ |). In other cases λ depends on all parameters of R matrix, however if R matrix is parametrized in one of specific ways, λ depends only on θ k angles and not on φ k angles. More details in Appendix B. Dependence of norms with respect to θ k and φ k parameters are shown in Figs. 1 and 2.

Extracting neutrino parameters
Elements of the O † L U matrix contain information about N M k and ν M j neutrino masses, R matrix and S matrix elements. In the case n ν = n N with fixed D ν one can simply obtain R and D N by properly diagonalizing the left-hand side of where In case n ν > n N projectors P on both sides should be added, but if n N > n ν the diagonalization becomes meaningless, since the S matrix appears in various places on the right-hand side of the equation.
Processes involving neutrinos, but not charged leptons in the initial and final state usually do not depend on matrices O † L U and O † L V , since masses of charged leptons (in loops) are small compared to the masses of W and Z bosons. Therefore observables in these processes depend only on neutrino masses, R matrix elements and if n N > n ν also onS matrix elements. In the heavy N k and pseudo-Dirac neutrino limits such processes are sensitive to (small) non-unitary corrections to the PMNS matrix, but not to PMNS matrix elements themselves. As explained above, outside of both limits the precise relation between the experimentally measured PMNS matrix elements and neutrino parameters is highly nontrivial and thus more difficult to interpret.

Neutral meson decays to two neutrinos and (un)resolved photons
In this section we consider the decays of pseudoscalar mesons (P) to neutrinos and possibly additional photons (γ ), as prospective venues to constrain low energy neutrino parameters. We consider both signatures of P → E miss as well as P → E miss γ , where E miss is an energy imbalance registered in the particle detector. Within our theoretical setup and assuming a 4π detector coverage with finite EM energy resolution, the decay products which contribute to P → E miss include stable enough neutrinos as well as unresolved (soft) photons where the dots denote additional multibody decay channels which are further suppressed. The energy of photons present in the final state (γ * ), should be less than the energy resolution E 0 of the EM detector. According to Refs. [36] and [37] the EM calorimeter of Belle II has the energy resolution in the range of 20-50 MeV. For concreteness, we use the value E 0 = 50 MeV throughout this paper. Light neutrinos are clearly invisible to the detector, while heavy neutrinos may decay into lighter and observable decay products before they escape. Thus in general non-trivial conditions which depend on neutrino parameters are imposed on the branching ratio Br (P → E miss (γ )). We aim to study this dependence and in particular to estimate theoretical upper bounds on Br (P → E miss (γ )) based on the relevant experimental constraints. To do so, we first briefly discuss the basic properties of P → νν and P → ννγ decays in the following framework. 5 We use the Majorana neutrino model with n ν = n N = 3. Details about the model and all relevant equations to reproduce the results are described in Appendix A. We use the following compact notation along with the U matrix to estimate the relevant branching ratios. 6 Numerical results and plots in this section are cal- Fig. 1 Examples of the ||U † U ||, ||U † V || and ||V † V || norms as functions of the θ 2 and θ 3 parameters in a model with n ν = n N = 3. The other model parameters are fixed as θ 1 = 0 and φ k = 0 for all k (implying θ k → −θ k for each k). Light neutrino masses are set to satisfy present experimental constraints with normal mass hierarchy [35] and m ν2 /m ν1 = 2. The N M k neutrino masses are set to m N = 1 GeV (upper plots) and m N1 = 1 MeV, m N2 = 1 GeV and m N3 = 100 GeV (lower plots) Fig. 2 The ||U † U || norm as a function of θ 3 and φ 2 at φ 1 = 0 (left), as well as φ 2 and φ 1 at θ 3 = 12 (right), in a n ν = n N = 3 scenario with m N1 = 1 MeV, m N2 = 1 GeV and m N3 = 100 GeV. Other θ j and φ j parameters are set to zero. Light neutrino masses are set to satisfy present experimental constraints with normal mass hierarchy [35] and m ν2 /m ν1 = 2 culated using expressions and numerical inputs described in Appendix A. For light neutrino masses m ν j we impose experimental constraints from [35]. For concreteness, in cases where observables significantly depend on light neutrino masses, we assume a normal mass hierarchy with m ν 2 /m ν 1 = 2.

P → νν
The P → ν C ν D decay is helicity suppressed and therefore highly sensitive to neutrino masses as can be seen from Eq. (A9). Assume for the moment that U and R matrices are purely real or have a negligibly small imaginary part. Then the dependence of the decay width on U matrix elements can be factorized such thatBr (P → ν C ν D ) ≡ Br (P → and/or Re(U 2 C D ). Both are unaffected by the change U C D → U * C D . In the case n N ≤ n ν , this implies a symmetry θ → −θ, since R(φ 1 , φ 2 , φ 3 , θ 1 , θ 2 , θ 3 ) * = R(φ 1 , φ 2 , φ 3 , −θ 1 , −θ 2 , −θ 3 ). Exchanging θ k with −θ k for a single k is in general not a symmetry of observables. ν C ν D )/|U C D | 2 becomes independent of |U C D | and its dependence on neutrino masses m ν C and m ν D is shown in Fig. 3.
Note the different kinematical behavior of P → N M k N M k and P → ν M j N M k due to the the sign difference sgn[Re , see Eq. (A9). Furthermore, assuming that neutrinos are (nearly) degenerate D N ≈ m N I , D ν ≈ m ν I with m N m ν , we can decompose the branching fraction of the P → νν decay into In this scenario, norms satisfy scaling relations described in Sect. 2.7. From here it is easy to see that within these assumptions the branching ratio Br (P → νν) is largest if m N is of a non-negligible fraction of the P mass and the ξ scale is sufficiently large. This can be seen from Fig. 4 for the case of B s decays. Due to scaling relations similar behavior like in Fig. 4 is expected in case θ 3 is exchanged by k θ k . In general the branching fractions Br (P → νν) are minimal in the Dirac limit. Precise values depend on the light neutrino masses, but are in any case negligibly small compared to experimental resolution of any currently foreseen experiments [27]. On the other hand, the maximal values of The maximal values of Br (P → νν), however, also crucially depend on the experimentally allowed values of U matrix elements. It turns out that currently the constraints are mildest for m N k in the range of a few GeV [38], in particular, there |(U † V ) jk | 2 < 5 × 10 −5 as reported by the DELPHI collaboration [39]. This bound implies ξ 1 and in turn ||V † V || ||U † V ||. If we thus assume approximately degenerate N k with masses m N k 0.6M B ∼ 3 GeV, and all (U † V ) jk matrix elements saturating the current experimental bound in that mass region we obtain

P → ννγ
In the limit m C , m D m P , the branching ratio Br (P → ν C ν D γ ) is approximately independent of neutrino masses From here the identity Br (P → ν C ν D γ ) |U C D | 2B r (P → ν C ν D γ ), whereB does not depend on U C D , holds. Therefore in the case ||D N || m P the branching ratio of the P → ννγ decay is given by which coincides with the SM prediction Br (P → ννγ ) SM (with three massless neutrinos). In the last step the identity in Eq. (32) was used. TheBr (P → ν C ν D γ ) function takes lower values if neutrino masses m ν C and m ν D increase as can be seen from the left-hand side plot in Fig. 5 for the case of the B s → N k N k γ decay where we plot the photon energy spectrum of the decay (normalized to the B s lifetime and the relevant |U| 2 matrix element) as a function of the neutrino mass.
Using properties of the ||UU † || norm described in Appedix B one finds that in cases n ν = 3 and n N ∈ {2, 3} with max(D ν )/ min(D N ) 1 and max(D ν ) m P , the branching ratio of P → ννγ decay lies in the interval: The maximal value of the branching ratio takes place in the SM. By taking ||D ν || = ||D N || → 0, form factors F A (q 2 ) and F V (q 2 ) from the most recent estimate [28] and integrating Eq. (A11) over the whole phase-space, we find the SM predictions for branching ratios Br (B s,d → ννγ ) to be where the O(30%) uncertainties are dominated by the relevant hadronic form factor estimates, see also Appendix A. As can be seen from Table 2, the most recent form factor inputs lead to somewhat reduced predictions compared to previous estimates.
In addition to the branching ratio, measuring the photon energy spectrum in P → ννγ decays would in principle  allow to infer on the mass spectrum of the neutrinos appearing in the final state. In particular, the average photon energy E γ , defined as is inversely correlated with final state neutrino masses, as can be seen from the left-hand side plot in Fig. 5. Assuming neutrinos are completely unobserved, the P → ννγ decay contributes also to the invisible P decay width effectively due to the finite resolution of any electromagnetic calorimeter, since photons with energies lower than some threshold energy E 0 of the detector are not registered. This Fig. 6 The branching fraction of the B s → ννγ * decay with a soft photon (E γ < E 0 ) as a function of θ 3 and (m) in scenario with degenerate N M k neutrino masses with R = R(θ 3 ) (right plot), and in scenario with m N1 = m N2 = 2 GeV, m = m N3 and R = R(θ 3 ) (left plot). In both scenarios parameters θ 1 , θ 2 , φ k for all k are zero contribution is simply given by where the other integrals are performed over the whole available phase space. On the right-hand side plot in Fig. 5 we show the threshold energy and neutrino mass dependence of the P → ννγ * decay branching fraction for the case B s → N M k N M k γ * . Specifically, in the SM and for a threshold energy of E 0 = 50 MeV we obtain the pre- We show in Fig. 6 a typical dependence of the B s → ννγ * branching fraction on the neutrino mixing parameters in two scenarios for N M k neutrino masses, E 0 = 50 MeV and with R = R(θ 3 ), θ 1 = θ 2 = 0. In comparison to Fig. 4 we observe that these contributions to B s → E miss are generically still more than three orders of magnitude smaller than the current upper limit on B s → νν and thus completely subleading. However, at the same time, P → ννγ * are always much bigger than P → νννν [29] and are thus expected to dominate P → E miss in the (pseudo)Dirac neutrino limit. In addition, P → ννγ might contribute effectively to P → E miss also due to other detector effects, such as non-perfect 4π coverage. Such contributions are however difficult to model without a detailed knowledge of the detector components and geometry and we leave such a study to the experimental collaborations.

Conclusions
In this work we have reconsidered Majorana neutrino mass models and derived a model-independent general parametrization of neutrino mass matrices with physically interpretable and irreducible set of parameters. The parametrization is valid for any number of left-handed (as in SM) and right-handed (gauge singlet Majorana) neutrinos and for all mass hierarchies. In particular, in the heavy Majorana neutrino limit we recover the standard Casas-Ibara parametrization [18], while the parametrization nicely interpolates also through the (pseudo)Dirac neutrino limit.
We have applied the new parametrization to the study of P → νν and P → ννγ decays within the SM extended by n N = 3 additional singlet neutrinos. 7 Along the way we have updated the SM predictions for the branching ratios of B s,d → ννγ decays and found almost an order of magnitude smaller values compared to previous estimates, mainly due to a recent reevaluation of the relevant hadronic form factors.
Finally, we have discussed the sensitivity of the B s,d → E miss (γ ) decays to neutrino mass and mixing parameters. In the case of B s,d → E miss , for typical EM calorimeter threshold energies and assuming 4π coverage, the dominant contribution could still come from B s,d → νν decays where one of the final state neutrinos is predominantly a SM gauge singlet of mass of the order a few GeV. However, the maximum allowed branching ratios, given by the current experimental bounds on the relevant neutrino mixing matrices, are at least four orders of magnitude below the direct limits from the B factories. It remains to be seen if Belle II can reach the required sensitivity to constrain the parameter space of neutrino mass models in this interesting region.
In the case of B s,d → E miss γ decays, additional light neutrinos in the final state could affect both the branching ratios as well as the photon energy spectra and thus in principle allow to extract information on the neutrino mass parameters. Unfortunately, however, possible deviations from SM predictions (i.e. the limit of massless neutrinos) are theoretically constrained and at most comparable to current uncer- 7 The parametrization can trivially be applied to other neutral current mediated processes involving pairs of neutrinos. In the case of charged current mediated processes commonly used to constrain the PMNS matrix (see e.g. Ref. [43]), or in searches for massive Majorana neutrinos with lepton number violating signatures [33,44], one instead needs to consider additional dependence on the unitary O L matrix. tainties due to the limited knowledge of the relevant form factors. Any relevant experimental sensitivity to neutrino mass parameters is thus conditional upon an improved understanding of the relevant hadronic parameters, which could possibly come from future Lattice QCD studies (see Refs. [45,46] for current prospects).

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The results presented in the paper are theoretical and no data beyond that shown in the figures (readily derivable using analytic formulae and inputs presented in the paper) has been generated.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Calculation of P → νν and P → ννγ decay widths
Below we summarize our calculation of the B s,d → νν(γ ) decay widths as discussed in the main text. We consider the physically relevant scenario with n ν = 3 neutrinos ν M j and n N neutrinos N M k . The Majorana fields are defined in a way such that the Dirac limit can be approached analytically with both D ν and D N matrices being positive semi-definite. For other choices of the relecant phase factors one should properly redefine the mixing matrix U. Note that our calculation can be applied also to the corresponding K , D meson decays, with suitable quark flavor replacements.
We calculate the decay widths B s,d → νν(γ ) using the relevant effective weak Lagrangian [47,48] where the leptonic current J μ C D is given by The relevant loop function X (x t ) can be written as where X 0 (x t ) is the Inami-Lim function [47] X 0 (x t ) = x t 8 and the leading QCD corrections are parametrized by X 1 whose explicit expression can be found in Refs. [48,49]. Here x μ = μ 2 /M 2 W , x t = m 2 t /M 2 W and the MS QCD renormalization scheme is assumed throughout. In the following we compress the common constant prefactors entering the Lagrangian into Above and in the following we have suppressed the light flavor (q = s, d) indices where the identification of the relevant B (q) meson flavor is unambiguous. For the P → νν decay we parametrize the relevant hadronic matrix elements in the standard way 0|bγ μ q|P( p) = 0, (A6a) where f P is the relevant P = B s,d meson decay constant. In particular, we use f B s = 224 MeV and f B d = 186 MeV from Ref. [50].
In the case of the radiative decay, only the emission of photons from the hadronic part is relevant and is parametrized by the relevant radiative form factors Here q = p C + p D = p −k and q 2 = m 2 P −2m P E γ . For the axial (A) and vectorial (V) form factors F A (q 2 ) and F V (q 2 ) we take the most recent estimate [28] parametrized by where F(0), σ 1 , σ 2 and M R parameters for P = B s,d are given in Ref. [28]. After a quick calcuation one finds the expression for the P → ν C ν D decay width where (A10) In the Dirac limit the process P → ν M j N M j is forbidden. Similarly, the triply differential P → ν C ν D γ decay width (for Majorana neutrinos) is given by Note that due to Majorana nature of neutrinos the full integral over the solid angle d gives 2π instead of the usual 4π . From kinematic constraints one furthermore obtains where is equal to while k · p C is given by Finally, in our numerical results we use μ = m Z = 91.2 GeV, α em = 1/137, m B s = 5.37 GeV, m B d = 5.28GeV, m W = 80.4 GeV, |V tb V * ts | = 0.0403, |V tb V * td | = 0.00875, sin 2 (θ W ) = 0.22, τ B s = 1.51 ps, τ B d = 1.52 ps, α s (m Z ) = 0.118 and m t (m Z ) = 172GeV [51].

Appendix B: Lower bounds on ||UU † ||
In this section we formally prove Eq. (45) and discuss additional properties of the ||UU † || norm. We assume a model with n ν light neutrinos (m ν j /m P 1) and n N heavy neutrinos (m ν j /m N k 1). From properties of P → ννγ decay follows that inequality: holds. Moreover from upper limit of the branching ratio Br (P → ννγ ) we find: This is a direct consequence of Eq. (32). From here, within assumed model we have: Value of the norm is due to derived parametrization of neutrino matrices directly related to the eigenvalues μ j of the Q Q † matrix through relation: This equation is obtained by diagonalizing Q Q † matrix inside Frobenius norm. For U V † and V V † matrices similar formulas can be found: We use Eq. (B4) as a starting point to determine theoretical lower bounds on the ||UU † || norm for different choices of n ν and n N .

Case n ν ≥ n N
In scenario with n ν ≥ n N we use following theorems. Notation 1 Let X be a n×n hermitian matrix, then we denote its eigenvalues as λ j (X ), where λ 1 (X ) < · · · < λ n (X ). We denote by a the largest eigenvalue of D ν and by b the smallest eigenvalue of D N , therefore D −1 N ≤ 1/bI . Using Theorem 1 we get:

Definition 1 (Loewner order) Let
By Theorem 3 we find λ j (Q Q † ) = 0 for j = 1, . . . , n ν − n N . Next for j > n ν − n N , we apply Theorems 2, 3 and 1 in this order to obtain: From here, we finally get: From the property (R † R) −1 = (R † R) * follows directly, that if λ(R R † ) is an eigenvalue of the R † R matrix, then 1/λ(R † R) is also eigenvalue of R † R matrix. A consequence of this property is that for any (2n + 1) × (2n + 1) orthogonal matrix R, the R R † matrix has at least one eigenvalue equal to 1. Therefore, if R is a general n × n complex orthogonal matrix, then R † R matrix have maximally n/2 (integer part of n/2) eigenvalues which are greater than 1. This implies: ||U † U || 2 ≥ n ν − n N /2 . (B10) Using property: a better lower bound can be obtained: Eigenvalues of R R † matrix depend only on θ k parameters if R matrix is parametrized in the following way: where R 1 , . . . , R n are real orthogonal matrices and therefore unitary. We can absorb them in the process of diagonalizing R R † matrix: From here, eigenvalues of R R † matrix must depend only on θ k parameters.

Case n ν < n N
In scenario with n ν < n N is more difficult to obtain eigenvalues of Q Q † matrix, since S matrix is present in it. In case n N ≥ 2n ν the lowest bound on ||UU † || norm is 0. This can be proven using R = I n ν ×n ν , D ν = a I n ν ×n ν , S = [I n ν ×n ν , ix I n ν ×n ν , 0 (n N −2n ν )×(n N −2n ν ) ], where x ∈ R and: where C is a positive semi-definite diagonal matrix. From here one gets: For x = ± √ b/c matrix S D N S T is invertible. When x approaches to √ b/c eigenvalues of Q Q † matrix become very large and therefore the lowest value for the ||UU † || norm is 0.