Dirac-Majorana neutrinos distinction in four-body decays

A novel method to differentiate the effects of Dirac and Majorana (D-M) neutrinos in four-body decays has been discussed in arXiv:2106.11785. There, it is concluded that the back-to-back kinematic scenario seems to avoid the constraint imposed by the"practical Dirac-Majorana confusion theorem", as one does not need to fully integrate over neutrino and antineutrino momenta. In this paper, we propose to analyse radiative leptonic lepton-decays ($\ell\to\ell'\nu\bar{\nu}\gamma$), as an independent alternative process to study the possible Majorana nature of neutrinos. Our approach demonstrates that, in the back-to-back kinematic configuration (for the $\ell'- \gamma$ and $\nu-\bar{\nu}$ systems, respectively), the distinction between Dirac and Majorana cases disappears when the inaccessible neutrino angle is integrated out, which might appear unexpected considering the claims in arXiv:2106.11785. Our results apply in absence of non-standard interactions, which can enhance generally the sensitivity to the neutrino nature.


Introduction
Nowadays the Dirac or Majorana nature of neutrinos, as well as the mechanism from which they acquire mass, is one of the unsolved puzzles whose solution seems to lie beyond the Standard Model (SM).There are many minimal extensions of it in order to account for nonzero masses and mixings for the active neutrinos; namely adding new gauge singlet fields, such as right-handed neutrinos.Some of them are the well-known νSM [2,3] and the seesaw-mechanisms [4][5][6][7][8].
If neutrinos were Majorana particles, lepton number would not be conserved and neutrinos would be their own antiparticles, i.e, no conserved quantum number would distinguish neutrino and anti-neutrino.There are some important proposals to probe the Majorana nature of the neutrino, such as the neutrino-less double beta decay (0νββ) of nuclei * [9][10][11][12][13], coherent scattering of neutrino on nucleus with bremsstrahlung radiation [14], and many other lepton number violating processes , most of them motivated by the propagation of a massive heavy neutral lepton.
However, the difference among Dirac and Majorana neutrinos, when the unobservable neutrinos momenta get fully integrated out, is proportional to some power of neutrino masses, assuming they enter left-handed charged weak currents, as in the SM.This is precisely stated by the "practical Dirac-Majorana confusion theorem" (DMCT) [49], and makes extremely challenging to distinguish experimentally between both possibilities, specially for tiny neutrino masses.
It is crucial that the difference between Dirac and Majorana nature can be settled independently of the mass of the neutrinos provided their momenta are not integrated out.Since neutrinos momenta are not experimentally accessible, we need a method to infer them, so that we can discern the neutrinos variables without the need of any explicit observation of them, which is extremely difficult any way at present.Kim, Murthy and Sahoo claimed [1] that we can deduce the neutrinos momenta working in the back-to-back (b2b from now on) kinematic configuration of a four-body decay with two final-state neutrinos, where the difference between Dirac and Majorana cases does survive irrespectively of the neutrino mass values, as long as neutrinos are not strictly massless.
This could be a really exciting and important strategy to distinguish the specific neutrino nature.Motivated by this method ("KMS method" from now on), we analyse radiative leptonic lepton-decays (ℓ − → ℓ ′− ν ℓ νℓ ′ γ) as an independent approach in order to distinguish the Dirac or Majorana nature of neutrinos.We emphasize that the quoted KMS method of analysis considers processes with ν ν final states with same-flavour neutrinos, in order to apply the quantum statistical properties of Majorana neutrinos.Nevertheless, the Majorana nature of neutrinos could be implemented even with ν ν final states with different flavours, working in the mass basis, where the quantum statistical properties could be applied to distinguish between Dirac and Majorana nature, as we shall discuss, see e.g.[50, 51, 55, 56]  * .This fact allows us to study many other processes as an alternative strategy (without hadronic transitions) with larger branching ratios (BR), such as the radiative leptonic decays of leptons, in which we focus here, that we will study neglecting possible non-standard interactions.
In Ref. [52], it was suggested that measurements of the Michel parameters in muon decay cannot help to distinguish the nature of neutrinos because they are of different flavours and the required anti-symmetrization of the amplitude for Majorana neutrinos does not apply.However, as we shall discuss, the reason why we cannot distinguish between Dirac and Majorana nature of neutrinos in this process was carefully explained by London and Langacker in [53]; they have shown that it is not possible, even in principle, to test lepton-number conservation in muon decay if the final neutrinos are massless and are not observed (which agrees with the DMCT theorem for the SM neutrinos interaction).Therefore, if we can infer some information on neutrinos, it becomes possible to get a non-vanishing difference between Dirac and Majorana cases, even for this type of leptonic decays.
As it was first introduced in [51], and discussed in several following works [50,[54][55][56], the Majorana nature could affect processes with ν ν final states involving different flavours.We emphasize that this effect is not a direct consequence of the presence of indistinguishable fermions in the final state but a result of the intrinsic Majorana fermion properties.Using the fact that flavoured neutrinos do not have a definite mass; i.e., the existence of two different basis, due to finite neutrino masses, this could lead to "internal interference effects" [53] when considering Majorana neutrinos.These effects can be analyzed with the usual implementation of Feynman diagrams and their corresponding Feynman rules for Majorana fermions.These interference effects can be explained in detail from first principles by studying the action of Majorana fields over the final asymptotic neutrino states in the QFT scheme, where now the operator Ψ can create either of the two states, and so can Ψ.
The main idea is now clear: we will study the possible effects of Majorana neutrinos in the radiative leptonic decay of µ and τ leptons assuming we can 'measure' them (i.e.infer their kinematic variables), avoiding integration over their momenta.In such cases, we shall prove that the difference between Dirac and Majorana nature of neutrinos is still present at the level of differential decay rate and explicitly depends on the neutrinos kinematics.Then, given the b2b configuration analysis as explained in Section 3.2, we will show that the D-M γ(p 5 ) Figure 1: Lowest order Feynman diagrams for radiative leptonic lepton-decay in the Dirac neutrino case.
difference vanishes upon integration on the neutrinos variables and discuss the origin of the discrepancy with the KMS result.This article is structured as follows: in section 2 we present the explicit radiative lepton-decay rate in the b2b configuration, taking into account Dirac (section 2.1) and Majorana (section 2.2) neutrinos.In section 3.1 we present the final energy and angular spectrum obtained from applying our approach.The main results are discussed: specifically showing that the differences for distinct neutrino nature do not survive independently of the non-vanishing neutrino mass, once the unobservable neutrino angle is integrated out.Then, in section 3.2 we track the origin of the absence of a difference between Dirac and Majorana cases in the radiative leptonic lepton decay, we elaborate on and clarify it, with the implementation of consistency tests and helpful discussions.Finally our conclusions are given in section 4. Very useful complementary computations of the phase space treatment, the b2b configuration and the specific branching ratio calculation can be found in Appendix A, B and C, respectively, where all the comments raised on Ref. [57] are corrected and clarified in detail.

Radiative leptonic ℓ-decay
Since m ℓ << M W , we can safely integrate out the W boson and use the Fermitype theory of weak interactions to describe the charged lepton decays.Then, the leading Feynman diagrams contributing to the radiative leptonic ℓ-decay are shown in figure 1.
From now on, we will be working in the basis where the mass matrix of charged leptons is already diagonalized and the flavor-eigenstate neutrino (ν L ) is taken to be the superposition of the mass-eigenstate neutrinos (N j ) with mass m j , that is, where j = {1, 2, 3, .., n} is tagging mass-eigenstate neutrinos.
In the mass basis, the ℓ − → ℓ ′− ν ℓ νℓ ′ γ decay consists of the subsets of all the n 2 separate decays of neutrino mass eigenstates allowed by phase space, i.e., the incoherent sum of the separate modes ℓ − → ℓ ′− N j N k γ [51]: Note that N represents an antineutrino for the Dirac neutrino case, but should be identified with N for the Majorana neutrino case (N =N c =CN T ).
Here it is important to remark that all the following analysis and results are valid only for massive neutrinos, no matter how light their masses could be, as long as they are nonzero.In this case, we have two different well defined basis (mass and flavor) and the corresponding neutrino nature will affect in a different way the differential decay rate computation as we shall see.In the case of massless neutrinos, only one basis exists and then eq.(2.1) does not hold, also the U (1) B−L symmetry of the Lagrangian is recovered and our results do not apply.For a helpful discussion see ergo [58].

Dirac Neutrino Case
In the Dirac case, the corresponding amplitude for the process ℓ − → ℓ ′− N j N k γ is given by: * where (see figure 1 for momenta convention) (2.4) Neglecting all final lepton masses, as a good approximation, the unpolarized * Where, when computing the decay rate, we must sum incoherently over the probabilities of all the allowed {j, k} channels.squared amplitude is: Finally, motivated by the KMS method, we need to work in the initial chargedlepton rest frame, specifically in the b2b kinematic configuration shown in figure 2, where the scalar products, neglecting the final lepton masses, are given by: * (2.6) Thus, in this kinematical configuration, the process is described in terms of only two variables, the final charged-lepton energy * , E 4 ≡ E ℓ ′ , and the angle between the neutrino and final charged-lepton directions, θ 24 = θ.With these considerations, the final Dirac amplitude is given by: The subindex '↔' denotes the b2b configuration.The 1/E ℓ ′ -dependence reflects the infrarred behaviour of the amplitude in the soft-photon limit.

Majorana Neutrino Case
Unlike the Dirac case, the properties of Majorana neutrinos have strong consequences in the amplitude.For Majorana neutrinos the decay modes ℓ − → ℓ ′− N j N k γ and ℓ − → ℓ ′− N k N j γ yield the same final states for k ̸ = j as well as for k = j (since N i = N i ), and hence their amplitudes must be added coherently.This is a result that can be obtained rigorously in the QFT formalism.
The possible first order Feynman diagrams for the ℓ − → ℓ ′− N j N k γ decay are shown in figure 3, leading to the total amplitude: where the relative sign between M(p 2 , p 3 ) and M(p 3 , p 2 ) stems from the application of Wick's theorem in presence of Majorana fermions (see, e.g.Ref. [60]).
It can be shown that, after summing over polarizations, Re (M(p 2 , p 3 )M * (p 3 , p 2 )) ∝ m 2 ν ≈ 0 due to the smallness of neutrino masses † .Thus (2.9) The computation, neglecting the final lepton masses, leads to the following squared amplitude for the b2b kinematic configuration: (2.10) Thus, the Majorana nature of neutrinos would generate a different behavior of the amplitude compared with the Dirac neutrinos case, eq.(2.7).

Energy and Angular Distributions
In this section we develop our own derivation, motivated by the KMS method, and compute the corresponding energy and angular distributions of the ν − ν and ℓ ′ − γ systems in the b2b configuration, respectively.We will finally discuss our results.

Our Approach
When we restrict ourselves to the b2b scenario, which is just a specific kinematic configuration of the general case, as we did before, we will adopt the following notation of the corresponding decay rate according to eq. (A.6) (see Appendix A for all the details): where ϵ = 1(2) for Dirac (Majorana) neutrinos, showing that it is just the standard differential decay rate evaluated in the specific b2b kinematics * , and the amplitudes involved were already quoted in the last section.
Using the |M D/M ↔ | 2 previously computed, we obtain the following neutrinos differential decay rate: where we already used the unitarity relations for the mixing matrix elements considering just the three active light neutrinos mass eigenstates.We can also consider the presence of a massive sector leaving the sum explicit, where the new sterile neutrinos could be produced if they are energetically allowed, or else could affect the unitarity relation if they are forbidden by kinematics.Nevertheless, as discussed in Ref. [56], the possible effect would be suppressed by the specific heavy mixing and can not be higher than 10 −4 , thus it should only be considered if that precision is needed.
Then the difference between Dirac and Majorana cases in the b2b scenario is: The difference is evident, unfortunately, the angle θ is not experimentally accessible, so we should rather integrate over it to get the final charged-lepton energy distributions for Dirac and Majorana cases.
In order to integrate over this inaccessible angle, we just need to rewrite the angle θ that appears in the right-hand side of the above equations in terms of our phase space variables θ ℓ ′ and ϕ, i.e., θ = θ(θ ℓ ′ , ϕ).
Since θ is the angle between ⃗ p ν and ⃗ p ℓ ′ , it is easy to get, as shown in eq.(B.9) that: Now, since this neutrinos differential decay rate does not depend explicitly on θ ν , we have the freedom to fix it to ease further computations * .Then, for example, we can choose the system in such a way that θ ν = 0 and thus we get cos θ = − cos θ ℓ ′ , which does not have a dependence on the ϕ angle, showing explicitly the azimuthal symmetry of this specific selection.Using this dependence, where 0 ≤ θ ℓ ′ ≤ π and 0 ≤ ϕ ≤ 2π so that all possible angular configurations between ℓ ′ and ν are accounted for, we have from eq.(3.4): Doing this, it is straightforward that the difference between Dirac and Majorana cases vanishes once the angular integration is made.We emphasize that this difference will vanish in any other selected frame, doing the angular integration properly, as derived in eq.(B.9).Specifically, for subsequent discussions, if we work in the system where the neutrinos define the x-axis (θ ν = π/2), as done in the KMS method, in such a way that cos θ = cos ϕ sin θ ℓ ′ , the difference between Dirac and Majorana cases after the angular integration is:  Again, the difference between Dirac and Majorana cases vanishes, consistently with the last computation, since the physics must not depend on the selected reference frame we are working on.Here we anticipate that the main discrepancy with the KMS method is that their analysis sets ϕ = 0, where in this example it is clear that setting ϕ = 0 would lead to different results depending on the selected system, which makes no physical sense.This reason, along with several other arguments, will lead us to conclude that ϕ is not fixed by any kinematic condition and must be integrated over its entire range, as we will discuss in detail later.
Actually, we can do the same computation for the ℓ ′ − γ decay rate (since, even if in the b2b configuration we can relate E ℓ ′ and E ν , the neutrinos and electronphoton distributions are not the same in this kinematic scenario, as shown in eqs.(3.1) and (3.2)) and plot the specific energy distribution for the ℓ ′ − γ and ν − ν pairs, integrating over the inaccessible angle θ; this is shown in figure 4 (just for a τ decay, for simplicity).As we pointed out before, both distributions are different and the Dirac and Majorana cases are, unfortunately, indistinguishable in each case.
This last point is specially important, since the explanation of the impossibility to distinguish the specific neutrino nature in each distribution is different.For the ℓ ′ −γ spectrum, it is a direct application of the DMCT, since we already integrated over all the neutrino variables and we are neglecting neutrino mass contributions.For the ν − ν spectrum, one might think that there should be a difference between Dirac and Majorana distributions, since we are keeping information about the neutrinos energies.This would be true if the neutrinos energies were not equal, so the change E ν ↔ E ν could lead to a difference while computing the Majorana neutrino case.Regrettably, in the b2b case, they are exactly the same (E ν = E ν ), so we can not distinguish between Dirac and Majorana nature from this energy variable alone.This statement is precisely the reason why the difference between Dirac and Majorana distributions can not be independent of the angular variable θ.

Discussions and Consistency Tests
In this section we elaborate on our results and focus on discussing the main reason for the lack of a difference between Dirac and Majorana nature of neutrinos.This outcome, which may differ from initial expectations based on the motivation provided by the KMS method, deserves careful analysis.We emphasize that the specific process under study is not the same as the one used in [1], and the Majorana nature affects in a distinct manner the two of them.Nevertheless this will only influence the dynamics and not the kinematics.So, from now on, we will be focusing on the kinematic analysis, which must be the same considering all our previous assumptions.
Particularly, working in the same system as done in the KMS method, we traced the main difference of this feature with the KMS approach and found that it primarily arises from the ϕ variable treatment.Therefore we want to delve further into this topic and specify the reasons why this variable should not be fixed in the b2b configuration.Additionally, we will outline the results that would be obtained if we do so in our specific process.
It is suggested in the KMS method that ϕ = 0 is a condition fixed by the b2b kinematics.As we mentioned, the b2b scenario is obtained by applying the condition ⃗ p 2 = −⃗ p 3 , or equivalently ⃗ p 4 = −⃗ p 5 .These restrictions affect three of the five phase-space kinematic independent variables as follows: E ν = E ν and Θ ν ν = π.Therefore, the remaining two angular variables must run over all their possible configurations, meaning that we must not set ϕ to any specific value.In other words, the condition ⃗ p 2 = −⃗ p 3 can be achieved for any value of ϕ and not just for ϕ = 0. Nevertheless, beyond these qualitative arguments, we did the quantitative derivation of this assertion in Appendix B, where it is shown explicitly that ϕ is not fixed at all by the b2b restriction.
Finally, the last argument given in the KMS approach is that in the b2b configuration the ν ν and ℓ ′ γ systems define a plane (since they are two independent vectors) and thus ϕ = 0. We also clarify this statement in Appendix B, where we fully agree with the ν ν and ℓ ′ γ systems defining a plane, but show that such plane is independent of the ϕ value, being ϕ = 0 just an allowed specific configuration.This can be seen directly from the plane equation and applied to any selected system.Then ϕ remains an independent variable.Also, as discussed in previous sections, the physics must not depend on the selected reference frame.We already showed that, at least, in two different systems (θ ν = 0 and θ ν = π/2) the difference between Dirac and Majorana distributions vanishes while doing the angular integration properly (not fixing ϕ = 0), which extends to any other.We shall discuss next what happens when assuming the condition ϕ = 0 and will obtain that this difference will be non-vanishing in the system defined by θ ν = π/2, while it will remain zero in the system where θ ν = 0, which is clearly a physical contradiction and again another argument against fixing ϕ in the b2b kinematic configuration.Now, for completeness, we will study the possible modifications to our main results in the case that we set ϕ = 0.In the system defined by θ ν = 0, it is straightforward that the difference in any expression will be just an overall 2π factor after the angular integration, due to the azimuthal symmetry.Also, for the Dirac-Majorana difference, the result will vanish again due to the θ ℓ ′ integration.For the following discussion, it is important to compare what happens now in the KMS system (θ ν = π/2).Replacing the ϕ = 0 condition into eq.(3.3),where now cos θ = cos ϕ sin θ ℓ ′ = sin θ ℓ ′ , the difference between Dirac and Majorana cases is precisely: and will be non-zero after the angular integration, as we shall see next, which is a physical contradiction as we just emphasized.For a complete discussion of this case (ϕ=0), we can now integrate over the energy range and get the angular distribution shown in figure 5a (5c) for the tau (muon) decay, where the difference between Dirac and Majorana cases is obvious.Again, for the experimental observable we must integrate over the inaccessible angle θ (cos θ = sin θ ℓ ′ ), to get the final neutrinos energy distributions for Dirac and Majorana cases: ) This final energy distribution, obtained after integrating over cos θ ℓ ′ , is shown in figure 5b (5d) for the tau (muon) decay, considering the Dirac and Majorana cases, with a clear difference between both of them, given explicitly by the following expression: Therefore, by setting ϕ = 0, it would be possible in principle to distinguish between Dirac and Majorana nature of neutrinos measuring the final energy distribution of the b2b configuration.Then, this could be a really promising result, motivating its search in current and future experiments, being an attractive alternative process to the one studied in [1], due to its larger branching ratio (BR) * .Finally, focusing on this last statement, we would also like to clarify the estimation process of the BR as done in the KMS method, which could lead to confusion for the b2b case.This is done in detail in Appendix C, where we find a BR for the b2b case of the following order: being too suppressed, as expected for such a specific kinematic configuration.
In conclusion, once the appropriate treatment of the phase-space in the b2b configuration is clarified, our approach remains consistent with all the tests carried out and allows a clearer interpretation of the results, leading to the fact that there is no difference produced by the Dirac or Majorana nature of neutrinos in ℓ → ℓ ′ ν νγ, independently of the neutrino mass, once the inaccessible neutrino variables are integrated out.Recalling our comments above eq.(2.9) and below eq.(3.3), any difference between these two cases will be heavily suppressed by squared sterile neutrino masses and mixings, likely preventing their soon observation.We remark that we were neglecting possible neutrino non-standard interactions.It remains to be studied how much these can enhance the sensitivity to the neutrino nature [66], through BR(ℓ = µ, τ ) b2b , while still being consistent with all other constraints.

