1 Introduction

Observations of neutrino oscillation [1, 2] have established the fact that neutrinos have non-zero but tiny masses and the flavor neutrinos (\(\nu _\ell \) with \(\ell =e, \mu , \tau \)), which only participate in the weak interactions, are linear combinations of mass eigenstates (\(\nu _i\) with \(i=1, 2, 3\) having masses \(m_i\)),

$$\begin{aligned} |\nu _\ell \rangle = \sum _i U_{\ell i}^* |\nu _i \rangle , \end{aligned}$$
(1)

where U denotes the unitary \(3\times 3\) lepton mixing matrix, also called the PMNS matrix [3, 4]. At present [5], we know neither the values of the individual masses of the neutrinos nor the mechanism that gives rise to these tiny masses (see Refs. [6, 7] for recent reviews of various neutrino mass models). In the Standard Model of particle physics (SM) neutrinos are massless and get produced via weak interaction such that all neutrinos are left-handed and all antineutrinos are right-handed. However, as the neutrinos are not massless, handedness is not a good quantum number for neutrinos, i.e. in a frame moving faster than the neutrino (or antineutrino) the handedness will be reversed. Nevertheless, direct production of right-handed neutrinos and left-handed antineutrinos has not yet been observed in any experiment, consistent with the SM allowed \(V-A\) interactions. Therefore, it is currently unknown whether the right-handed neutrinos have the same mass as the left-handed neutrinos or not. Most certainly, detection of the right-handed neutrinos and the left-handed antineutrinos would require invoking some sort of new physics (NP) which we currently do not know. To complicate the matter further, neutrinos being devoid of any electric or colour charge are the only known fermions which can possibly be their own antiparticles, i.e. they can be Majorana fermions [8,9,10] unlike all the others which are Dirac fermions.Footnote 1 In most of the interesting mass models for neutrinos, all of which involve different NP setups, the Majorana nature is assumed. Like with many other aspects of neutrinos, this is yet to be decisively decided by experiments.

Various methodologies that have been proposed to probe the Majorana nature of neutrinos can be broadly grouped into the following two categories, depending on the type of process being considered.

  1. 1.

    Propagator probes or neutrino mediated probes: Well known examples of such processes are neutrinoless double-beta decay (\(0\nu \beta \beta \)) [11], neutrinoless double-electron capture [12], neutrinoless muon (\(\mu ^-\)) to positron (\(e^+\)) conversion [4, 13], \(\nu \) exchange force (\(\nu \) Casimir force) [14,15,16,17,18,19,20,21] or lepton number violating massive sterile neutrino exchange [22], which involve one or more neutrino(s) and/or antineutrino(s) as propagator(s). The lepton number violating processes involve helicity flip and are, therefore, always proportional to the unknown mass of the propagating Majorana neutrino. In addition, all the nuclear processes explored as propagator probes to infer the nature of sub-eV active neutrinos suffer from large theoretical uncertainty in the nuclear matrix elements. Although \(0\nu \beta \beta \) searches have thus far not yielded any signal [23], the experimental searches for it are going on. Nevertheless, the associated issues mentioned above underline the importance of alternative methods.

  2. 2.

    Initial / final state probes: These involve one (or more) neutrino(s) and/or antineutrino(s) in the initial state or the final state. Some previously explored processes under this category are the neutrino-electron elastic scattering [24,25,26,27], neutrino-nucleon scattering [24], neutrino-nuclei scattering [28], 2-body processes \(\gamma ^* \rightarrow \nu \,{\overline{\nu }}\) [29] that try to explore the electromagnetic properties of neutrinos [30], \(Z \rightarrow \nu \, {\overline{\nu }}\) [31], \(e^+ \, e^- \rightarrow \nu \, {\overline{\nu }}\) [32], 3-body processes \(K^+ \rightarrow \pi ^+\,\nu \,{\overline{\nu }}\) [33], \(e^+ \, e^- \rightarrow \gamma \,\nu \,{\overline{\nu }}\) [34], \(e\,\gamma \rightarrow e \, \nu \, {\overline{\nu }}\) [35], radiative emission of neutrino pair [36], and 4-body B meson decay \(B^0 \rightarrow \mu ^+ \, \mu ^- \, \nu _\mu \, {\overline{\nu }}_\mu \) [37]. In contrast with the propagator probes which essentially directly test the Majorana mass term, the initial/final state probes try to examine the quantum mechanical identicalness of Majorana neutrino and antineutrino via quantum statistics.

    Unlike the neutrino mass suppressed processes for propagator probes, the initial/final state probes involve SM allowed processes and therefore can be produced relatively easily in experiments. The challenge in probing the Dirac or Majorana nature of neutrinos using such probes is to infer the momenta or energies of the missing neutrinos which is always experimentally challenging.

The quantum statistics of Majorana neutrino and antineutrino stems from their quantum mechanical identical nature, and it is an intrinsic property of the Majorana neutrino just like its spin and charge. For a final state containing a pair of Majorana neutrino and antineutrino of the same flavor, the Fermi–Dirac statistics asserts that the corresponding transition amplitude is always anti-symmetrized with respect to their exchange. This anti-symmetrization does not depend on the size of the mass of neutrino (\(m_\nu \)). Therefore, there is no reason a priori to think that the difference between Dirac and Majorana neutrinos via such initial/final state probes would necessarily be dependent on \(m_\nu \), a result that apparently (but not necessarily) contradicts the “practical Dirac Majorana confusion theorem” (DMCT) [24, 29]. A brief overview of the literature, in this context, comparing two-body decays, three-body decays as well as four-body decays is given in Ref. [37].

The various initial/final state probes can, in general, be related to one another via crossing symmetry [25]. For ease of discussion and without loosing any generality we concentrate in this paper on providing a generalised framework for final state probes, i.e. using processes in which a pair of neutrino and antineutrino would appear in the final state. Furthermore, we relax the helicity considerations, i.e. we do not confine ourselves to left-handed neutrinos and right-handed antineutrinos alone, since the involved interactions will take care of the helicities by default. Unlike considering only SM allowed interactions as done in Ref. [37] here we consider NP effects in a model-independent manner and point out their effects as well.

Before providing the layout of our paper, we would like to emphasise that the active neutrinos and antineutrinos in all the final states considered in this paper are invisible. In many NP models it is possible to have other invisible particles such as the heavy sterile neutrinos, dark matter particles or some long-lived light particles (e.g. axion-like) that decays outside the detector volume. Therefore, observation of decays into invisible states with rates modified comparing to the SM predictions usually is not sufficient by itself to determine the source of the modification, being it different types of neutrinos or entirely novel particles. The analysis we present in this paper is limited by the assumptions that neutrinos are only invisible final state but it is still useful for several reasons.

  • If new invisible particles are massive (even relatively light but with masses significantly exceeding neutrino masses), they could have noticeable effects such as significant reduction in phase space, or presence of some resonance peaks corresponding to on-shell production of the long-lived particles. Corresponding modifications of visible final state particle spectra could help in distinguishing such cases from the neutrino production. Such possibilities are well explored in the literature, see e.g. [38] for a recent detailed review of invisible K meson decays in various SM extensions such as the 2-body decay \(K \rightarrow \pi \, X_{\text {inv}}\) where the invisible particle \(X_{\text {inv}}\) could be a Higgs mixed scalar, dark photon, or axion-like particle.

  • Current experiments are in general in agreement with the SM predictions, therefore the bounds on non-standard neutrino couplings which we discuss in the paper come from comparing their effects with the size of experimental errors of relevant observables. If there are more invisible states in a given NP model, even very light ones, barring fine tuning they compete with neutrino interactions in saturating the errors, tightening the bounds. Thus the limits on NP neutrino couplings which we discuss can be considered as (again barring fine tuning and potential cancellations) being on the conservative side.

  • Finally, in many SM extensions neutrinos remain the only ultra-light and weakly interacting (thus “invisible”) particles. In such cases methodology proposed in the paper can be directly used to constrain NP neutrino couplings for a chosen BSM model.

Our paper is organised as follows. Following Introduction, in Sect. 2 we present a model-independent general formalism for processes involving a neutrino and antineutrino pair of the same flavor in the final state and we provide a generalised statement of DMCT in context of these processes. In Sect. 3 we discuss how the effects of quantum statistics and symmetry properties of transition amplitudes can be utilised to probe NP interactions in neutrino sector and potentially lead to substantial difference between Dirac and Majorana neutrinos. In the following Sect. 4 we illustrate the general results with two suitably chosen two-body and three-body neutral current decays, showing how the bounds on the NP couplings depend on the assumed neutrino nature. The possibility of distinguishing the Dirac and Majorana neutrinos by measuring the visible particles energy spectrum in the three-body decays is discussed in Sect. 5. Finally we conclude in Sect. 6 highlighting the various salient aspects of our study.

2 Model-independent formalism for processes containing \(\varvec{\nu \, {\overline{\nu }}}\) in the final state

2.1 Choice of process

Let us consider a general process with a neutrino and an antineutrino of the same flavor in the final state, say

$$\begin{aligned} X (p_X) \rightarrow Y(p_Y) \, \nu (p_1) \, {\overline{\nu }}(p_2), \end{aligned}$$

where X, Y can be single or multi-particle states, Y can also be null, the contents of X and Y (if it exists) can only be any visible particle and the 4-momenta \(p_X, p_Y\) are assumed to be well measured so that one can unambiguously infer the total missing 4-momentum of \(\nu \, {\overline{\nu }}\), \(p_{\text {miss}} = p_1 + p_2\). So the 4-momentum of X must either be fixed by design of the experiment (e.g. X might be a particle produced at rest in the laboratory or be the constituent of a collimated beam of known energy or it could consist of two colliding particles of known 4-momenta) or the 4-momentum of X be inferred from the fully-tagged partner particle with which it is pair-produced. The final state Y should not contain any additional neutrinos or antineutrinos. The process could be a decay or scattering depending on whether X is a single particle state or two particle state. Some actual processes that satisfy such criteria are \(e^+\, e^- \rightarrow \nu \, {\overline{\nu }}\), \(Z \rightarrow \nu \, {\overline{\nu }}\), \(e^+\, e^- \rightarrow \gamma \, \nu \, {\overline{\nu }}\), \(K \rightarrow \pi \, \nu \, {\overline{\nu }}\), \(B \rightarrow K \, \nu \, {\overline{\nu }}\), \(B^0 \rightarrow \mu ^+ \, \mu ^- \, \nu _\mu \, {\overline{\nu }}_\mu \), \(J/\psi \rightarrow \mu ^+ \, \mu ^- \, \nu _\mu \, {\overline{\nu }}_\mu \), \(H \rightarrow \tau ^+ \, \tau ^- \, \nu _\tau \, {\overline{\nu }}_\tau \) etc.