Summary and Conclusions
In this work we have studied the radiative leptonic lepton-decay ℓ → ℓ ′ ννγ, in the neutrinos mass basis.We have developed our own approach, inspired in the * Unfortunately, if we consider the full angular dependence (without setting ϕ = 0) as we just discussed, the difference between Dirac and Majorana cases vanishes once the angular integration is made, consistently with the result in the last section.method put forward in [1], extending its application to final-state neutrinos of different flavour, in order to distinguish between Dirac and Majorana neutrinos; we have computed its matrix element for the b2b configuration in the decaying lepton rest frame for both cases.
Presumably, this b2b kinematic scenario avoids the constraint imposed by the 'practical Dirac-Majorana confusion theorem' (DMCT) [49] leading to striking differences between the Dirac and Majorana cases, that are not proportional to tiny neutrino masses.Instead, we found that there is no difference between Dirac and Majorana distribution in ℓ → ℓ ′ ννγ once the inaccessible neutrino angle is integrated out.
We discussed in detail the angular treatment, with quantitative and qualitative arguments favoring our conclusions.This turns out to be crucial, since its inaccurate interpretation could lead to very attractive results, even observable in current and near future experiments.However, after careful study, this leads to no difference between Dirac and Majorana energy distributions, for the process under consideration, remaining consistent with all the tests done.
Finally, we wish to emphasize that the idea proposed by [1] is very appealing in order to avoid the DMCT.This fact highlights the necessity to study other types of decays and specific kinematic scenarios within this approach, where angular or energy dependencies could lead to a non-zero difference between the Dirac and Majorana distributions, hopefully observable in current or forthcoming experiments.
• ϕ, the azimuthal angle described between the two planes defined by the ν − ν and ℓ ′ − γ systems, in the rest frame of ℓ.
We can write the differential decay width as 2) where X is the magnitude of three-momentum of ν − ν or ℓ ′ − γ system in the rest frame of ℓ, while β ν (β ℓ ′ ) refers to the magnitude of three-momentum of the ν (ℓ ′ ) in the center-of-momentum frame of the ν ν (ℓ ′ γ) pair.
It is convenient to rewrite the differential width in terms of E ν and E ν in order to obtain the energy spectrum of the neutrinos.Thus, we can do the variables change: where Θ ν ν is the angle between the three-momenta of both neutrinos, in the rest frame of ℓ.Also, we can obtain the energy spectrum for the final charged lepton and the photon with the following change of variables: where Θ ℓ ′ γ is the angle between ℓ ′ and γ, in the rest frame of ℓ.
Finally, neglecting all the final-state masses, as a good approximation, we get for the neutrinos differential decay distribution: where ϵ = 1(2) for Dirac (Majorana) neutrinos.Here E ℓ ′ and E γ must be written in terms of E ν and E ν , according to the energy-momentum conservation law.Note also that we are taking into account all the possible neutrino mass final states and the sum extends over all energetically allowed neutrino pairs.The 1/2 factor that appears in the Majorana case has two different origins.For the k = j case, it is due to the presence of indistinguishable fermions in the final state; while for k ̸ = j, it arises because of double counting, since the sum is not restricted to j ≤ k.
Meanwhile, for the differential decay distribution involving the charged lepton and photon energies we obtain: where E ν and E ν must be written in terms of E ℓ ′ and E γ according to the energymomentum conservation law.These, in principle, are two different spectra and will be so in any specific kinematic configuration.

B Back-to-back configuration
As written above, it is convenient to describe our phase space variables in the rest frame of the decaying particle to avoid any confusion.First, we are going to denote p ′ i as the momentum of the i−particle in the corresponding center of mass frame for the relevant particle pair (ν − ν or ℓ ′ − γ) and p i the corresponding momentum in the rest frame of the decaying particle.Now, following Fig. [4] of ref. [1], both momenta, p ′ i and p i , are related by a Lorentz boost in the ẑ direction.We define our boost in the ẑ direction for the 4-momentum p ′ ν and p ′ ν as follows: For p ′ 4 and p ′ 5 we use the Lorentz transformation: where we employ the definition of X, s ν ν and s ℓ ′ γ from the previous appendix.
In general, we can write the corresponding 4-momentum in the rest frame of ℓ as follows: Finally, we can apply the b2b constraint ⃗ p 2 = −⃗ p 3 or equivalently, due to energy momentum conservation, ⃗ p 4 = −⃗ p 5 .In this kinematic scenario it is easy to show that X = 0, which is consistent with the fact that, in the b2b configuration, the boosts in eqs.(B.1,B.2) are exactly the identity matrix, i.e., the center of mass frame of the ν − ν and ℓ ′ − γ systems coincides with the decaying lepton rest frame.Then, in the b2b case the three-momentum of the final-state particles can be written as follows Here it is essential to emphasize that the above equations fulfill the b2b constraint (⃗ p 2 = −⃗ p 3 and ⃗ p 4 = −⃗ p 5 ) for all possible values of (θ ℓ ′ , θ ν , ϕ) and not only when ϕ = 0. Thus, we showed cleverly that ϕ = 0 is not a constraint imposed by the b2b kinematics, as suggested in [1] (see also [57]), and it needs to be integrated over its full range.
Another important result comes from the definition of the angle θ that appears in the amplitude, which is the angle between the neutrino and the charged lepton.From that definition, it is straightforward to compute its explicit form, in terms of the angles θ ℓ ′ , θ ν and ϕ, for the b2b configuration.Using the above expressions for the three-momentum in the b2b case we obtained: that shows the specific dependence of θ on θ ν , θ ℓ ′ and ϕ.Other relations resulting from ⃗ p 2 = −⃗ p 3 and ⃗ p 4 = −⃗ p 5 (neglecting the mass of the final-state particles) are Finally, since in this b2b configuration the ⃗ p 2 (−⃗ p 3 ) and ⃗ p 4 (−⃗ p 5 ) are two independent vectors, they can always span a plane, i.e. they can always form a basis of a two-dimensional space.This argument is used in the KMS method to claim that ϕ = 0.For completeness we work on this subject below and demonstrate that, even it is certainly true that these two vectors span a plane, this condition does not fix ϕ = 0, as we just showed before.
How can we describe that plane?Since in the ℓ rest frame, these vectors start from the origin of the coordinate system (x 0 , y 0 , z 0 ) = (0, 0, 0), then the plane spanned by the vectors ⃗ p 2 and ⃗ p 4 is given by the well-known equation: where λ and ν are just the parameters of the plane-equation (−∞ < λ, ν < +∞).
The following computation can be done in any selected reference frame, but for this specific discussion we keep working in the system where θ ν = π/2 as done in the KMS method, where the plane equation can be put in the general form, giving rise to (after a fast simplification): Then it is completely clear that the b2b condition is satisfied for each value of ϕ and the corresponding plane spanned by the vectors ⃗ p 2 and ⃗ p 4 is given by eq.(B.12).This reaffirms that ϕ = 0 is not a condition fixed by the b2b scenario, and instead ϕ remains as an independent variable that must be integrated over its full range.To illustrate this plane condition, we show in Fig. 6 various examples with different ϕ and θ ν : • Diagram (a): For ϕ = 0 we get y = 0, which means that the plane is precisely the x − z plane, which is in agreement with the KMS method, but is not the only option allowed by the kinematics.
• Diagram (b): For θ ℓ ′ = π/2 we get z = 0, which means that the plane is precisely the x − y one, a configuration completely allowed by the kinematics.
And you can keep going with all the possible configurations, as diagram (c), etc.All of them are in agreement with the b2b constraints and with the fact that the two final independent vectors span a plane, checking again that ϕ = 0 is just an allowed configuration but not a condition fixed by the b2b kinematics.