Schematically the process \(X \rightarrow Y \, \nu \, {\overline{\nu }}\) is drawn in Fig. 1. Since Majorana neutrino and antineutrino are quantum mechanically identical, there is an additional diagram with exchanged 4-momenta \(p_1 \leftrightarrow p_2\) in the Majorana case as shown in Fig. 1. The effect of this exchange diagram in Majorana case can, in many instances but not always, be absorbed into a single diagram like the one for the Dirac case but with a different effective vertex. In Sect. 4 we shall encounter a few examples where the effective vertex for Majorana case is different from Dirac case. An example process where we can not absorb the direct and exchange diagram contributions to the effective vertex factor is \(B^0 \rightarrow \mu ^+ \, \mu ^- \, \nu _\mu \, {\overline{\nu }}_\mu \) as discussed in Ref. [37].

Additionally, it should be noted that in this work we assume that in any of the processes under our consideration, the effect of measurements should not destroy the identical nature of Majorana neutrino and antineutrino. This is akin to putting the constrain that in a double-slit experiment meant to observe the interference of light, no measurement should identify through which slit the photon has passed.

Fig. 1
figure 1

A cartoon (pseudo-Feynman diagrams) for the processes \(X \rightarrow Y \, \nu \, {\overline{\nu }}\), with both Dirac and Majorana neutrino possibilities. The blob in these diagrams represents the effective vertex and it includes both SM and NP contributions. Although Majorana antineutrino is indistinguishable from Majorana neutrino, we keep using the notation of \({\overline{\nu }}\) for antineutrino and \(\nu \) for neutrino simply as a book-keeping device

Full size image

2.2 Origin of observable difference between Dirac and Majorana neutrinos

The transition amplitude is, in general, dependent on all the 4-momenta of participating particles. For brevity of expression and without loss of generality, we denote the transition amplitude by only mentioning the \(p_1\), \(p_2\) dependence. For Dirac neutrinos, the transition amplitude can be written as,

$$\begin{aligned} {\mathscr {M}}^D = {\mathscr {M}}(p_1, p_2), \end{aligned}$$
(2)

while for Majorana case the amplitude is anti-symmetrized with respect to the exchange of the Majorana neutrino and antineutrino which are quantum mechanically identical fermions,

$$\begin{aligned} {\mathscr {M}}^M = \dfrac{1}{\sqrt{2}} \Big (\underbrace{{\mathscr {M}}(p_1, p_2)}_{\text {Direct amplitude}} - \underbrace{{\mathscr {M}}(p_2, p_1)}_{\text {Exchange amplitude}} \Big ), \end{aligned}$$
(3)

where \(1/\sqrt{2}\) takes care of the symmetry factor. Note that the amplitudes of Eqs. (2) and (3) do not necessarily assume SM interactions, they can involve NP effects as well, and hence they include the most general structures of the amplitude that are allowed by Lorentz invariance. A general analysis that specifically exploits the symmetry properties of the amplitude in context of generic NP possibilities is given in Sect. 3.

The difference between Dirac and Majorana cases that can possibly be probed is obtained after squaring the amplitudes and taking their difference, which is given by,

$$\begin{aligned} \left| {\mathscr {M}}^D\right| ^2 -\left| {\mathscr {M}}^M\right| ^2&= \frac{1}{2} \Bigg (\underbrace{\left| {\mathscr {M}}(p_1, p_2)\right| ^2}_{\text {Direct term}} - \underbrace{\left| {\mathscr {M}}(p_2, p_1)\right| ^2}_{\text {Exchange term}}\Bigg ) \nonumber \\&\quad + \underbrace{\text {Re}\left( {\mathscr {M}}(p_1, p_2)^* \, {\mathscr {M}}(p_2, p_1)\right) }_{\text {Interference term}}. \end{aligned}$$
(4)

From Eq. (4) it is easy to conclude that there are essentially two sources of difference between Dirac and Majorana cases,

  1. 1.

    unequal contributions from “Direct term” and “Exchange term”, i.e.

    $$\begin{aligned} \underbrace{\left| {\mathscr {M}}(p_1, p_2)\right| ^2}_{\text {Direct term}} \ne \underbrace{\left| {\mathscr {M}}(p_2, p_1)\right| ^2}_{\text {Exchange term}}, \end{aligned}$$
    (5)
  2. 2.

    non-zero contribution from the “Interference term”, i.e.

    $$\begin{aligned} \underbrace{\text {Re}\left( {\mathscr {M}}(p_1, p_2)^* \, {\mathscr {M}}(p_2, p_1)\right) }_{\text {Interference term}} \ne 0. \end{aligned}$$
    (6)

However, note that in the case when no individual information about \(\nu \, {\overline{\nu }}\) are either known or deducible, the only difference between Dirac and Majorana cases that can be experimentally accessed is obtained after full integration over \(p_1\) and \(p_2\) which gives,

$$\begin{aligned}&\iint \left( \left| {\mathscr {M}}^D\right| ^2 -\left| {\mathscr {M}}^M\right| ^2 \right) \textrm{d}^4 p_1 \, \textrm{d}^4 p_2 \nonumber \\&\quad = \iint \underbrace{\text {Re}\left( {\mathscr {M}}(p_1, p_2)^* \, {\mathscr {M}}(p_2, p_1)\right) }_{\text {Interference term}} \textrm{d}^4 p_1 \, \textrm{d}^4 p_2, \end{aligned}$$
(7)

where we have used the fact that although, in general, the “Direct” and “Exchange” terms may differ, when we fully integrate over the 4-momenta of neutrino and antineutrino we do get

$$\begin{aligned} \iint \underbrace{\left| {\mathscr {M}}(p_1, p_2)\right| ^2}_{\text {Direct term}} \textrm{d}^4p_1 \, \textrm{d}^4 p_2 = \iint \underbrace{\left| {\mathscr {M}}(p_2, p_1)\right| ^2}_{\text {Exchange term}} \textrm{d}^4p_1 \, \textrm{d}^4 p_2,\nonumber \\ \end{aligned}$$
(8)

as \(p_1\) and \(p_2\) act as dummy variables in this case.

2.3 “Dirac–Majorana Confusion Theorem” for the neutrino interactions in the SM

In most of the experimental scenarios, especially true for processes of the form \(X \rightarrow Y \, \nu \, {\overline{\nu }}\), information about individual neutrino momenta is not available. In such a case the difference between the Dirac and the Majorana neutrinos is given only by the integrated interference term in Eq. (7). If we further consider only SM interactions without any NP contributions, then the \(V-A\) interaction of SM would always produce left-handed neutrino and right-handed antineutrino, even if the mass of the neutrino and antineutrino is considered to be non-zero. In such a case the evaluation of the squared Feynman diagram for the “Interference term” would, in general, necessarily involve two helicity flips which would make it proportional to \(m_\nu ^2\), i.e.

$$\begin{aligned} \underbrace{\text {Re}\left( {\mathscr {M}}^\text {SM}(p_1, p_2)^* \, {\mathscr {M}}^\text {SM}(p_2, p_1)\right) }_{\text {Interference term in the SM}} \propto m_\nu ^2. \end{aligned}$$

Thus, if only SM interactions are considered and one fully integrates over the neutrino and antineutrino 4-momenta, the difference between Dirac and Majorana cases would be proportional to \(m_\nu ^2\). This can be considered as a generalised statement of the “practical Dirac Majorana confusion theorem” (DMCT).

The above statement still holds taking into account that in many cases experiments can measure the total 4-momentum of the neutrino and antineutrino pair in the final state, \(p_1+p_2\), equal to the difference between initial and final 4-momenta of the visible particles. The knowledge of \(p_1+p_2\) only (not \(p_1\) and \(p_2\) individually), gives access to that part of the full phase space for the neutrino pair which is obviously symmetric under the exchange \(p_1 \leftrightarrow p_2\). Thus argument about symmetric integration and corresponding cancellation between “Direct” and “Exchange” terms in Eq. (8) still holds, even if the phase space integral is not done over the full range of \(p_1,p_2\). In such a scenario, all the cross sections or decay rate distributions expressed in terms of parameters related to the 4-momenta of visible particles will always be identical if pure SM interactions are assumed (up to the neutrino mass suppressed effects). However, as we will show in the following sections, the distributions of visible particles may differ if one also allows non-standard neutrino interactions.

The remaining possibility of experimental observation of the difference between the Dirac and Majorana neutrinos in the SM is to reconstruct or at least constrain the neutrino 4-momenta or some of their components such as the energies, so that full integration over the neutrino and antineutrino 4-momenta is not necessary and one can directly use Eq. (4), without worrying about the cancellation in Eq. (8). This, in general, is difficult as neutrinos are invisible in detectors and their 4-momenta can be only inferred from the kinematics of visible particles. Still, it may be possible for some particular configuration(s) of them, however occurring only for a small fraction of all kinematically allowed final states.

Relevant example of such “special” kinematic scenario has been proposed in Ref. [37], where authors consider the 4-body decay of the form \(X \rightarrow y_1\, y_2\, \nu \, {\overline{\nu }}\) with \(y_1\) and \(y_2\) flying away back-to-back with equal magnitudes of 3-momenta in the rest frame of X. In such case, the \(\nu \) and \({\overline{\nu }}\) also fly away back-to-back with equal energies \(E_\nu = (m_X - E_{y_1} - E_{y_2})/2\), where \(E_{y_1}\), \(E_{y_2}\) denote the energies of \(y_1\) and \(y_2\) respectively and \(m_X\) is the mass of parent particle X. This is discussed in detail in Ref. [37] taking the example of \(B^0 \rightarrow \mu ^- \, \mu ^+ \, \nu _\mu \, {\overline{\nu }}_\mu \). As the authors point out, other meson decays such as those of neutral K, or D, or \(J/\psi \), or \(\varUpsilon (n\text {S})\) decays as well as Higgs decays might also be useful for such studies. The difference between Dirac and Majorana cases does not necessarily depend on the size of the neutrino mass, and hence might be feasible to discover in future experiments.

3 New physics scenarios and practical DMCT

There is no reason a priori for the “practical DMCT” to hold, if NP contributions in the neutrino interactions are allowed, as in this case Eq. (4) “Direct” and “Exchange” terms in general do not need to cancel each other. To illustrate it more clearly using symmetry properties of the transition amplitude, let us assume that some (yet unknown) NP at high energy modifies the low energy effective neutrino interactions we want to explore in processes of the form \(X \rightarrow Y \, \nu \, {\overline{\nu }}\), and the SM singlet right-handed neutrinos and left-handed antineutrinos might as well participate in such interactions.

3.1 DMCT from the perspective of exchange symmetry

Although process specific analysis using exchange symmetry is found in the literature, a process-independent general analysis highlighting the difference between Dirac and Majorana neutrinos would be useful in context of NP interactions. Below we provide such an analysis.

The amplitude for \(X(p_X) \rightarrow Y(p_Y) \, \nu (p_1) \, {\overline{\nu }}(p_2)\) for Majorana neutrinos must be still anti-symmetric with respect to the 4-momentum exchange \(p_1 \leftrightarrow p_2\), as was noted in Eq. (3). The Dirac amplitude of Eq. (2) could, in principle, be split into two parts, one symmetric and the other anti-symmetric under the exchange \(p_1 \leftrightarrow p_2\), i.e.

$$\begin{aligned} {\mathscr {M}}^D = {\mathscr {M}}(p_1, p_2) = {\mathscr {M}}_{\text {symm}}(p_1, p_2) + {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2), \nonumber \\ \end{aligned}$$
(9)

where by definition we have

$$\begin{aligned} {\mathscr {M}}_{\text {symm}}(p_1, p_2)&= \frac{1}{2} \Big ({\mathscr {M}}(p_1, p_2) + {\mathscr {M}}(p_2, p_1) \Big ) \nonumber \\&= {\mathscr {M}}_{\text {symm}}(p_2, p_1), \end{aligned}$$
(10a)
$$\begin{aligned} {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2)&= \frac{1}{2} \Big ({\mathscr {M}}(p_1, p_2) - {\mathscr {M}}(p_2, p_1) \Big ) \nonumber \\&= -\,{\mathscr {M}}_{\text {anti-symm}}(p_2, p_1), \end{aligned}$$
(10b)

so that the Majorana amplitude of Eq. (3) is automatically given by,

$$\begin{aligned} {\mathscr {M}}^M = \sqrt{2} \, {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2). \end{aligned}$$
(11)

With the decomposition of amplitudes as shown in Eqs. (9) and (11), the difference between Dirac and Majorana neutrinos is given by,

$$\begin{aligned} \left| {\mathscr {M}}^D\right| ^2 -\left| {\mathscr {M}}^M\right| ^2&= \left| {\mathscr {M}}_{\text {symm}}(p_1, p_2) + {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2)\right| ^2 \nonumber \\&\quad - 2 \left| {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2)\right| ^2\,,\nonumber \\&= \left| {\mathscr {M}}_{\text {symm}}(p_1, p_2)\right| ^2 - \left| {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2)\right| ^2 \nonumber \\&\quad + 2\, \text {Re}\left( {\mathscr {M}}_{\text {symm}}(p_1, p_2)^* \, {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2)\right) \,. \end{aligned}$$
(12)

If one were to do full integration over the neutrino and antineutrino 4-momenta, the interference term in Eq. (12) vanishes as it is anti-symmetric under the \(p_1\leftrightarrow p_2\) exchange and one is left with

$$\begin{aligned}&\iint \left( \left| {\mathscr {M}}^D\right| ^2 -\left| {\mathscr {M}}^M\right| ^2 \right) \textrm{d}^4p_1 \, \textrm{d}^4 p_2 \nonumber \\&\quad =\iint \left( \left| {\mathscr {M}}_{\text {symm}}(p_1, p_2)\right| ^2 - \left| {\mathscr {M}}_{\text {anti-symm}}(p_1, p_2)\right| ^2\right) \textrm{d}^4p_1 \, \textrm{d}^4 p_2, \end{aligned}$$
(13)

which does not necessarily vanish in presence of NP interactions. Note that Eqs. (7) and (13) are equivalent to one another, but Eq. (13) highlights the symmetry properties that were not explicitly evident in Eq. (7).

3.2 Processes with \(\varvec{\nu \,{\overline{\nu }}}\) produced via neutral current interactions

As an example, let us consider NP effects on neutrino-antineutrino neutral current interactions. In such a case, the amplitude would have the generic form,

$$\begin{aligned} {\mathscr {M}}(p_1,p_2) = \sum _i {\mathcal {S}}_i \, \Big [ {\overline{u}}(p_1) \, \varGamma _i \, v(p_2) \Big ], \end{aligned}$$

where \({\mathcal {S}}_i\) denotes the parts of the amplitude which do not depend on \(p_1\) and \(p_2\) or contain algebraic expressions such as \((p_1 + p_2)^2\) which are symmetric under \(p_1 \leftrightarrow p_2\) exchange, and \(\varGamma _i\) could be \({\textbf{1}}\), \(\gamma ^5\), \(\gamma ^\alpha \), \(\gamma ^\alpha \, \gamma ^5\), \(\sigma ^{\alpha \beta }\) or \(\sigma ^{\alpha \beta } \, \gamma ^5\), which respectively correspond to scalar, pseudo-scalar, vector, axial-vector, tensor and axial-tensor interactions. For Majorana neutrinos it is well known that \({\overline{u}}(p_1) \, \varGamma \, v(p_2) = - \, {\overline{u}}(p_2) \, C \, \varGamma ^T \, C^{-1} \, v(p_1)\) [39] where C denotes the charge conjugation operator. Using this it is easy to show that \({\mathscr {M}}_{\text {symm}}(p_1,p_2)\) gets contributions only from vector, tensor and axial-tensor interactions while \({\mathscr {M}}_{\text {anti-symm}}(p_1,p_2)\) gets contributions from scalar, pseudo-scalar and axial-vector interactions alone, i.e.