C Branching ratio computation
First we will compute the BR of this b2b kinematic scenario in a way that might seem correct at first glance, but that -without the right considerations-can lead to erroneous conclusions.Finally we will discuss this problem in detail and the correct way to estimate this observable.
For a first consistency test, the total BR(ℓ → ℓ ′ ν νγ) for the general kinematic configuration was computed, being in perfect agreement with those reported by [61], giving us a corroboration that our procedure was correct.Now, as a first attempt, one might be tempted to estimate the b2b BR by integrating over all the energy and angular configurations of the differential decay rates evaluated in this kinematic case.Essentially, for the case ϕ = 0, integrating over the remaining energy dependence of eqs.(3.9, 3.10): Doing this, and cutting off photons below an energy threshold of 10 MeV in the decaying-lepton rest frame, given by the experimental resolution at Belle [62] (achievable at Belle-II [63]) * , we obtain: which, in principle, is a result that could motivate even more its search and reflect the advantages of this specific process (ℓ → ℓ ′ ν νγ), since we do not have to deal with hadronic form factors and the computed BR is much larger than the ones estimated in [1] for B decays (B ↔ ≈ 10 −12 ) and related processes.
(a) Nevertheless, the estimated BRs (C.2, C.3) seem troublesome.First of all, they are different for Dirac and Majorana cases (because we used the condition ϕ = 0, as stressed), which disagrees with the DMCT theorem, as we must integrate over the neutrinos variables to calculate them.Even so, this can be easily corrected, just considering the full range of variation of ϕ.
The main problem is quite clear: Since the total BR(ℓ → ℓ ′ ν νγ) is of the order O(10 −2 ), it is hard to believe that a specific kinematic configuration, such as the b2b, is only two orders of magnitude more suppressed than the general case.Then, for a complete discussion of this problem we will comment on the specific reason why this BR estimation is wrong, and also do the right computation for this case, which can be applied for any other specific kinematic configuration in which one (or more) of the continuous phase space variables is(are) fixed to specific value(s).
The width of an N-body decay can be seen as the hypervolume of the phase space weighted by the dynamic condition (squared amplitude) of the specific process.This hypervolume is determined by the specific range of all the continuous phase space independent variables, which is specified by the minimum and maximum values they could take according to energy-momentum conservation.If one (or more) of these variables take(s) a fixed value in a specific kinematic configuration, an integration over a zero-range variable has to be done in order to compute the theoretical BR, leading to a vanishing contribution for this specific case.
In other words, once a continuous phase space variable is fixed, the phase-space hypervolume is reduced to a phase-space hypersurface, meaning that in that case the purely theoretical BR estimation will be zero for that configuration.We note this is congruent with the intuitive notion of obtaining a null probability for a unique configuration among all the continuous (infinite) possibilities.
This does not mean that the differential decay rate is zero in that case.Actually, we can compute without further difficulties any differential distribution as long as the fixed variables are not integrated.Then, the main problem of the estimated BR is that we integrated over the already evaluated differential decay rate, leading to a number that does not have a probability interpretation, since for the correct theoretical BR computation, we first need the differential decay rate for the general case and then to integrate over all phase-space variables, which range will be fixed by the specific kinematic scenario and the energy-momentum conservation.
In particular our notation first introduced in eqs.(3.1), was precisely motivated to avoid this possible confusion.It provides evidence that, once the decay rate is already evaluated in a specific kinematic scenario, we can not integrate over the kinematic variables fixed by the b2b condition and interpret the result as a probability, specifically as the BR of the b2b case, that could lead to a larger value than the real one.
Finally, we know that the experimental resolution is finite and thus experimentally we can not have a strictly fixed variable.Then, it is well-defined to estimate a non-zero BR for the b2b configuration, considering a small range of variation for the theoretically fixed variables, according to the experimental resolution.
Then, to estimate the real BR for the experimental b2b case, without fixing ϕ and using a proper method for this estimation, we have to do the following: First we integrate over the remaining neutrino's variables (θ ν and ϕ) in eq.(A.6), which are not fixed by any kinematic condition, leading to a decay rate that depends only on the electron and photon energy, together with the angle between them dΓ dE ℓ ′ dEγd cos Θ ℓ ′ γ , in the general kinematic case.Then, using this differential decay rate, we can apply the experimental energy and angular resolution * to integrate over the "pseudo" b2b case, i.e. an infinitesimal phase-space region that will be indistinguishable from the theoretical b2b case by the experiment, as shown in figure 7, where the energy of the final charged-lepton and the photon are equal, up to the energy resolution (E ℓ ′ −∆E ≤ E γ ≤ E ℓ ′ +∆E), and the angle between them is π, up to the angular resolution (π −∆θ ≤ θ ℓ ′ γ ≤ π).
Since we are not evaluating the differential decay rate in a specific kinematic configuration in any moment, we can safely interpret this result as an occurrence probability.Furthermore, since the integration region contains the theoretical b2b case, the corresponding BR obtained after the integration must therefore be larger than the theoretical b2b case alone.Doing these, cutting off photons below an energy threshold of 10 MeV in the decaying-lepton rest frame, we estimate the following BR for the experimental b2b case as follows: As we can see, these branching ratios are many orders of magnitude smaller than the ones obtained, using the KMS method, in (C.2) and (C.3).Also, the BR calculated is the same for Dirac and Majorana neutrinos, in agreement with the DMCT theorem, being all consistent with our previous discussions.

Figure 2 :
Figure 2: b2b kinematic configuration in the initial charged-lepton rest frame.

Figure 3 :
Figure 3: Lowest order Feynman diagrams for radiative leptonic lepton-decay in the Majorana neutrino case.

Figure 4 :
Figure 4: Energy spectra for charged lepton-photon and neutrinos in τ → ℓ ′ ν νγ decays, restricted to the b2b case.They are identical for Dirac and Majorana cases.
Comparison of angular distributions between Dirac and Majorana cases in the b2b scenario (ℓ = τ ).Comparison of energy distributions between Dirac and Majorana cases in the b2b scenario (ℓ = τ ).Comparison of angular distributions between Dirac and Majorana cases in the b2b scenario (ℓ = µ).
Comparison of energy distributions between Dirac andMajorana cases in the b2b scenario (ℓ = µ).

Figure 5 :
Figure 5: Comparison of Dirac and Majorana distributions.

Figure 7 :
Figure 7: Experimental b2b kinematic configuration in the initial charged-lepton rest frame.