$$\begin{aligned}&{\mathscr {M}} (p_1,p_2) = \sum _i {\mathcal {S}}_i \, \Big [ {\overline{u}}(p_1) \, \varGamma _i \, v(p_2) \Big ] \\&\quad = {\left\{ \begin{array}{ll} {\mathscr {M}}_{\text {symm}}(p_1,p_2) &{} \text { for } \varGamma _i=\gamma ^\alpha , \, \sigma ^{\alpha \beta }, \, \sigma ^{\alpha \beta } \,\gamma ^5,\\ {\mathscr {M}}_{\text {anti-symm}}(p_1,p_2) &{} \text { for } \varGamma _i= {\textbf{1}}, \, \gamma ^5, \, \gamma ^\alpha \, \gamma ^5. \end{array}\right. } \end{aligned}$$

Therefore, the NP effects could lead to observable differences between the Dirac and Majorana neutrino interactions, even after integrating over the unknown neutrino and antineutrino 4-momenta, if the symmetric and anti-symmetric parts of the amplitude do not cancel either due to some imposed symmetry (like in the case of vector \(V-A\) interactions in the SM) or due to some accidental fine-tuning and numerical cancellation between different types of couplings. In Sect. 4 we illustrate the general considerations presented above discussing a chosen physical process of the form \(X \rightarrow Y \, \nu \, {\overline{\nu }}\).

4 Example 2-body and 3-body neutral current decays with \(\varvec{\nu \,{\overline{\nu }}}\) in the final state

To illustrate the results arising from the general formalism presented in previous sections, we study here in detail two example processes, in both cases assuming the non-vanishing NP effects, namely the 2-body decay of \(Z^0\) boson, \(Z \rightarrow \nu _\ell \, {\overline{\nu }}_\ell \) and the 3-body decay of an initial pseudo-scalar meson to a final pseudo-scalar meson and \(\nu \, {\overline{\nu }}\). A recent model-independent analysis of 3-body leptonic decays, but in context of heavy sterile neutrinos and facilitated by charged current interaction, can be found in Ref. [40].

4.1 2-Body decay of \(\varvec{Z \rightarrow \nu _\ell \, {\overline{\nu }}_\ell }\)

Considering Lorentz invariance alone, it is possible to write down the most general \(Z(p) \rightarrow \nu (p_1) \, {\overline{\nu }}(p_2)\) decay amplitude as follows (for Dirac neutrinos),

$$\begin{aligned} {\mathscr {M}}^D&= -i\,\epsilon _\alpha (p) \, {\overline{u}}(p_1) \, \Big [ \Big ( g_S^+ + g_P^+ \, \gamma ^5 \Big ) \, p^\alpha + \Big ( g_S^- + g_P^- \, \gamma ^5 \Big ) \, q^\alpha \nonumber \\&\quad + \gamma ^\alpha \Big ( g_V + g_A \, \gamma ^5 \Big ) + \sigma ^{\alpha \beta } \Big (g_{T_{\text {md}}}^+ + g_{T_{\text {ed}}}^+ \gamma ^5 \Big ) \, p_\beta \nonumber \\&\quad + \sigma ^{\alpha \beta } \Big (g_{T_{\text {md}}}^- + g_{T_{\text {ed}}}^- \gamma ^5 \Big ) \, q_\beta \Big ]\, v(p_2), \end{aligned}$$
(14)

where \(\epsilon _\alpha (p)\) denotes the polarisation 4-vector of Z, \(q=p_1 - p_2\), \(p=p_1 + p_2\), \(g_X^{(\pm )}\) with \(X = S, P, V, A, T_{\text {md}}, T_{\text {ed}}\) denote the various possible coupling constants corresponding to scalar, pseudo-scalar, vector, axial-vector, tensor (magnetic dipole) and axial-tensor (electric dipole) kind of interactions. In the SM, \(g_S^\pm = g_P^\pm = g_{T_{\text {md}}}^\pm = g_{T_{\text {ed}}}^\pm = 0\) and \(g_V = - g_A = \dfrac{g_Z}{4}\) where \(g_Z = e/(\sin \theta _W\,\cos \theta _W)\) with \(\theta _W\) being the weak mixing angle and e being the electric charge of positron. In presence of NP all the above coupling constants could, in principle, differ from the SM values. Noting that all terms proportional to \(p^\alpha \) vanish since \(p^\alpha \, \epsilon _\alpha (p) = 0\), and utilising Gordon identities as well as neglecting terms proportional to \(m_\nu \), the tensorial components get eliminated to yield,

$$\begin{aligned} {\mathscr {M}}^D= & {} -i\,\epsilon _\alpha (p) \, {\overline{u}}(p_1) \, \Big [ \Big ( g_S + g_P \, \gamma ^5 \Big ) \, q^\alpha \nonumber \\ {}{} & {} + \gamma ^\alpha \Big ( g_V + g_A \, \gamma ^5 \Big ) \Big ]\, v(p_2), \end{aligned}$$
(15)

where \(g_S = g_S^- + i\,g_{T_{\text {md}}}^+\) and \(g_P = g_P^- + i\,g_{T_{\text {ed}}}^+\). Further, if we consider CP and CPT conservation in the decay \(Z \rightarrow \nu _\ell \, {\overline{\nu }}_\ell \), \(g_S =0\) and \(g_P=0\) ought to be satisfied respectively. This implies, only vector and axial-vector couplings are allowed. The most general expression for the decay amplitude for \(Z \rightarrow \nu _\ell \, {\overline{\nu }}_\ell \) with Dirac neutrinos is given by,

$$\begin{aligned} {\mathscr {M}}^D= & {} {\mathscr {M}}(p_1, p_2)\nonumber \\ {}= & {} - \frac{i\, g_Z}{2} \upepsilon _\alpha (p) \, \left[ {\overline{u}}(p_1) \, \gamma ^\alpha \left( C_V^\ell -C_A^\ell \, \gamma ^5\right) \, v(p_2)\right] ,\quad \end{aligned}$$
(16)

where, to keep our analysis more general, we have elevated the parameters \(g_V\) and \(g_A\) to possibly be different for every the lepton family by adding the superscript \(\ell =e,\mu ,\tau \) such that the vector and axial-vector couplings are expressed as \(g_V^\ell = \dfrac{g_Z}{2} C_V^\ell \) and \(g_A^\ell = - \dfrac{g_Z}{2} C_A^\ell \) respectively, with

$$\begin{aligned} C_{V, A}^\ell = \frac{1}{2} + \varepsilon _{V, A}^\ell , \end{aligned}$$
(17)

where \(\varepsilon _V^\ell \), \(\varepsilon _A^\ell \) parameterise the NP effects, vanishing in the SM case. The amplitude for Majorana case is given by

$$\begin{aligned} {\mathscr {M}}^M&= \dfrac{1}{\sqrt{2}} \Big ({\mathscr {M}}(p_1, p_2) - {\mathscr {M}}(p_2, p_1) \Big ) \nonumber \\&= \frac{i\, g_Z\, C_A^\ell }{\sqrt{2}} \upepsilon _\alpha (p) \, \left[ {\overline{u}}(p_1) \, \gamma ^\alpha \, \gamma ^5 \, v(p_2)\right] . \end{aligned}$$
(18)

It is clear that we can combine the direct and exchange amplitudes in this case and effectively redefine the vertex structure for \(Z\rightarrow \nu _\ell \, {\overline{\nu }}_\ell \) when Majorana neutrinos are considered.

Keeping neutrino mass dependent terms in the amplitude squares, we get different results for Dirac and Majorana neutrinos:

$$\begin{aligned} \left| {\mathscr {M}}^D\right| ^2&= \frac{g_Z^2}{3} \Bigg (\left( (C_V^\ell )^2 + (C_A^\ell )^2\right) \left( m_Z^2 - m_\nu ^2\right) \nonumber \\&\quad + 3 \left( (C_V^\ell )^2 - (C_A^\ell )^2\right) m_\nu ^2\Bigg ), \end{aligned}$$
(19)
$$\begin{aligned} \left| {\mathscr {M}}^M\right| ^2&= \frac{2\, g_Z^2\, (C_A^\ell )^2}{3} \left( m_Z^2 - 4\, m_\nu ^2\right) , \end{aligned}$$
(20)

such that

$$\begin{aligned} \left| {\mathscr {M}}^D\right| ^2 - \left| {\mathscr {M}}^M\right| ^2&= \frac{g_Z^2}{3} \bigg (\left( (C_V^\ell )^2-(C_A^\ell )^2\right) \left( m_Z^2+2\, m_\nu ^2\right) \nonumber \\&\quad + 6\, (C_A^\ell )^2\, m_\nu ^2\bigg ) \nonumber \\&= {\left\{ \begin{array}{ll} \dfrac{g_Z^2}{2} m_\nu ^2,&\left( \begin{array}{c} \text {for SM alone}\end{array} \right) \\ \dfrac{g_Z^2}{3} \left( \varepsilon _V^\ell -\varepsilon _A^\ell \right) m_Z^2, &{} \left( \begin{array}{c}\text {with NP but}\\ \text {neglecting}\ m_\nu \end{array}\right) \end{array}\right. } \end{aligned}$$
(21)

where we have kept only the leading order contributions of \(\varepsilon _{V, A}^\ell \) while considering NP effects. It is clear that the SM result is fully in agreement with “practical Dirac Majorana confusion theorem”. However, the difference between Dirac and Majorana neutrinos appears in context of NP contributions even when one neglects \(m_\nu \) dependent terms (unless, of course, \(\varepsilon _V^\ell = \varepsilon _A^\ell \) in which case the additional NP contributions effectively rescale the SM allowed \(V-A\) coupling). Possible example of NP effects in this Z boson decay could arise from kinetic mixing of Z with the neutral gauge bosons from extra gauge groups like additional U(1) or \(SU(2)_R\). In this paper we are not concerned with any specific model of NP to keep our results very general.

Fig. 2
figure 2

Constraints on the NP parameters \(\varepsilon _V\) and \(\varepsilon _A\) from the experimentally measured number of light neutrino species \(N_\nu \). The point denoted by ‘\(+\)’ corresponds to the SM values \(\varepsilon _V = 0 = \varepsilon _A\). In the plots here \(\varepsilon _{V,A}^2\) contributions have not been neglected

Full size image

The decay width of Z boson into invisible final states is well measured and can constrain NP contributions to neutrino couplings. Neglecting small neutrino masses, the decay rates for \(Z\rightarrow \nu _\ell \, {\overline{\nu }}_\ell \) for Dirac and Majorana neutrino possibilities are given by

$$\begin{aligned} \varGamma ^D (Z\rightarrow \nu _\ell \, {\overline{\nu }}_\ell )&= \varGamma _Z^{0} \left( 1 + 2\, \varepsilon _V^\ell + 2\, \varepsilon _A^\ell \right) , \end{aligned}$$
(22a)
$$\begin{aligned} \varGamma ^M (Z\rightarrow \nu _\ell \, {\overline{\nu }}_\ell )&= \varGamma _Z^{0} \left( 1 + 4\, \varepsilon _A^\ell \right) , \end{aligned}$$
(22b)

where \(\varGamma _Z^{0}\) denotes the SM decay rate for \(m_\nu =0\),

$$\begin{aligned} \varGamma _Z^{0} = \dfrac{G_F\, m_Z^3}{12\sqrt{2}\, \pi }, \end{aligned}$$
(23)

with \(G_F\) being the Fermi constant. Thus, in presence of NP effects and neglecting \(m_\nu \) dependent terms, the total invisible width of Z boson is given by,

$$\begin{aligned} \varGamma _{Z, \text {inv}} = {\left\{ \begin{array}{ll} \displaystyle \varGamma _Z^0 \, \Bigg (3 + 2 \sum _{\ell =e, \mu , \tau } \left( \varepsilon _V^\ell + \varepsilon _A^\ell \right) \Bigg ), &{} \left( \begin{array}{c}\text {for Dirac}\\ \text {neutrinos}\end{array} \right) \\ \displaystyle \varGamma _Z^0 \, \Bigg (3 + 4 \sum _{\ell =e, \mu , \tau } \varepsilon _A^\ell \Bigg ). &{} \left( \begin{array}{c}\text {for Majorana}\\ \text {neutrinos}\end{array} \right) \end{array}\right. }\nonumber \\ \end{aligned}$$
(24)

The ratio \(\varGamma _{Z, \text {inv}}/\varGamma _Z^0\) is a measure of the number of light neutrino species \(N_\nu \) and its experimental estimate is [5, 41],

$$\begin{aligned} N_\nu = \varGamma _{Z, \text {inv}}/\varGamma _Z^0 = 2.9963 \pm 0.0074. \end{aligned}$$
(25)

Using Eq. (24) in Eq. (25) we get the following constraints on the NP parameters,

$$\begin{aligned}&\sum _{\ell =e, \mu , \tau } \left( \varepsilon _V^\ell + \varepsilon _A^\ell \right) = -0.0018 \pm 0.0037 ,\quad \left( \begin{array}{c}\text {for Dirac}\\ \text {neutrinos}\end{array} \right) \end{aligned}$$
(26a)
$$\begin{aligned}&\sum _{\ell =e, \mu , \tau } \varepsilon _A^\ell = -0.0009 \pm 0.0018,\quad \left( \begin{array}{c}\text {for Majorana}\\ \text {neutrinos}\end{array} \right) \end{aligned}$$
(26b)

which are perfectly consistent with zero, however obviously different for Majorana and Dirac neutrinos. In the latter case, it is in principle even possible that deviations from the SM \(V-A\)-type couplings are individually larger but cancel between the vector and axial contributions.

Assuming that the NP contributions do not depend on the flavor of the neutrino, i.e. \(\varepsilon _{V,A}^e = \varepsilon _{V,A}^\mu = \varepsilon _{V,A}^\tau \equiv \varepsilon _{V,A}\) (say), and keeping the contributions from \(\varepsilon _{V,A}^2\) terms as well, we find that \(\varepsilon _V\) and \(\varepsilon _A\) get constrained by Eq. (25) as shown in Fig. 2. It is clear from Fig. 2 that the SM values \(\varepsilon _V = 0 = \varepsilon _A\) (corresponding to the point indicated by ‘\(+\)’) are within the allowed \(1\sigma \) region.

The results above can be interpreted in two-fold way. First, one can assume some specific NP model that affects the weak neutral current processes involving neutrinos, with well defined Lagrangian and field interactions. In this case nature of the neutrinos follows from the model construction. As always, to predict the phenomenological consequences of the given SM extension, one needs first to estimate the experimentally allowed size of NP couplings. An example above illustrates how such bounds differ for a particular decay depending on the type of neutrinos in a model. Obviously, the realistic NP models usually contain a number of new parameters and constraining them requires combining bounds from various measurements, often not only from the neutrino sector. The invisible Z boson decay discussed in this section, as a quantity known to a good accuracy, should always be one of observables contributing to such more general analysis and fitting NP model to data.

Second, it may be not known in advance what is the nature of the neutrino fields, apart from the assumption that they can have interactions going beyond the SM-like \(V-A\) couplings. In this case single measurement like the invisible Z boson decay cannot help distinguish the type of neutrino fields if some deviations from the SM predictions are observed (or even if neutrinos are responsible for such deviations). However, again when more measurements are combined (especially for the final states with more than 2 particles, which allow to construct additional observables like shown in the next section) it may be possible to constrain individual couplings and, in an optimistic scenario, understand the nature of “invisible” sector of the theory. Again, reconstructing the model couplings must take into account the fact that it depends on the nature of neutrino fields.

Finally, apart from specific decay discussed in this section, the general formalism presented before can be applied to other 2-body neutrino decays which can be considered in searches for various channels of decays of BSM particles, e.g. new neutral heavy gauge bosons or scalars (here possible example are the sneutrino decays in R-parity violating MSSM).

4.2 3-Body decay of an initial pseudo-scalar meson to a final pseudo-scalar meson and \(\varvec{\nu \,{\overline{\nu }}}\)

As a second example let us consider the general decay \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu _\ell \, {\overline{\nu }}_\ell \) where \({\mathcal {P}}_i\) is the parent pseudo-scalar meson with mass \(M_i\) and \({\mathcal {P}}_f\) is the daughter pseudo-scalar meson with mass \(M_f\). For example, \({\mathcal {P}}_i\) could be B or K meson, then \({\mathcal {P}}_f\) would be either K or \(\pi \) meson respectively.

Considering Lorentz invariance, the effective Lagrangian for \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu _\ell \, {\overline{\nu }}_\ell \) decay can be written as follows,

$$\begin{aligned} {\mathscr {L}}&= J_{SL}^\ell \Big ({\overline{\psi }}_\nu \,P_L\, \psi _{{\overline{\nu }}} \Big ) + J_{SR}^\ell \Big ({\overline{\psi }}_\nu \, P_R \, \psi _{{\overline{\nu }}} \Big ) \nonumber \\&\quad + \left( J_{VL}^\ell \right) _{\alpha } \Big ({\overline{\psi }}_\nu \, \gamma ^\alpha P_L \, \psi _{{\overline{\nu }}} \Big ) + \left( J_{VR}^\ell \right) _\alpha \Big ({\overline{\psi }}_\nu \, \gamma ^\alpha \, P_R \, \psi _{{\overline{\nu }}} \Big ) \nonumber \\&\quad + \left( J_{TL}^\ell \right) _{\alpha \beta } \Big ({\overline{\psi }}_\nu \,P_L\, \sigma ^{\alpha \beta } \, \psi _{{\overline{\nu }}} \Big ) \nonumber \\ {}&\quad + \left( J_{TR}^\ell \right) _{\alpha \beta } \Big ({\overline{\psi }}_\nu \, \sigma ^{\alpha \beta }\, P_R \, \psi _{{\overline{\nu }}} \Big ) + \text {h.c.}, \end{aligned}$$
(27)

where \(\psi _j\) denotes the fermionic field of \(j=\nu _\ell , {\overline{\nu }}_\ell \) and \(J_{SL}^\ell \), \(J_{SR}^\ell \), \(\left( J_{VL}^\ell \right) _\alpha \), \(\left( J_{VR}^\ell \right) _\alpha \), \(\left( J_{TL}^\ell \right) _{\alpha \beta }\) \(\left( J_{TR}^\ell \right) _{\alpha \beta }\) denote the different hadronic currents describing the quark level transitions from \({\mathcal {P}}_i\) to \({\mathcal {P}}_f\) meson. We have used the superscript \(\ell \) to accommodate any possible differences in the effective Lagrangian when different neutrino flavors are considered. The subscripts S, V and T in the hadronic currents denote the fact that the associated external leptonic currents are scalar, vector and tensor type, respectively, and in addition they carry also the chirality index. In the SM, at tree-level, only the \(V-A\) interaction is present, while the scalar and tensor interactions involving left-handed neutrinos might be generated from quantum loops and are hence expected to be suppressed. Additionally, in the SM the right-handed neutrinos do not take part in any interaction. In this analysis we do not consider any specific NP model to keep our results general, thus in Eq. (27) we have included all possible interactions disregarding the helicity of neutrino.

The decay amplitude for \({\mathcal {P}}_i (p_i) \rightarrow {\mathcal {P}}_f (p_f) \, \nu _\ell (p_1) \, {\overline{\nu }}_\ell (p_2)\), considering the neutrino and the antineutrino to be Dirac fermions, is given by,

$$\begin{aligned} {\mathscr {M}}^D = {\mathscr {M}}(p_1, p_2)&= {\overline{u}}(p_1) \Big [ F_{SL}^\ell \, P_L\, + F_{SR}^\ell \, P_R \nonumber \\&\quad + \left( F_{VL}^{\ell +} \, p_\alpha + F_{VL}^{\ell -} q_\alpha \right) \, \gamma ^\alpha \, P_L \nonumber \\&\quad + \left( F_{VR}^{\ell +} \, p_\alpha + F_{VR}^{\ell -} q_\alpha \right) \, \gamma ^\alpha \, P_R \, \nonumber \\&\quad + F_{TL}^\ell \, p_\alpha \, q_\beta \, \sigma ^{\alpha \beta } \,P_L \nonumber \\&\quad + F_{TR}^\ell \, p_\alpha \, q_\beta \, \sigma ^{\alpha \beta } \, P_R \Big ] v(p_2), \end{aligned}$$
(28)

where \(p=p_i + p_f\), \(q=p_i-p_f = p_1 + p_2\) and the various form factors \(F_{SX}^\ell \), \(F_{VX}^{\ell \pm }\), \(F_{TX}^\ell \) are defined as follows,

$$\begin{aligned} \left\langle {\mathcal {P}}_f \Big |J_{SX}^\ell \Big |{\mathcal {P}}_i \right\rangle&= F_{SX}^\ell , \end{aligned}$$
(29a)
$$\begin{aligned} \left\langle {\mathcal {P}}_f \Big |\left( J_{VX}^\ell \right) _\alpha \Big |{\mathcal {P}}_i \right\rangle&= F_{VX}^{\ell +} \, p_\alpha + F_{VX}^{\ell -} \, q_\alpha , \end{aligned}$$
(29b)
$$\begin{aligned} \left\langle {\mathcal {P}}_f \Big |\left( J_{TX}^\ell \right) _{\alpha \beta } \Big |{\mathcal {P}}_i \right\rangle&= F_{TX}^\ell \, p_\alpha \, q_\beta , \end{aligned}$$
(29c)

where by X we denote the chirality index, \(X=L,R\).

The form factors defined in Eq. (29) are complex and, in general, functions of the square of the momentum transferred, i.e. \(q^2\equiv s\). Here the form factors also include the CKM matrix elements and their mass dimensions are different since the explicit dependencies on quark/meson masses are hidden for simplicity of expressions.

The decay amplitude for Majorana neutrinos is obtained by anti-symmetrization with respect to exchange of \(p_1\) and \(p_2\) and is given by,

$$\begin{aligned} {\mathscr {M}}^M&= \frac{1}{\sqrt{2}} \Big ({\mathscr {M}}(p_1, p_2) - {\mathscr {M}}(p_2, p_1)\Big ) \nonumber \\&= \sqrt{2} \, {\overline{u}}(p_1) \,\Bigg [ F_{SL}^\ell \,P_L + F_{SR}^\ell \, P_R \nonumber \\&\quad + \bigg (\frac{F_{VR}^{\ell +} - F_{VL}^{\ell +}}{2} \, p_\alpha + \frac{F_{VR}^{\ell -} - F_{VL}^{\ell -}}{2} q_\alpha \bigg ) \, \gamma ^\alpha \, \gamma ^5 \Bigg ]\, v(p_2). \end{aligned}$$
(30)

Once again we find that the direct and exchange amplitudes in this case can be combined to effectively redefine the vertex factors for each possible interaction. Notably, for Majorana neutrinos no vector or tensor neutral currents are possible.

Since the scalar products that appear in the amplitude squares are Lorentz invariant, they can be evaluated in any chosen frame of reference. It is well known that any 3-body decay can be fully described by two independent parameters. For simplicity we choose them as \(s = \big (p_1 + p_2\big )^2\) and the angle \(\theta \) between the 3-momenta of the neutrino and the final state meson in the center-of-momentum frame of the \(\nu \, {\overline{\nu }}\) pair, as shown in Fig. 3.

Fig. 3
figure 3

Kinematic configuration of the decay \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu \, {\overline{\nu }}\) in center-of-momentum frame of the \(\nu \, {\overline{\nu }}\) pair

Full size image

The amplitude squares for both Dirac and Majorana cases can be written in the following form in terms of the variables s and \(\cos \theta \),

$$\begin{aligned} \left|{\mathscr {M}}^{D/M} \right|^2 = C_0^{D/M} + C_1^{D/M} \, \cos \theta + C_2^{D/M} \, \cos ^2\theta , \end{aligned}$$
(31)

where, neglecting the small terms proportional to neutrino mass \(m_\nu \) we have,

$$\begin{aligned} C_0^D&= s\,\Big ( \left| F_{SL}^\ell \right| ^2 + \left| F_{SR}^\ell \right| ^2 \Big ) + \lambda \, \Big (\left| F_{VL}^{\ell +} \right| ^2 + \left| F_{VR}^{\ell +} \right| ^2\Big ) , \end{aligned}$$
(32a)
$$\begin{aligned} C_1^D&= 2\, s \,\sqrt{\lambda }\, \left( \text {Im}\left( F_{SL}^\ell F_{TL}^{\ell *}\right) + \text {Im}\left( F_{SR}^\ell F_{TR}^{\ell *} \right) \right) , \end{aligned}$$
(32b)
$$\begin{aligned} C_2^D&= - \lambda \, \Big ( \left| F_{VL}^{\ell +} \right| ^2 + \left| F_{VR}^{\ell +} \right| ^2 - s \, \left( \left| F_{TL}^\ell \right| ^2 + \left| F_{TR}^\ell \right| ^2 \right) \Big ), \end{aligned}$$
(32c)
$$\begin{aligned} C_0^M&= 2 \,s\, \Big (\left| F_{SL}^\ell \right| ^2 + \left| F_{SR}^\ell \right| ^2\Big ) + \lambda \, \left| F_{VL}^{\ell +} - F_{VR}^{\ell +} \right| ^2 ,\end{aligned}$$
(32d)
$$\begin{aligned} C_1^M&= 0, \end{aligned}$$
(32e)
$$\begin{aligned} C_2^M&= - \lambda \, \left| F_{VL}^{\ell +} - F_{VR}^{\ell +} \right| ^2 , \end{aligned}$$
(32f)

with

$$\begin{aligned} \lambda = M_i^4 + M_f^4 + s^2 - 2 \left( M_i^2 \, M_f^2 + s \, M_i^2 + s \, M_f^2\right) . \end{aligned}$$
(33)

The Dirac amplitude square may contain a term odd in \(\cos \theta \), while it is always absent in Majorana case. This, in principle, could be probed if the angle \(\theta \) could be measured experimentally. However, at present the angle \(\theta \) can not be observed as both the neutrino and antineutrino remain undetected in the detector near their point of production. In addition, such a term also requires simultaneously existence of both tensor and scalar NP interactions (in addition, with complex couplings) and is likely to be small.

The differential decay rate in the rest frame of the parent meson \({\mathcal {P}}_i\) after integration over the unobservable \(\cos \theta \) is given by,

$$\begin{aligned} \frac{\textrm{d}\varGamma ^{D/M}}{\textrm{d}s} = \frac{1}{(2\, \pi )^3} \frac{b}{16\, M_i^3} \left( C_0^{D/M} + \frac{1}{3} C_2^{D/M}\right) , \end{aligned}$$
(34)

where

$$\begin{aligned} b = \frac{\sqrt{\lambda }}{2} \, \sqrt{1- \frac{4\,m_\nu ^2}{s}}. \end{aligned}$$
(35)

In terms of the form factors, we get

$$\begin{aligned} \frac{\textrm{d}\varGamma ^{D}}{\textrm{d}s}&= \frac{1}{(2\, \pi )^3} \frac{\lambda \,b}{24\, M_i^3} \Bigg ( \left| F_{VL}^{\ell +} \right| ^2 + \left| F_{VR}^{\ell +} \right| ^2 \nonumber \\&\quad + \frac{3s}{2\lambda } \, \left( \left| F_{SL}^\ell \right| ^2 + \left| F_{SR}^\ell \right| ^2 \right) \nonumber \\&\quad + \frac{s}{2} \, \left( \left| F_{TL}^\ell \right| ^2 + \left| F_{TR}^\ell \right| ^2 \right) \Bigg )\, , \end{aligned}$$
(36a)
$$\begin{aligned} \frac{\textrm{d}\varGamma ^{M}}{\textrm{d}s}&= \frac{1}{(2\, \pi )^3} \frac{\lambda \,b}{24\, M_i^3} \Bigg ( \left| F_{VL}^{\ell +} - F_{VR}^{\ell +} \right| ^2 + \frac{3s}{\lambda } \, \left( \left| F_{SL}^\ell \right| ^2 + \left| F_{SR}^\ell \right| ^2 \right) \Bigg ). \end{aligned}$$
(36b)

From Eqs. (36a) and (36b) it is clear that the presence of additional forms of NP interactions on top of the usual SM allowed \(V-A\) interaction can lead to observable differences between the Dirac and Majorana neutrinos.

To show the impact of presence of such additional non-standard interactions we consider a simple effective parametrization of NP effects (using again the chirality index \(X=L,R\)),

$$\begin{aligned} F_{SX}^\ell&= M_i \, F_{\text {SM}} \, \varepsilon _{SX}^\ell , \end{aligned}$$
(37a)
$$\begin{aligned} F_{VL}^{\ell +}&= F_{\text {SM}} \,\Big (1+\varepsilon _{VL}^\ell \Big ), \end{aligned}$$
(37b)
$$\begin{aligned} F_{VR}^{\ell +}&= F_{\text {SM}}\, \varepsilon _{VR}^\ell , \end{aligned}$$
(37c)
$$\begin{aligned} F_{TX}^\ell&= \frac{1}{M_i} \, F_{\text {SM}} \, \varepsilon _{TX}^\ell , \end{aligned}$$
(37d)

where \(\varepsilon _{QX}^\ell \), \(Q=S,V,T\), denote the relative size of the NP contributions with respect to the SM contributionFootnote 2 and are assumed to be small, \( \left| \varepsilon _{QX}^\ell \right| \ll 1\). We include \(M_i\) factors in Eq. (37) to accommodate for the difference in mass dimensions of various form factors so as to keep the NP parameters \(\varepsilon _{QX}^\ell \) dimensionless. Substituting Eq. (37) in Eq. (36) we obtain,

$$\begin{aligned} \frac{\textrm{d}\varGamma ^{D}}{\textrm{d}s}&= \frac{\textrm{d}\varGamma ^{\text {SM}}}{\textrm{d}s} \Bigg (1+ 2\, \text {Re}\, \varepsilon _{VL}^\ell + \left| \varepsilon _{VL}^\ell \right| ^2 + \left| \varepsilon _{VR}^\ell \right| ^2 \nonumber \\&\quad + \frac{3 s M_i^2}{2\lambda } \, \left( \left| \varepsilon _{SL}^\ell \right| ^2 + \left| \varepsilon _{SR}^\ell \right| ^2 \right) \nonumber \\&\quad + \frac{s}{2 M_i^2} \, \left( \left| \varepsilon _{TL}^\ell \right| ^2 + \left| \varepsilon _{TR}^\ell \right| ^2 \right) \Bigg ) \nonumber \\&\approx \frac{\textrm{d}\varGamma ^{\text {SM}}}{\textrm{d}s} \Bigg (1+ 2\, \text {Re}\, \varepsilon _{VL}^\ell \Bigg )\, , \end{aligned}$$
(38a)
$$\begin{aligned} \frac{\textrm{d}\varGamma ^{M}}{\textrm{d}s}&= \frac{\textrm{d}\varGamma ^{\text {SM}}}{\textrm{d}s} \Bigg (1 + 2\, \text {Re}\, \varepsilon _{VL}^\ell - 2\, \text {Re}\, \varepsilon _{VR}^\ell + \left| \varepsilon _{VL}^\ell - \varepsilon _{VR}^\ell \right| ^2 \nonumber \\&\quad + \frac{3 s M_i^2}{\lambda } \, \left( \left| \varepsilon _{SL}^\ell \right| ^2 + \left| \varepsilon _{SR}^\ell \right| ^2 \right) \Bigg ) \nonumber \\&\approx \frac{\textrm{d}\varGamma ^{\text {SM}}}{\textrm{d}s} \Bigg (1 + 2\, \text {Re}\, \varepsilon _{VL}^\ell - 2\, \text {Re}\, \varepsilon _{VR}^\ell \Bigg ), \end{aligned}$$
(38b)

where we defined the SM decay rate as

$$\begin{aligned} \frac{\textrm{d}\varGamma ^{\text {SM}}}{\textrm{d}s} = \frac{1}{(2\, \pi )^3} \frac{\lambda \, b\,\left| F_{\text {SM}} \right| ^2}{24\, M_i^3}\,. \end{aligned}$$
(39)

Therefore, the difference between Dirac and Majorana cases, to the leading order in the NP parameters, is given by,

$$\begin{aligned} \frac{\textrm{d}\varGamma ^D}{\textrm{d}s} - \frac{\textrm{d}\varGamma ^M}{\textrm{d}s} \approx 2 \, \frac{\textrm{d}\varGamma ^{\text {SM}}}{\textrm{d}s} \, \text {Re}\,\varepsilon _{VR}^\ell , \end{aligned}$$
(40)

which implies that non-zero difference can, in principle, arise only if \(\text {Re}\,\varepsilon _{VR}^\ell \ne 0\).

Similar to the case of \(Z\rightarrow \nu \,{\overline{\nu }}\) decay, only the NP contributions to the vector currents contribute to both decay rates in the lowest linear order. In addition, presence of only \(\varepsilon _{VL}\), i.e. change of normalisation of the SM \(V-A\) interaction, modifies the decay into Dirac neutrinos, while both \(\varepsilon _{VL}\) and \(\varepsilon _{VR}\) affect the expression for the Majorana neutrino decay. Scalar and tensor interactions lead to higher order corrections, quadratic in NP effects, and are likely to be small and difficult to observe taking into experimental accuracy of corresponding measurement, both current and in foreseeable future. Moreover, if only constant vector or scalar NP contributions are present, the shape of the decay spectrum \(\textrm{d}\varGamma ^{D/M}/\textrm{d}s\) is identical for Majorana and Dirac neutrinos and, therefore, when used alone it can not distinguish between the two possibilities (unless in some specific NP model both \(\varepsilon _{VL}\) and \(\varepsilon _{VR}\) are known and thus the difference in normalisation of the spectrum is also known). However, it is interesting to note that if tensor NP contributions are substantial enough to be observed, it would point out to the Dirac nature of neutrino, since tensor interactions do not affect the s distribution for Majorana neutrinos at all. We elaborate on these aspects in Sect. 5 where we study the s distributions while considering each of the possible interactions one at a time.

Fig. 4
figure 4

Constraints on the NP parameters \(\varepsilon _{VL}\) and \(\varepsilon _{VR}\) (assuming them to be real and including quadratic term contributions) using experimentally measured branching ratio of \(K^+ \rightarrow \pi ^+ \, \nu \, {\overline{\nu }}\) as given in Eq. (41). The point denoted by ‘\(+\)’ corresponds to the SM values \(\varepsilon _{VL} = 0 = \varepsilon _{VR}\)

Full size image

As an example of experimental constraints which can be imposed on the NP neutrino couplings by the three-body decays, let us consider the decay \(K^+ \rightarrow \pi ^+ \, \nu \, {\overline{\nu }}\) whose branching ratio as reported by NA62 Collaboration [42] is,

$$\begin{aligned} \text {Br}^{\text {Exp}}\left( K^+ \rightarrow \pi ^+ \, \nu \, {\overline{\nu }}\right) = \left( \left. 10.6^{+4.0}_{-3.4}\right| _{\text {stat}} \pm 0.9_{\text {syst}}\right) \times 10^{-11}, \end{aligned}$$
(41)

which is consistent with the SM prediction [43],

$$\begin{aligned} \text {Br}^{\text {SM}}\left( K^+ \rightarrow \pi ^+ \, \nu \, {\overline{\nu }}\right) = \left( 9.11 \pm 0.72\right) \times 10^{-11}. \end{aligned}$$
(42)

Neglecting the terms quadratic in NP contributions and assuming that terms that are linear in the parameters \(\varepsilon _{QX}\) are approximately constant with respect to \(q^2\equiv s\) (which could affect the integration of \(\textrm{d}\varGamma /\textrm{d}s\) to the total branching ratio), we get the following constraints from the above mentioned experimental measurement,

$$\begin{aligned}&\sum _{\ell =e, \mu , \tau } \text {Re}\,\varepsilon _{VL}^\ell = 0.082 \pm 0.214, \quad \left( \begin{array}{c}\text {for Dirac}\\ \text {neutrinos}\end{array}\right) \end{aligned}$$
(43a)
$$\begin{aligned}&\sum _{\ell =e, \mu , \tau } \left( \text {Re}\,\varepsilon _{VL}^\ell - \text {Re}\,\varepsilon _{VR}^\ell \right) \nonumber \\ {}&\quad = 0.082 \pm 0.214, \quad \left( \begin{array}{c}\text {for Majorana}\\ \text {neutrinos}\end{array}\right) \end{aligned}$$
(43b)

which are consistent with zero.

If we assume that only vector NP contributions are dominant and \(\varepsilon _{VL}^\ell \), \(\varepsilon _{VR}^\ell \) are both real and identical for all neutrino flavors, i.e. \(\varepsilon _{VX}^e = \varepsilon _{VX}^\mu = \varepsilon _{VX}^\tau \equiv \varepsilon _{VX}\) (with \(X=L,R\)), then the experimental measurement of Eq. (41) constrains the NP parameters \(\varepsilon _{VL}\) and \(\varepsilon _{VR}\) as shown in Fig. 4 (there we have kept the quadratic terms in \(\varepsilon _{VL}\), \(\varepsilon _{VR}\) as well, assuming that in principle they can be large compared to the SM terms but cancel out due to some kind of accidental or symmetry-related fine-tuning). It is clear from Fig. 4 that if we allow for some amount of fine tuning and cancellations, large NP contributions are still possible in both Dirac and Majorana neutrino possibilities. One can also note that the contributions from scalar and tensor form factors are always positive, so even if they are also non-negligible, they can only tighten the bounds on \(\varepsilon _{VL}\) and \(\varepsilon _{VR}\) plotted in Fig. 4.

With more precise experimental decays data in future, the estimate of NP contribution will also improve. Another interesting decay in the same category as \(K^+\rightarrow \pi ^+\,\nu \,{\overline{\nu }}\) is the decay \(B^+ \rightarrow K^+ \, \nu \, {\overline{\nu }}\) which is being pursued at Belle II and awaiting its branching ratio measurement. The SM prediction for the branching ratio for this decay [44] is,

$$\begin{aligned} \text {Br}^{\text {SM}} \left( B^+ \rightarrow K^+ \, \nu \, {\overline{\nu }} \right) = \left( 4.6 \pm 0.5 \right) \times 10^{-6}, \end{aligned}$$
(44)

while the current experimental upper limit by Belle II [45] is \(4.1 \times 10^{-5}\) at \(90\%\) confidence level.

Again, as already discussed in the previous section, constraining individual parameters of the NP model (or in the best case scenario even determining the nature of invisible sector of the theory, including distinguishing the Dirac and Majorana character of neutrinos) requires combining bounds from many measurements, in addition to the example discussed above.

5 Distinguishing Dirac and Majorana neutrinos via NP effects

Although experimentally measured branching ratios of the two- and three-body decays can only constrain the probable NP couplings of neutrinos, from Eq. (36) and the non-approximated part of Eq. (38) it seems that by analysing observed s (missing mass-square) distribution in the three-body decays one can also distinguish between Dirac and Majorana neutrino possibilities, provided it can be measured with sufficient accuracy. This could be very challenging and may not possible in the near future as meson decays into neutrino pair are rare and current experimental statistics are not sufficient to probe the energy spectra of the visible final state particles (equivalent to the s distribution). However, it might eventually be possible when more data get accumulated. Alternatively, similar analysis could be applied to collision experiments with the neutrino antineutrino pair in the final state, where statistics are much higher but one has to deal with significant backgrounds and related uncertainties.

As mentioned in the Introduction, one should note that modification of the spectra of the visible particle in the final state of the decay can also be caused if some invisible particles other than neutrinos are present in the final state. For distinguishing such NP possibilities, i.e. the presence of other invisible particles, we encourage the reader to follow other existing proposals such as the ones discussed in Refs. [38, 46]. Specifically, in context of 3-body semi-hadronic decays of mesons such as \(K \rightarrow \pi \, f \, {\overline{f}}\) where f could be fermionic dark matter, heavy neutral lepton or some long-lived fermion in addition to the usual active neutrino possibility, the methods proposed in Ref. [46] using angular distributions and Dalitz plot asymmetries could be useful in distinguishing these NP possibilities, without worrying about hadronic uncertainties. In the discussion below, we focus only on active neutrinos as the invisible particles and study the effect of NP interactions over the experimentally observable distribution.

Before proceeding further we note that in Eq. (38), \(\varepsilon _{SL}\) and \(\varepsilon _{SR}\) as well as \(\varepsilon _{TL}\) and \(\varepsilon _{TR}\) contribute in a similar manner. Therefore, we can discuss their impact in terms of two effective NP parameters \(\varepsilon _S\) and \(\varepsilon _T\), defined as follows,

$$\begin{aligned} \varepsilon _S&= \sqrt{|\varepsilon _{SL}|^2 +|\varepsilon _{SR}|^2}, \end{aligned}$$
(45a)
$$\begin{aligned} \varepsilon _T&= \sqrt{|\varepsilon _{TL}|^2 +|\varepsilon _{TR}|^2}. \end{aligned}$$
(45b)

Further we also consider, for simplicity, that all the NP parameters are real quantities and agnostic to the flavor of the neutrino in the final state. Below we shall first look at the generic trends in the various s distributions for various NP possibilities without considering any specific decay, and then we show how these generic features can affect the spectrum of a specific decay, on an example chosen to be \(B \rightarrow K \, \nu \, {\overline{\nu }}\).

5.1 The generic trends in \({\varvec{s}}\) (missing mass-square) distributions

Fig. 5
figure 5

Comparison of various generic NP contributions in \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu _\ell \, {\overline{\nu }}_\ell \) decay assuming that the mass of the final meson \({\mathcal {P}}_f\) can be neglected in comparison with the mass of the initial meson \({\mathcal {P}}_i\) and assuming constant form factors. As noted in the main text, all NP possibilities, except the left-handed vector NP contribution, can in principle distinguish between Dirac and Majorana neutrinos

Full size image

It is possible to infer the overall trends and order of magnitude effects of the NP contributions in the s or missing mass-square distributions, without considering any specific \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu \, {\overline{\nu }}\) decay, and by imposing the following simplifying assumptions in Eqs. (38a) and (38b):

  1. 1.

    the form factor \(F_{\text {SM}}\) has no s dependence,

  2. 2.

    neglect \({\mathcal {O}}\left( M_f^2/M_i^2\right) \) terms, and

  3. 3.

    neglect the \(m_\nu \) dependent term.

It is clear that except the last assumption, the other two need not even be good approximations in a real world example. We will relax these assumptions later while considering a specific decay. For the time being, we note that the above assumptions simplify the analysis of Eq. (38) which can now be cast into the discussion of the dimensionless and normalised distributions \(X^{D/M}\) defined as,

$$\begin{aligned} X^{D/M} = \frac{M_i^2}{\varGamma ^{\text {SM}}} \frac{\textrm{d}\varGamma ^{D/M}}{\textrm{d}s} = \frac{1}{\varGamma ^{\text {SM}}} \frac{\textrm{d}\varGamma ^{D/M}}{\textrm{d}{\tilde{s}}}, \end{aligned}$$
(46)

with,

$$\begin{aligned} \varGamma ^{\text {SM}} \approx \frac{1}{(2\, \pi )^3} \frac{M_i^5\,\left| F_{\text {SM}} \right| ^2}{192}\,. \end{aligned}$$
(47)

Using the simplifying assumptions, \(X^{D/M}\) can be rewritten in the following universal form,

$$\begin{aligned} X^{D}&= 4 \left( 1-{\tilde{s}}\right) \Bigg (\left( 1-{\tilde{s}}\right) ^2 \left( 1+ 2\, \varepsilon _{VL} + \varepsilon _{VL}^2 + \varepsilon _{VR}^2 + \frac{{\tilde{s}}}{2} \, \varepsilon _{T}^2 \right) \nonumber \\&\quad + \frac{3}{2} {\tilde{s}} \, \varepsilon _{S}^2 \Bigg ), \end{aligned}$$
(48a)
$$\begin{aligned} X^{M}&= 4 \left( 1-{\tilde{s}}\right) \Bigg (\left( 1-{\tilde{s}}\right) ^2 \left( 1 + 2\, \varepsilon _{VL} - 2\, \varepsilon _{VR} + \left( \varepsilon _{VL} - \varepsilon _{VR} \right) ^2 \right) \nonumber \\&\quad + 3\,{\tilde{s}} \, \varepsilon _{S}^2 \Bigg ) , \end{aligned}$$
(48b)

where \({\tilde{s}}=s/M_i^2\) is a dimensionless quantity which varies in the range [0, 1].

To illustrate the order of the expected magnitude of effects of NP neutrino couplings in the decay distributions, let us consider contribution from each of the NP parameters individually and one at a time. In Fig. 5 we plot, in presence of individual NP contributions, the distributions of the quantities \(X^{D/M}\). The main features of these distributions are as follows.

  1. 1.

    The effect of left-handed vector NP contributions is identical for both Dirac and Majorana neutrinos. Therefore, such a NP contribution would not lead to any difference between the Dirac and Majorana possibilities.

  2. 2.

    The right-handed vector NP contribution has opposite effects on Dirac and Majorana neutrinos. For Dirac case the differential decay rate is enhanced, while for Majorana case it is reduced. This opposite behaviour can clearly distinguish between the two possibilities.

  3. 3.

    As is also clear from Eqs. (48a) and (48b) (also Eqs. (38a) and (38b)), the effect of scalar NP contributions is enhanced in case of Majorana neutrinos by a factor 2 compared to that in case of Dirac neutrinos.

  4. 4.

    The tensor NP contributions do not contribute to the differential decay rate for Majorana neutrinos, while for Dirac neutrinos such rate gets enhanced. However, from Figs. 5 it is clear that such enhancement is not substantial when compared with effects from other NP possibilities. Therefore, tensor NP interaction would be harder to probe than the other NP possibilities.

5.2 Patterns in \({\varvec{s}}\) distributions of \(\varvec{B \rightarrow K \, \nu \, {\overline{\nu }}}\) decay in presence of NP

Fig. 6
figure 6

Comparison of various generic NP contributions in \(B \rightarrow K \, \nu _\ell \, {\overline{\nu }}_\ell \) decay. We use the scalar and tensor NP parameters \(\varepsilon _{S}\) and \(\varepsilon _{T}\) as defined in Eq. (45). \(\varGamma _{\text {SM}}\) denotes the partial decay rate of \(B \rightarrow K \, \nu _\ell \, {\overline{\nu }}_\ell \) in the SM. Also note the different overall normalisation factors used to show the distributions more clearly

Full size image

The discussion in Sect. 5.1 is applicable to all \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu _\ell \, {\overline{\nu }}_\ell \) decays but under the set of strong assumptions which may not be satisfied in real processes. In particular, the form factor \(F_{\text {SM}}\) is known to be in general s-dependent, as shown by the QCD lattice calculations. To study the effect of s dependence of \(F_{\text {SM}}\) on the distinguishing features of Dirac and Majorana neutrinos, let us analyse as an example the decay \(B \rightarrow K \nu {\overline{\nu }}\), relaxing the assumptions we had made in Sect. 5.1. In particular, we use here the functional form of \(B \rightarrow K\) form factors in the SM as given in Ref. [47] and consider the full form of Eqs. (38a) and (38b) (although we have continued to assume the unknown NP form factors \(\varepsilon _{X}^{B \rightarrow K}\) to be independent of s for simplicity).

Without the simplifying assumptions of Sect. 5.1, the detailed shape of s distributions changes as shown in Fig. 6, but they remain visibly different for Dirac and Majorana neutrino possibilities. We note the following salient features noticeable in Fig. 6.

  1. 1.

    The general trends noted before in context of Fig. 5 are also clearly seen in Fig. 6. The presence of NP contributions, except the left-handed vector NP contribution, do indeed lead to different and distinguishable distributions for Dirac and Majorana neutrinos.

  2. 2.

    The distributions in Fig. 6 are significantly altered comparing to the ones shown in Fig. 5 for larger values of s, in concordance with the enhancement of \(F_\mathrm {\text {SM}}\) with increasing s.

  3. 3.

    We note also that if \(m_\nu \) is not neglected, the distributions plotted in Fig. 6 actually drop to zero at the minimal value of \(s = 4\,m_\nu ^2\) (when the factor b in the numerator of Eq. (39) vanishes). However, as \(m_\nu \) is extremely tiny in comparison with other masses, such drop of distributions is very steep and so close to \(s=0\) that it is unlikely to be observed experimentally.

We would like to also emphasise the fact that in Figs. 5 and 6, the NP parameters are shown to be positive just for the purpose of illustration. In Eqs. (38) and (48), it is clear that the scalar and tensor NP contributions are insensitive to the sign of the NP parameter. For left-handed vector NP contribution alone, both Dirac and Majorana cases yield identical s distributions, where as for right-handed vector NP contribution alone, the Dirac and Majorana cases would still have opposite behaviour. Therefore, our conclusions regarding distinguishability of Dirac and Majorana nature from observed s distribution still holds true, in general, even for negative NP parameters.

Study of the energy spectra of the decays with neutrino antineutrino pair in the final state is very challenging as such decays have very low branching ratios. Notably also the example considered in this section, \(B \rightarrow K \nu {\overline{\nu }}\), has yet to be observed, although current experimental bound for it is already less than order of magnitude worse than the prediction of its SM decay rate [44, 45].

In order to estimate the accuracy of spectra measurement required to see the difference between the Dirac and Majorana neutrinos, let us define the quantity

$$\begin{aligned} \varDelta X = \frac{\left| X^D-X^M \right| }{X^D}. \end{aligned}$$
(49)

One can assume two scenarios. First, one can estimate that the maximal size of neutrino NP couplings to be of the order of limit given in Eq. (43), so \({{{\mathcal {O}}}}(0.1)\). The second scenario assume that there is a fine tuning between NP couplings and they may be large, \({{{\mathcal {O}}}}(1)\), as shown in Fig. 4. This is possible only for \(\varepsilon _{VL}\) and \(\varepsilon _{VR}\) as scalar and tensor contributions to total decay branching ratios are always positive and cannot cancel out with other terms. It turns out that the relative modification of spectrum shape, given by \(\varDelta X\), is of the similar order as the size of the (dimensionless) coupling \(\varepsilon _X\) themselves. This can be e.g. seen assuming that only \(\varepsilon _{VR}\) does not vanish. Then one has for every \(\tilde{s}\) the fixed ratio

$$\begin{aligned} \varDelta X = \frac{2\left| \varepsilon _{VR} \right| }{1+ \varepsilon _{VR}^2}. \end{aligned}$$
(50)

Given the rarity of considered decays and statistics necessary to obtain its energy spectrum, the first scenario \(\varepsilon _{VR} \sim \varDelta X\sim {{{\mathcal {O}}}}(0.1)\) is unlikely to be observed at least in the near future. However, in an (optimistic!) scenario of fine tuned large NP couplings, \(\varepsilon _{VR} \sim \varDelta X\sim \mathcal{O}(1)\), the nature of neutrinos can become clear if only the decay spectrum is at all measurable, even with small accuracy. If it happens, it also points out to the type of the neutrino NP interactions, as only the right-handed vector \(\varepsilon _{VR}\) couplings can simultaneously be large and lead to difference between the Dirac and Majorana case.

One should also note that Figs. 5 and 6 were plotted by taking individual non-vanishing NP neutrino couplings one at a time. In general, the real shape of spectrum can depend on several parameters and differ from those displayed in Figs. 5 and 6. In any specific NP model the effective shape is calculable and if experiment is able to measure it, fitting procedures can give the information on model parameters. However, as usual just discovery of discrepancy with the SM predictions does not immediately disclose what kind of NP interactions may be involved and further studies with combining more observables are necessary.

As the experimental accuracy of these s distributions gets improved, one could be optimistic that in future with enough of detected events, such a study might become feasible. We hope that the possibility of discovering NP as well as distinguishing the Dirac and Majorana nature of neutrinos based on their different NP behaviours, would lend support to any proposal for further experimental exploration of such rare decays as \(B \rightarrow K \, \nu \, {\overline{\nu }}\) and \(K \rightarrow \pi \, \nu \, {\overline{\nu }}\).

6 Conclusions

In this paper, we have provided a general framework for analysing the processes containing a neutrino antineutrino pair in the final state, considering neutrinos to be either Dirac or Majorana fermions. We emphasise the fact that the quantum mechanically identical nature of Majorana neutrino and antineutrino does not depend on the size of their mass. Thus, the different statistical properties of the Dirac and Majorana neutrinos can lead, at least in the New Physics models with non-standard neutrino couplings, to significantly different predictions for the decay rates and for the spectra of visible particles. Such differences are (at least in principle) not proportional to the tiny neutrino masses, unlike the huge suppression inherently present in the neutrino-mediated Lepton Number Violating processes.

We argue that within the SM due to the particular structure of the neutrino vector-like left-handed only couplings, such not-suppressed effects of different statistics unfortunately still disappear when integrating over the phase space of final state neutrinos (excluding eventually some special kinematic scenarios where the neutrino momenta can be reconstructed or inferred, see Ref. [37] for more details). However, we point out that, when New Physics effects are taken into consideration, the difference between Dirac and Majorana neutrino possibilities could survive even if one neglects neutrino mass dependent terms and integrates over the phase space of final state neutrinos.

We explicitly illustrate the general analysis by examples of neutral current processes with the 2-body or 3-body final states. Using our example processes \(Z \rightarrow \nu \, {\overline{\nu }}\) and \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu \, {\overline{\nu }}\) decays, we show how the experimental bounds on the NP couplings differ between Dirac neutrino and Majorana neutrino possibilities. We also point out that if the NP terms are substantial comparing to SM couplings and the differential decay rate or equivalently the missing mass-square distribution for \({\mathcal {P}}_i \rightarrow {\mathcal {P}}_f \, \nu \, {\overline{\nu }}\) decay can be measured accurately enough, the shape of the spectrum can directly help discern whether neutrinos are Dirac or Majorana fermions, at least in NP models were there are no new light exotic states which could also escape detection.