Spacetime evolution of lepton number densities and wave packet-like effects for neutrino flavor and chiral oscillations in quantum field theory

Apriadi Salim Adam [email protected] Research Center for Quantum Physics, National Research and Innovation Agency (BRIN), South Tangerang 15314, Indonesia Nicholas J. Benoit [email protected] Physics Program, Graduate School of Advanced Science and Engineering,
Hiroshima University, Higashi-Hiroshima 739-8526, Japan Yuta Kawamura [email protected] Yokkaichi city, Mie, Japan Yamato Matsuo [email protected] Akashi city, Hyogo, Japan Takuya Morozumi [email protected] Physics Program, Graduate School of Advanced Science and Engineering,
Hiroshima University, Higashi-Hiroshima 739-8526, Japan Core of Research for the Energetic Universe, Hiroshima University,
Higashi-Hiroshima 739-8526, Japan Yusuke Shimizu [email protected] Physics Program, Graduate School of Advanced Science and Engineering,
Hiroshima University, Higashi-Hiroshima 739-8526, Japan Core of Research for the Energetic Universe, Hiroshima University,
Higashi-Hiroshima 739-8526, Japan Naoya Toyota [email protected] Hiroshima city, Hiroshima, Japan

Abstract

We present a formulation of lepton family numbers, based on quantum field theory, for neutrino oscillation phenomenology that can be applied to nonrelativistic and relativistic energies for neutrinos. It is formulated for both types of neutrinos, Dirac and Majorana. The formulation is constructed as the time evolution of a lepton family number density operator. Then, the time evolution of the lepton family number density operator becomes dependent on the mass and new features appear. The expectation value of the density operator is evaluated for the initial state with a Gaussian distribution for the momentum amplitude. This enables us to study wave packet-like decoherence effects. We show in the nonrelativistic regime, the type of neutrino mass are distinguishable even under the presence of wave packet-like decoherence effects.

^†^†preprint: HUPD-2102

I Introduction

Neutrinos were originally formulated, within the weak interactions of the Standard Model (SM), to be massless fermions. However, hints of massive neutrinos appeared in the second half of the 20th century with the solar neutrino problem[1]. Eventually, updated solar models and experiments from the Kamiokande laboratory and the Sudbury Neutrino Observatory (SNO) proved the disappearance is caused by flavor oscillations[2, 3]. For neutrinos the existence of flavor oscillations mean there is a mixing between mass and flavor eigenstates. The mixing between mass and flavor is governed by the unitary PMNS matrix, which has six (four) free parameters in the $3\times 3$ instance[4, 5]. Thus, the mixing caused by flavor oscillations is direct evidence that neutrinos are massive and, consequently, require physics beyond the SM.

In this decade, the six parameters related to the PMNS matrix; the oscillation angles $\theta_{12},\theta_{23},\text{and }\theta_{13}$ , the mass squared differences $\Delta m_{21}^{2}\text{ and }\Delta m^{2}_{31}$ , and the Dirac CP violating phase $\delta_{cp}$ , are expected to be measured within a few percentages by neutrino oscillation experiments[6, 7, 8, 9, 10, 11]. In addition, the absolute mass of the lightest neutrino, $m_{1}$ or $m_{3}$ , is expected to be directly limited down to the sub-eV range by the Karlsruhe Tritium Neutrino Experiment (KATRIN)[12]. However, even with these precision measurements, questions remain in neutrino phenomenology. Specifically, the question of what type of mass neutrinos possess, Dirac or Majorana remains. This question could be answered by future neutrinoless double- $\beta$ $(0\nu\beta\beta)$ decay experiments[13, 14, 15, 16, 17]. Arguably, additional approaches to answer those questions would be ideal.

For this work, we present a unified description of neutrino phenomenology that leads to an additional approach for investigations of the previously mentioned questions. We define unified to mean a description that includes neutrino flavor oscillations and, particle-anti-particle or chiral oscillations. In addition, our description is different from the usual neutrino oscillation theory that assumes neutrinos are relativistic. The relativistic assumption is often taken, because cosmological limits on the neutrinos absolute masses place them in the sub-eV range[18]; whereas, experiments are performed on neutrinos with energies in the keV, MeV, GeV, and recently the PeV ranges[19]. Nonrelativistic neutrinos are predicted to exist in nature as the cosmic neutrino background (C $\nu$ B). Experiments to detect the C $\nu$ B are under preparation (see for example, [20]), and we plan future studies of the C $\nu$ B as an application of our description.

For studies of the C $\nu$ B we need to build a framework that uses a momentum distribution [21]. Since the state of the C $\nu$ B can be described by a mixed state at their decoupling time, their momentum distribution should be given by a density matrix. For oscillation experiments, if the neutrino is produced with a localized space, it should be described by a pure state with a momentum distribution. Descriptions of that type have been done using wave packet formulations [22, 23, 24, 25, 26]. To incorporate a momentum distribution in our framework, we consider densities of lepton family numbers that are localized in space as pure states. We have previously discussed details of the formulation for the Majorana mass case using plane waves [27, 28, 29].

In future work we can extend this study to the mixed state required by C $\nu$ B, and further consider the time evolution of lepton numbers under the expansion of the universe. Then, we may be able to clarify when the C $\nu$ B transits from the relativistic regime to the nonrelativistic regime. One may also find how the coherence for a given momentum distribution will continue or will be lost as the C $\nu$ B is redshifted due to the expansion. Therefore, our present work can be useful to clarify those natures of the C $\nu$ B.

II Brief outline of related studies

Consider the decay of a positively charged pion at rest to an anti-muon and neutrino,

\pi^{+}\rightarrow\mu^{+}+\nu_{\alpha}.

(1)

The Standard Model defines the neutrino flavor, $\alpha$ , to be a muon-neutrino $\nu_{\mu}$ . This definition is from connecting the neutrino and charged lepton with the electroweak doublets, i.e., $(\nu_{e},e)$ , $(\nu_{\mu},\mu)$ , and $(\nu_{\tau},\tau)$ . Then, theory introduced lepton family numbers as a way to express those connections; $L_{e}$ , $L_{\mu}$ , and $L_{\tau}$ [30, 31]. A consequence is the Standard Model conserves lepton family numbers in all processes.

Yet, neutrino masses imply lepton family numbers are not necessarily conserved [32]. This is evident in neutrino oscillation experiments. Consider the example of Eq.(2), which uses the produced anti-muon to identify a muon neutrino at the production point. Then after the neutrino propagates over a macroscopic distance, either a muon neutrino or electron neutrino is identified at a detector.

\pi^{+}\rightarrow\mu^{+}+\nu_{\mu}\xrightarrow{\text{propagation}}\begin{cases}\nu_{\mu}+X\rightarrow\mu^{-}+Y&\text{disappearance,}\\ \nu_{e}+X\rightarrow e^{-}+Y&\text{appearance.}\end{cases}

(2)

If we compare the family of the charged leptons from production to the appearance detection, a change of $L_{\mu}$ to $L_{e}$ occurs. Models explain this change by treating the interacting and propagating neutrinos differently.

In the quantum mechanics model, the interacting neutrinos are from a flavor basis and the propagating neutrinos are from a mass basis [33, 34]. Those bases are related to each other by a unitary transformation with a unitary matrix called the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, $U^{\ast}_{\alpha i}$ . Then, a neutrino flavor state $|\nu_{\alpha}\rangle$ is written as a coherent superposition of massive states weighted by the PMNS matrix,

|\nu_{\alpha}\rangle=\sum^{3}_{i=1}U^{\ast}_{\alpha i}|\nu_{i}\rangle.

(3)

Each forms an orthogonal basis $\langle\nu_{\alpha}|\nu_{\beta}\rangle=\delta_{\alpha\beta}$ and $\langle\nu_{i}|\nu_{j}\rangle=\delta_{i,j}$ . Sometimes these flavor states are called Pontecorvo states in the literature. However, it was shown that the states of Eq.(3) can not be applied to quantum field theory [35, 36, 37]. Consequently, effort has been made to understand the nature of flavor in quantum field theory. We introduce some quantum field theory models that are relevant to our work.

A common quantum field theory model treats the produced or detected neutrino flavor state is different from the Pontecorvo states in Eq.(3). The difference appears in the amplitudes of the superposition to each mass eigenstate. They are not given by the PMNS matrix elements alone. They depend on the additional factor coming from either the production or detection amplitude of each mass eigenstate. Then, that model derives the neutrino flavor states based on the production or detection processes [38, 39]. As a result, the neutrino flavor states are process dependent and modify the amplitudes of neutrino flavor oscillation probabilities. In the appearance experiment of Eq.(2), it also results in the detected neutrino flavor state and the produced neutrino flavor state being non-orthogonal to each other even before oscillations start. Nevertheless, if the neutrinos are treated as ultra-relativistic particles, the modifications are lost.

In a different model, wave packets for the external particles are introduced, and the neutrinos appear as intermediate states between production and detection [40, 41, 42, 43, 44]. For this model, no neutrino flavor states are considered or calculated. Only the source current of the neutrino and the detector current appear as physical particles. If the method is applied to Majorana neutrinos, two processes are required. The first process is to detect the chirality and lepton number conserving propagation. The second process is for chirality and lepton number violating propagation of the intermediate neutrinos. In our framework, both of the chirality conserving and violating effects are included in the lepton family numbers.

Lastly, another model considers the neutrino flavor states to be from a physical Fock space. The resulting neutrino flavor states are related to massive neutrino states by a Bogoliubov transformation, which leads to inequivalent vacua [45, 46, 47, 48, 49]. This model introduces non-zero Dirac mass to each flavor neutrino field [50, 51, 52]. A flavor charge is defined through vector current of the flavor neutrino. In our framework, the flavor charge for active neutrinos is defined based on the weak eigenstates that form SU(2) doublets with the mass eigenstates of the charged leptons. Therefore, the charge is defined through the left-handed current of the weak eigenstates. We call this charge lepton family number. Then in the formula for the lepton family number, only the neutrino masses for the mass eigenstates appear and additional mass parameters for the flavor neutrinos are not introduced.

III Lepton Number Density

We take an approach to define lepton family numbers for neutrinos based on the isospin doublet used for charged current weak interactions. Meaning, in the diagonal basis for charged lepton mass, the lepton family number is assigned to each lepton doublet $\Psi_{L\alpha}=\begin{pmatrix}\nu_{L\alpha}&e_{L\alpha}\end{pmatrix}^{T}$ as well as the right-handed charged lepton $e_{R\alpha}$ . Definitely, under a $U_{\alpha}$ (1) transformation of the lepton family number $L_{\alpha}$ , they transform as,

\displaystyle\Psi^{\prime}_{L\alpha}=\begin{pmatrix}\nu^{\prime}_{L\alpha}\\ e^{\prime}_{L\alpha}\end{pmatrix}=e^{i\theta_{\alpha}}\Psi_{L\alpha}=e^{i\theta_{\alpha}}\begin{pmatrix}\nu_{L\alpha}\\ e_{L\alpha}\end{pmatrix},

\displaystyle e^{\prime}_{R\alpha}=e^{i\theta\alpha}e_{R\alpha}.

(4)

Where the subscript $L$ and $R$ denote the left-handed projection operator $P_{L}=\frac{1-\gamma^{5}}{2}$ and $P_{R}=\frac{1+\gamma^{5}}{2}$ , respectively. The associated Noether current is defined by

l^{Noether}_{\mu\alpha}=:\overline{e_{\alpha}}\gamma_{\mu}e_{\alpha}:+:\overline{\nu_{\alpha}}\gamma_{\mu}P_{L}\nu_{\alpha}:,

(5)

where the colon : implies normal ordering. Note, the lepton family current consists of a vector current of the charged lepton and a left-handed current of the neutrino. We describe the effects of mass of neutrinos on lepton number in this particular weak basis where the charged lepton mass matrix is real and diagonal. Through the detection of a charged lepton flavor, the state of the neutrino is projected to the flavor associated to the charged lepton. In the process, the lepton family number is assumed to be conserved. For the example, $\pi^{+}\rightarrow\mu^{+}+\nu$ , the produced neutrino is assumed to have a muon number $+1$ . The projected state has a left-handed chirality with a definite lepton family number. We will see that the neutrino field with definite chirality can be expanded by the plane wave solution and spinors with helicity $-\frac{1}{2}$ . The effect of the non-zero mass of neutrino is encoded in non-trivial time dependence of creation and annihilation operators which we will evaluate below.

III.1 Majorana neutrino formulation

In the flavor basis where the charged lepton mass matrix is real and diagonal, the free part of the Lagrangian for Majorana neutrinos is,

\mathcal{L}^{M}=\overline{\nu_{L\alpha}}i\gamma^{\mu}\partial_{\mu}\nu_{L\alpha}-\frac{1}{2}\left(\overline{(\nu_{L\alpha})^{C}}m_{\alpha\beta}\nu_{L\beta}+\overline{\nu_{L\alpha}}(m^{\ast}_{\alpha\beta})(\nu_{L\beta})^{C}\right).

(6)

The equation of motion for the Majorana neutrino in the flavor basis is written as;

\begin{split}i\not{\partial}\nu_{L\alpha}&=m^{\ast}_{\alpha\beta}(\nu_{L\beta})^{C},\\ i\not{\partial}(\nu_{L\alpha})^{C}&=m_{\alpha\beta}(\nu_{L\beta}),\end{split}

(7)

where $m_{\alpha\beta}$ is a Majorana mass matrix. The Greek subscripts represent the flavors $e$ , $\mu$ , and $\tau$ . In the flavor basis, the Majorana mass matrix is complex and symmetric. In addition, we have used the notation $(\nu_{L\alpha})^{C}$ for charge conjugation. The solution of Eq.(7) can be written with the form of the expansion with the massless spinors plus a zero-mode contribution. The details on how we quantize this situation is written in appendix A. The remainder of the main text will focus on the massless spinor, i.e., the non-zero-mode contributions,

\nu_{L\alpha}(t,\mathbf{x})=\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2|\mathbf{p}|}\left(a_{\alpha}(\mathbf{p},t)u_{L}(\mathbf{p})e^{i\mathbf{p}\cdot\mathbf{x}}+b^{\dagger}_{\alpha}(\mathbf{p},t)v_{L}(\mathbf{p})e^{-i\mathbf{p}\cdot\mathbf{x}}\right),

(8)

where $\int^{\prime}$ implies the exclusion of zero momentum mode from the integration. The massless spinors $u_{L}(\mathbf{p})$ and $v_{L}(\mathbf{p})$ are given by [27],

	$\displaystyle u_{L}(\mathbf{p})=-v_{L}(\mathbf{p})=\sqrt{\lvert 2\mathbf{p}\rvert}\begin{pmatrix}0\\ \phi_{-}(\mathbf{n})\end{pmatrix},$		(9)
	$\displaystyle\mathbf{n}\cdot\bm{\sigma}\phi_{\pm}(\mathbf{n})=\pm\phi_{\pm}(\mathbf{n}),\quad\mathbf{n}=\frac{\mathbf{p}}{\lvert\mathbf{p}\rvert},$		(10)

where the two component spinors $\phi_{\pm}(\pm\mathbf{n})$ are written by the polar angle $\theta$ and the azimuthal angle $\phi$ specifying the direction of the momentum $\mathbf{p}$ ,

	$\displaystyle\phi_{+}(\mathbf{n})$	$\displaystyle=\begin{pmatrix}e^{-i\frac{\phi}{2}}\cos\frac{\theta}{2}\\ e^{i\frac{\phi}{2}}\sin\frac{\theta}{2}\end{pmatrix},$	$\displaystyle\phi_{-}(\mathbf{n})$	$\displaystyle=\begin{pmatrix}-e^{-i\frac{\phi}{2}}\sin\frac{\theta}{2}\\ e^{i\frac{\phi}{2}}\cos\frac{\theta}{2}\end{pmatrix},$		(11)
	$\displaystyle\phi_{+}(-\mathbf{n})$	$\displaystyle=i\phi_{-}(\mathbf{n}),$	$\displaystyle\phi_{-}(-\mathbf{n})$	$\displaystyle=i\phi_{+}(\mathbf{n}).$		(12)

We then formulate a Heisenberg operator for $L^{M}_{\alpha}$ where $\alpha=(e,\mu,\tau)$ , from the lepton number density $l^{M}_{\alpha}(t,\mathbf{x})$ by integrating over space,

	$\displaystyle L^{M}_{\alpha}(t)=\int d^{3}x\,l^{M}_{\alpha}(t,\mathbf{x})=\int^{\prime}\frac{d^{3}k}{(2\pi)^{3}2\lvert\mathbf{k}\rvert}\left[a^{\dagger}_{\alpha}(\mathbf{k},t)a_{\alpha}(\mathbf{k},t)-b^{\dagger}_{\alpha}(\mathbf{k},t)b_{\alpha}(\mathbf{k},t)\right],$		(13)
	$\displaystyle l^{M}_{\alpha}(t,\mathbf{x})=\,:\overline{\nu_{L\alpha}}(t,\mathbf{x})\gamma^{0}\nu_{L\alpha}(t,\mathbf{x}):\,,$		(14)

where the colon : denotes normal ordering. In the absence of mass matrix, the time evolution of creation and annihilation operators are

\displaystyle a_{\alpha}(\mathbf{p},t)=a_{\alpha}(\mathbf{p},t_{0})e^{-i|\mathbf{p}|(t-t_{0})},

\displaystyle b^{\dagger}_{\alpha}(\mathbf{p},t)=b^{\dagger}_{\alpha}(\mathbf{p},t_{0})e^{i|\mathbf{p}|(t-t_{0})}.

(15)

In this case, all three lepton family number operators are separately conserved,

L^{M}_{\alpha}(t)=L^{M}_{\alpha}(t_{0})

(16)

However as can be seen from Eq.(7), their time dependence is governed by non-diagonal mass matrix. To obtain the time dependence of the lepton number operator in Eq.(13), we need to know the time dependence of the flavor operators $a_{\alpha}(\mathbf{k},t)$ and $b_{\alpha}(\mathbf{k},t)$ . In order to evolve the operators to time $t$ from the initial time $t_{0}$ , one needs to express them at $t_{0}$ as a superposition of operators with definite masses. This is because the operators with definite masses $m_{i}$ evolve like $e^{\pm iE_{i}t}$ , with $E_{i}=\sqrt{|\mathbf{p}|^{2}+m_{i}^{2}}$ . Then the time dependence of the flavor operator can be derived.

We begin by diagonalizing the mass matrix in Eq.(7) with a unitary matrix $V$ as,

	$\displaystyle m_{i}\delta_{ij}=\left(V^{T}\right)_{i\alpha}m_{\alpha\beta}V_{\beta j},$		(17)
	$\displaystyle\nu_{L\alpha}=V_{\alpha i}\nu_{Li},$		(18)

where $\nu_{Li}$ denotes an operator of the mass basis. $\nu_{Li}$ and its charge conjugation form the operator for the Majorana field as,

\begin{split}\psi_{Mi}({\bf x},t)&=\nu_{Li}({\bf x},t)+(\nu_{Li}({\bf x},t))^{c}\\ &=V_{\alpha i}^{\ast}\nu_{\alpha L}({\bf x},t)+V_{\alpha i}(\nu_{\alpha L}({\bf x},t))^{c}.\end{split}

(19)

Since Eq.(19) satisfies the Majorana condition, it can be expanded as

\psi_{Mi}({\bf x},t)=\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2E_{i}(\mathbf{p})}\sum_{\lambda=\pm}\Bigl{(}a_{Mi}(\mathbf{p},\lambda)u_{i}(\mathbf{p},\lambda)e^{-i(E_{i}t-\mathbf{p}\cdot\mathbf{x})}+a^{\dagger}_{Mi}(\mathbf{p},\lambda)v_{i}(\mathbf{p},\lambda)e^{i(E_{i}t-\mathbf{p}\cdot\mathbf{x})}\Bigr{)}.

(20)

where $E_{i}=\sqrt{|\mathbf{p}|^{2}+m_{i}^{2}}$ denotes the energy of the mass eigenstate and $u_{i}(\mathbf{p},\lambda)$ and $v_{i}(\mathbf{p},\lambda)$ denote the massive Dirac spinors with the definite helicity $\lambda$ and mass $m_{i}$ . For the details of the definitions, see appendix A in Ref.[27]. The operators $a_{Mi}(\mathbf{p},\lambda)$ of Eq.(20) are distinct from the operators with the massless spinors in Eq.(8). We will call $a_{Mi}(\mathbf{p},\lambda)$ the Majorana operators, and they obey the anti-commutation relation,

\{a_{Mi}(\mathbf{p},\lambda),a^{\dagger}_{Mj}(\mathbf{q},\lambda^{\prime})\}=2E_{i}(\mathbf{p})(2\pi)^{3}\delta^{(3)}(\mathbf{p}-\mathbf{q})\delta_{ij}\delta_{\lambda\lambda^{\prime}}\,,

(21)

with all others being zero.

Using Eq.(19) with the expansions Eq.(20) and Eq.(8), one can derive the relation between the operators $a_{\alpha}(\mathbf{p},t),b_{\alpha}(\mathbf{p},t)$ and the operators for massive fields $a_{Mi}(\mathbf{p},\lambda)$ for arbitrary time $t$ .

	$\displaystyle a_{\alpha}(\pm\mathbf{p},t)=V_{\alpha i}\frac{\sqrt{2\lvert\mathbf{p}\rvert(E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert)}}{2E_{i}(\mathbf{p})}\left(a_{Mi}(\pm\mathbf{p},-)e^{-iE_{i}(\mathbf{p})t}\pm\frac{im_{i}}{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}a^{\dagger}_{Mi}(\mp\mathbf{p},-)e^{iE_{i}(\mathbf{p})t}\right),$		(22)
	$\displaystyle b_{\alpha}(\pm\mathbf{p},t)=V_{\alpha i}\frac{\sqrt{2\lvert\mathbf{p}\rvert(E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert)}}{2E_{i}(\mathbf{p})}\left(a_{Mi}(\pm\mathbf{p},+)e^{-iE_{i}(\mathbf{p})t}\pm\frac{im_{i}}{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}a^{\dagger}_{Mi}(\mp\mathbf{p},+)e^{iE_{i}(\mathbf{p})t}\right),$		(23)

where $+\mathbf{p}$ denotes the momentum directed toward the positive y-axis hemisphere $(0\leq\phi<\pi)$ and $-\mathbf{p}$ denotes the momentum directed toward the negative y-axis hemisphere $(\pi\leq\phi<2\pi)$ . Notice, the relations are a non-trivial mixing of the Majorana annihilation and creation operators forming $a_{\alpha}(\pm\mathbf{p},t)$ and $b_{\alpha}(\pm\mathbf{p},t)$ . Furthermore, we can identify the relations as a Bogoliubov transformation that occurs between the operators [53]. The operators $a_{\alpha}(\pm\mathbf{p},t)$ and $b_{\alpha}(\pm\mathbf{p},t)$ , in Eq.(22) and Eq.(23), satisfy the anti-commutation relations,

\{a_{\alpha}(\pm\mathbf{p},t),a^{\dagger}_{\beta}(\pm\mathbf{q},t)\}=\{b_{\alpha}(\pm\mathbf{p},t),b^{\dagger}_{\beta}(\pm\mathbf{q},t)\}=2\lvert\mathbf{p}\rvert(2\pi)^{3}\delta^{(3)}(\mathbf{p}-\mathbf{q})\delta_{\alpha\beta}\,,

(24)

with all other relations being zero. The inverse relations at $t=t_{0}$ are also obtained as,

	$\displaystyle a_{Mi}(\pm\mathbf{p},-)=e^{iE_{i}(\mathbf{p})t_{0}}\sqrt{\frac{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}{2\lvert\mathbf{p}\rvert}}\left(V^{\ast}_{\beta i}a_{\beta}(\pm\mathbf{p},t_{0})\mp V_{\beta i}\frac{im_{i}}{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}a^{\dagger}_{\beta}(\mp\mathbf{p},t_{0})\right),$		(25)
	$\displaystyle a^{\dagger}_{Mi}(\pm\mathbf{p},-)=e^{-iE_{i}(\mathbf{p})t_{0}}\sqrt{\frac{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}{2\lvert\mathbf{p}\rvert}}\left(V_{\beta i}a^{\dagger}_{\beta}(\pm\mathbf{p},t_{0})\pm V^{\ast}_{\beta i}\frac{im_{i}}{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}a_{\beta}(\mp\mathbf{p},t_{0})\right),$		(26)
	$\displaystyle a_{Mi}(\pm\mathbf{p},+)=e^{iE_{i}(\mathbf{p})t_{0}}\sqrt{\frac{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}{2\lvert\mathbf{p}\rvert}}\left(V_{\beta i}{b}_{\beta}(\pm\mathbf{p},t_{0})\mp V^{\ast}_{\beta i}\frac{im_{i}}{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}b^{\dagger}_{\beta}(\mp\mathbf{p},t_{0})\right),$		(27)
	$\displaystyle a^{\dagger}_{Mi}(\pm\mathbf{p},+)=e^{-iE_{i}(\mathbf{p})t_{0}}\sqrt{\frac{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}{2\lvert\mathbf{p}\rvert}}\left(V^{\ast}_{\beta i}b^{\dagger}_{\beta}(\pm\mathbf{p},t_{0})\pm V_{\beta i}\frac{im_{i}}{E_{i}(\mathbf{p})+\lvert\mathbf{p}\rvert}b_{\beta}(\mp\mathbf{p},t_{0})\right).$		(28)

By substituting the inverse relations Eqs.(25-28) into Eq.(22) and Eq.(23), we obtain relations between the arbitrary time operators $(a_{\alpha}(\pm\mathbf{p},t),b_{\alpha}(\pm\mathbf{p},t))$ and the initial time operators $(a_{\alpha}(\pm\mathbf{p},t_{0}),b_{\alpha}(\pm\mathbf{p},t_{0}))$ .

	$\displaystyle\begin{split}a_{\alpha}(\pm\mathbf{p},t)=\sum_{\beta=e}^{\tau}\sum_{j}\biggl{(}V_{\alpha j}V^{\ast}_{\beta j}&\left[\cos[{E_{j}(\mathbf{p})(t-t_{0})}]-i\frac{\lvert\mathbf{p}\rvert}{E_{j}(\mathbf{p})}\sin[{E_{j}(\mathbf{p})(t-t_{0})}]\right]a_{\beta}(\pm\mathbf{p},t_{0})\biggr{.}\\ &\biggl{.}\mp V_{\alpha j}V_{\beta j}\frac{m_{j}}{E_{j}(\mathbf{p})}\sin[{E_{j}(\mathbf{p})(t-t_{0})}]a^{\dagger}_{\beta}(\mp\mathbf{p},t_{0})\biggr{)},\end{split}$		(29)
	$\displaystyle\begin{split}b_{\alpha}(\pm\mathbf{p},t)=\sum_{\gamma=e}^{\tau}\sum_{j}\biggl{(}V^{\ast}_{\alpha j}V_{\gamma j}&\left[\cos[{E_{j}(\mathbf{p})(t-t_{0})}]-i\frac{\lvert\mathbf{p}\rvert}{E_{j}(\mathbf{p})}\sin[{E_{j}(\mathbf{p})(t-t_{0})}]\right]b_{\gamma}(\pm\mathbf{p},t_{0})\biggr{.}\\ &\biggl{.}\mp V^{\ast}_{\alpha j}V^{\ast}_{\gamma j}\frac{m_{j}}{E_{j}(\mathbf{p})}\sin[{E_{j}(\mathbf{p})(t-t_{0})}]b^{\dagger}_{\gamma}(\mp\mathbf{p},t_{0})\biggr{)}.\end{split}$		(30)

The equations (29) and (30) tell us the relations between the operators $(a_{\alpha}(\pm\mathbf{p},t),b_{\alpha}(\pm\mathbf{p},t))$ and $(a_{\alpha}(\pm\mathbf{p},t_{0}),b_{\alpha}(\pm\mathbf{p},t_{0}))$ depend on time through the difference $T\equiv t-t_{0}$ .

We rewrite the Majorana density operator of Eq.(14) using the expansion with massless spinors in Eq.(8);

l^{M}_{\alpha}(t,\mathbf{x})=\int^{\prime}\frac{d^{3}\mathbf{k}}{(2\pi)^{3}2|\mathbf{k}|}\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2|\mathbf{p}|}\\ \times\left[a^{\dagger}_{\alpha}(\mathbf{k},t)a_{\alpha}(\mathbf{p},t)\overline{u_{L}}(\mathbf{k})\gamma^{0}u_{L}(\mathbf{p})e^{-i(\mathbf{k}-\mathbf{p})\cdot\mathbf{x}}+b_{\alpha}(\mathbf{k},t)a_{\alpha}(\mathbf{p},t)\overline{v_{L}}(\mathbf{k})\gamma^{0}u_{L}(\mathbf{p})e^{i(\mathbf{k}+\mathbf{p})\cdot\mathbf{x}}\right.\\ +\left.a^{\dagger}_{\alpha}(\mathbf{k},t)b^{\dagger}_{\alpha}(\mathbf{p},t)\overline{u_{L}}(\mathbf{k})\gamma^{0}v_{L}(\mathbf{p})e^{-i(\mathbf{k}+\mathbf{p})\cdot\mathbf{x}}-b^{\dagger}_{\alpha}(\mathbf{p},t)b_{\alpha}(\mathbf{k},t)\overline{v_{L}}(\mathbf{k})\gamma^{0}v_{L}(\mathbf{p})e^{i(\mathbf{k}-\mathbf{p})\cdot\mathbf{x}}\right],

(31)

where the integration region $\int^{\prime}$ does not include the zero momentum mode. Next, we can write the time evolution the Majorana density, solely in terms of the operators $a_{\alpha}(\mathbf{p},t_{0})$ and $b_{\alpha}(\mathbf{p},t_{0})$ by substituting the relations of Eq.(29) and Eq.(30). The result for the time evolution of Eq.(13) was first presented in Ref.[27].

To study the spacetime evolution, we take the expectation value of the Majorana density operator $\langle\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}|l^{M}_{\alpha}(t,\mathbf{x})|\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ . We assume the initial momentum state of $\lvert\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ is sharply peaked at zero in the first and third components, but has a Gaussian distribution in the second, i.e., $\mathbf{q}=(0,q,0)$ . This leads to a shape similar to a 1-D Gaussian wave packet,

\lvert\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle=\frac{1}{\sqrt{\sigma_{q}}(2\pi)^{3/4}}\int^{\prime}\frac{dq}{\sqrt{A}\sqrt{2|q|}}e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}}a_{\sigma}^{\dagger}(\mathbf{q},t_{0})\lvert 0(t_{0})\rangle;,

(32)

where $A=(2\pi)^{2}\delta^{2}(0)$ denotes the area for the two-dimensional space $(x_{1},x_{3})$ perpendicular to the direction of the momentum. This is a normalization factor that leads to $\langle\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rvert L_{\sigma}(t=t_{0})\lvert\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle=1$ . The vacuum $\lvert 0(t_{0})\rangle$ satisfies,

a(\mathbf{p},t_{0})\lvert 0(t_{0})\rangle=b(\mathbf{p},t_{0})\lvert 0(t_{0})\rangle=0,

(33)

for all non-zero $\mathbf{p}$ . In Eq.(32), $\sigma_{q}$ is the width of the Gaussian distribution in the second component of the momentum. The second component of the mean momentum $q^{0}$ is positive for the Gaussian distribution. Sandwiching the Majorana density operator of Eq.(31) with Eq.(32) and taking the integration over $\mathbf{k}$ and $\mathbf{p}$ results in,

\langle\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}|l^{M}_{\alpha}(t,\mathbf{x})|\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle=\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}\frac{dq^{\prime}dq}{A}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)\mathbf{e}_{2}\cdot\mathbf{x}}\\ \times\left[\sum_{i,j}V^{\ast}_{\alpha i}V_{\sigma i}V_{\alpha j}V^{\ast}_{\sigma j}\left(\cos{E_{i}(q^{\prime})T}+i\frac{|q^{\prime}|}{E_{i}(q^{\prime})}\sin{E_{i}(q^{\prime})T}\right)\left(\cos{E_{j}(q)T}-i\frac{|q|}{E_{j}(q)}\sin{E_{j}(q)T}\right)\right.\\ \left.-\sum_{i,j}V^{\ast}_{\alpha i}V^{\ast}_{\sigma i}V_{\alpha j}V_{\sigma j}\frac{m_{j}}{E_{j}(q^{\prime})}\sin{E_{j}(q^{\prime})T}\frac{m_{i}}{E_{i}(q)}\sin{E_{i}(q)T}\right];

(34)

where $T=t-t_{0}$ is the time difference and $\mathbf{e}^{T}_{2}=(0,1,0)$ is the unit vector. The integration over $x_{1}$ and $x_{3}$ , results in a linear density, which we denote as

\lambda^{M}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})\equiv\iint dx_{1}dx_{3}\langle\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}|l^{M}_{\alpha}(t,\mathbf{x})|\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle.

(35)

The result of the integration is,

\lambda^{M}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})=\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}{dq^{\prime}dq}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)x_{2}}\\ \times\left[\sum_{i,j}V^{\ast}_{\alpha i}V_{\sigma i}V_{\alpha j}V^{\ast}_{\sigma j}\left(\cos{E_{i}(q^{\prime})T}+i\frac{|q^{\prime}|}{E_{i}(q^{\prime})}\sin{E_{i}(q^{\prime})T}\right)\left(\cos{E_{j}(q)T}-i\frac{|q|}{E_{j}(q)}\sin{E_{j}(q)T}\right)\right.\\ \left.-\sum_{i,j}V^{\ast}_{\alpha i}V^{\ast}_{\sigma i}V_{\alpha j}V_{\sigma j}\frac{m_{j}}{E_{j}(q^{\prime})}\sin{E_{j}(q^{\prime})T}\frac{m_{i}}{E_{i}(q)}\sin{E_{i}(q)T}\right];

(36)

where $x_{2}=\mathbf{e}_{2}\cdot\mathbf{x}$ . We note that a factor of $1/A$ is absent in the linear density because $\iint dx_{1}dx_{3}(1/A)=1$ .

To perform the integration over $q^{\prime}$ and $q$ in Eq.(36) we must assume two properties about the Gaussian distributions. First, the distributions are sharply peaked around the mean momentum value $q^{0}$ , i.e., $\sigma_{q}\ll q^{0}$ . Second, the width (variance) of the distribution $\sigma_{q}$ does not change in spacetime. Those two assumptions allow us to approximate the $q^{\prime}$ and $q$ integration as Gaussian; because,

E_{i,j}(q^{(\prime)})\simeq E_{i,j}(q^{0})+\frac{q^{0}}{E_{i,j}(q^{0})}(q^{(\prime)}-q^{0}).

(37)

The Gaussian integration of Eq.(36) leads to,

\lambda^{M}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})\simeq\frac{\sigma_{q}}{(2\pi)^{1/2}}\sum_{i,j}V^{\ast}_{\alpha i}V_{\sigma i}V_{\alpha j}V^{\ast}_{\sigma j}\\ \times\frac{1}{2}\left[(v_{i0}+v_{j0}+1+v_{i0}v_{j0})e^{i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}-v_{i0}T)^{2}+(x_{2}-v_{j0}T)^{2}]}\right.\\ \quad\quad-(v_{i0}+v_{j0}-1-v_{i0}v_{j0})e^{-i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}+v_{i0}T)^{2}+(x_{2}+v_{j0}T)^{2}]}\\ \quad\quad+(v_{i0}-v_{j0}+1-v_{i0}v_{j0})e^{i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}-v_{i0}T)^{2}+(x_{2}+v_{j0}T)^{2}]}\\ \left.\quad\quad-(v_{i0}-v_{j0}-1+v_{i0}v_{j0})e^{-i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}+v_{i0}T)^{2}+(x_{2}-v_{j0}T)^{2}]}\right]\\ -\frac{\sigma_{q}}{(2\pi)^{1/2}}\sum_{i,j}V^{\ast}_{\alpha i}V^{\ast}_{\sigma i}V_{\alpha j}V_{\sigma j}\sqrt{1-v^{2}_{i0}}\sqrt{1-v^{2}_{j0}}\\ \times\frac{1}{2}\left[e^{i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}-v_{i0}T)^{2}+(x_{2}-v_{j0}T)^{2}]}\right.\\ \quad+e^{-i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}+v_{i0}T)^{2}+(x_{2}+v_{j0}T)^{2}]}\\ \quad-e^{i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}-v_{i0}T)^{2}+(x_{2}+v_{j0}T)^{2}]}\\ \left.-e^{-i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}+v_{i0}T)^{2}+(x_{2}-v_{j0}T)^{2}]}\right],

(38)

where $v_{i,j0}=q^{0}/E_{i,j}(q^{0})$ are the group velocities of the distributions. The terms with the PMNS matrix combination $V^{\ast}_{\alpha i}V^{\ast}_{\sigma i}V_{\alpha j}V_{\sigma j}$ are dependent on the Majorana phases and are suppressed by the small masses of neutrinos $\sqrt{1-v^{2}_{i,j0}}=m_{i,j}/E_{i,j}(q^{0})$ . Lastly, to compare with our previous work’s result that considered only the time evolution of plane waves, we integrate the linear density of Eq.(38) over all space,

\int\lambda^{M}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})dx_{2}\simeq\frac{1}{4}\sum_{i,j}V^{\ast}_{\alpha i}V_{\sigma i}V_{\alpha j}V^{\ast}_{\sigma j}\\ \times\left[(v_{i0}+v_{j0}+1+v_{i0}v_{j0})e^{i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}\frac{(v_{i0}+v_{j0})^{2}T^{2}}{2}}\right.\\ \quad\quad-(v_{i0}+v_{j0}-1-v_{i0}v_{j0})e^{-i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}\frac{(v_{i0}+v_{j0})^{2}T^{2}}{2}}\\ \quad\quad+(v_{i0}-v_{j0}+1-v_{i0}v_{j0})e^{i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}\frac{(v_{i0}-v_{j0})^{2}T^{2}}{2}}\\ \left.\quad\quad-(v_{i0}-v_{j0}-1+v_{i0}v_{j0})e^{-i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}\frac{(v_{i0}-v_{j0})^{2}T^{2}}{2}}\right]\\ -\frac{1}{4}\sum_{i,j}V^{\ast}_{\alpha i}V^{\ast}_{\sigma i}V_{\alpha j}V_{\sigma j}\sqrt{1-v^{2}_{i0}}\sqrt{1-v^{2}_{j0}}\\ \times\left[(e^{i(E_{i}(q^{0})-E_{j}(q^{0}))T}+e^{-i(E_{i}(q^{0})-E_{j}(q^{0}))T})e^{-\sigma^{2}_{q}\frac{(v_{i0}+v_{j0})^{2}T^{2}}{2}}\right.\\ \left.-(e^{i(E_{i}(q^{0})+E_{j}(q^{0}))T}+e^{-i(E_{i}(q^{0})+E_{j}(q^{0}))T})e^{-\sigma^{2}_{q}\frac{(v_{i0}-v_{j0})^{2}T^{2}}{2}}\right].

(39)

The result of our previous work [27] is recovered in the plane wave limit, for which we take the width $\sigma_{q}$ of the density in momentum space to zero.

III.2 Dirac neutrino formulation

For a Dirac neutrino, we start from the following Lagrangian,

\mathcal{L}^{D}=\overline{\nu_{L\alpha}}i\gamma^{\mu}\partial_{\mu}\nu_{L\alpha}+\overline{\nu_{R\alpha}}i\gamma^{\mu}\partial_{\mu}\nu_{R\alpha}-\left(\overline{\nu_{R\alpha}}m_{\alpha\beta}\nu_{L\beta}+\overline{\nu_{L\alpha}}(m^{\dagger})_{\alpha\beta}\nu_{R\beta}\right).

(40)

The first and second terms are kinetic ones and third term is Dirac mass one. The equation of motion is derived as,

	$\displaystyle i\gamma^{\mu}\partial_{\mu}\nu_{L\alpha}=(m^{\dagger})_{\alpha\beta}\nu_{R\beta}$		(41)
	$\displaystyle i\gamma^{\mu}\partial_{\mu}\nu_{R\alpha}=m_{\alpha\beta}\nu_{L\beta}.$		(42)

We expand left and right chiral fields by massless spinors

	$\displaystyle\nu_{L\alpha}(t,\mathbf{x})=\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2\|\mathbf{p}\|}\left(u_{L}(\mathbf{p})a_{L\alpha}(\mathbf{p},t)e^{i\mathbf{p}\cdot\mathbf{x}}+v_{L}(\mathbf{p})b_{L\alpha}^{\dagger}(\mathbf{p},t)e^{-i\mathbf{p}\cdot\mathbf{x}}\right),$		(43)
	$\displaystyle\nu_{R\alpha}(t,\mathbf{x})=\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2\|\mathbf{p}\|}\left(u_{R}(\mathbf{p})a_{R\alpha}(\mathbf{p},t)e^{i\mathbf{p}\cdot\mathbf{x}}+v_{R}(\mathbf{p})b_{R\alpha}^{\dagger}(\mathbf{p},t)e^{-i\mathbf{p}\cdot\mathbf{x}}\right).$		(44)

We have introduced a creation and annihilation operators for the fields for each chirality. The massless spinors for left chirality field $\nu_{L\alpha}$ are the same as Eq.(9), while the following massless spinors are used for right chirality field,

u_{R}(\mathbf{p})=-v_{R}(\mathbf{p})=\sqrt{\lvert 2\mathbf{p}\rvert}\begin{pmatrix}\phi_{+}(\mathbf{n})\\ 0\end{pmatrix}.

(45)

In contrast to the Majorana case, the Dirac mass matrix is diagonalized by two mixing matrices,

	$\displaystyle\nu_{L\beta}=V_{\beta j}\nu_{Lj},$		(46)
	$\displaystyle\nu_{R\alpha}=U_{\alpha i}\nu_{Ri},$		(47)
	$\displaystyle(U^{\dagger})_{i\alpha}m_{\alpha\beta}V_{\beta j}=m_{i}\delta_{ij}.$		(48)

Then the equations of motion, Eq.(41) and Eq.(42), for the mass basis turn into the following form,

	$\displaystyle i\gamma^{\mu}\partial_{\mu}\nu_{Li}=m_{i}\nu_{Ri},$		(49)
	$\displaystyle i\gamma^{\mu}\partial_{\mu}\nu_{Ri}=m_{i}\nu_{Li}.$		(50)

The Dirac field with the definite mass $m_{i}$ is constructed as

	$\displaystyle\psi_{Di}({\mathbf{x}},t)$	$\displaystyle=\nu_{Li}({\bf x},t)+\nu_{Ri}({\bf x},t)$		(51)
		$\displaystyle=V_{\alpha i}^{\ast}\nu_{L\alpha}({\bf x},t)+U_{\alpha i}^{\ast}\nu_{R\alpha}({\bf x},t).$		(52)

We expand the Dirac field in Eq.(51) using massive spinors,

\psi_{Di}({\mathbf{x}},t)=\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2E_{i}(\mathbf{p})}\sum_{\lambda}\left(u_{i}(\mathbf{p},\lambda)a_{i}(\mathbf{p},\lambda)e^{-i(E_{i}t-\mathbf{p}\cdot\mathbf{x})}+v_{i}(\mathbf{p},\lambda)b_{i}^{\dagger}(\mathbf{p},\lambda)e^{i(E_{i}t-\mathbf{p}\cdot\mathbf{x})}\right),

(53)

where $\lambda$ is the helicity. The mass operators $a_{i}(\mathbf{p},\lambda)$ and $b_{i}(\mathbf{p},\lambda)$ obey the usual anti-commutation relations,

	$\displaystyle\{a_{i}(\mathbf{p},\lambda),a_{j}^{\dagger}(\mathbf{q},\lambda)\}$	$\displaystyle=(2\pi)^{3}2E_{i}(\mathbf{p})\delta_{ij}\delta^{(3)}(\mathbf{p}-\mathbf{q}),$		(54)
	$\displaystyle\{b_{i}(\mathbf{p},\lambda),b_{j}^{\dagger}(\mathbf{q},\lambda)\}$	$\displaystyle=(2\pi)^{3}2E_{i}(\mathbf{p})\delta_{ij}\delta^{(3)}(\mathbf{p}-\mathbf{q}).$		(55)

All other anti-commutation relations are zero. Using Eq.(43), Eq.(44), and Eqs.(51-53), the operators for the flavor basis can be related to the mass basis.

	$\displaystyle\frac{1}{\sqrt{2\|\mathbf{p}\|}}a_{L\alpha}(\pm\mathbf{p},t)$	$\displaystyle=\sum_{i}^{3}V_{\alpha i}\frac{\sqrt{E_{i}(\mathbf{p})+\|\mathbf{p}\|}}{2E_{i}(\mathbf{p})}\left(a_{i}(\pm\mathbf{p},-)e^{-iE_{i}(\mathbf{p})t}\pm i\frac{m_{i}}{E_{i}(\mathbf{p})+\|\mathbf{p}\|}b_{i}^{\dagger}(\mp\mathbf{p},-)e^{iE_{i}(\mathbf{p})t}\right),$		(56)
	$\displaystyle\frac{1}{\sqrt{2\|\mathbf{p}\|}}b_{L\alpha}^{\dagger}(\pm\mathbf{p},t)$	$\displaystyle=\sum_{i}^{3}V_{\alpha i}\frac{\sqrt{E_{i}(\mathbf{p})+\|\mathbf{p}\|}}{2E_{i}(\mathbf{p})}\left(b_{i}^{\dagger}(\pm\mathbf{p},+)e^{iE_{i}(\mathbf{p})t}\mp i\frac{m_{i}}{E_{i}(\mathbf{p})+\|\mathbf{p}\|}a_{i}(\mp\mathbf{p},+)e^{-iE_{i}(\mathbf{p})t}\right).$		(57)

Where the relations for the right-handed operators $a_{R\alpha},\,b_{R\alpha}$ are found by replacing the PMNS matrices $V\rightarrow U$ and flipping the operator helicity $a_{i}(\pm\mathbf{p},\pm)\rightarrow a_{i}(\pm\mathbf{p},\mp)$ , $b_{i}(\pm\mathbf{p},\pm)\rightarrow b_{i}(\pm\mathbf{p},\mp)$ ¹¹1The Dirac case operator relations first appeared in the thesis of [29].. Again, the operators obey the usual anti-commutation relations of,

	$\displaystyle\{a_{L\alpha}(\pm\mathbf{p},t),a_{L\beta}^{\dagger}(\pm\mathbf{q},t)\}$	$\displaystyle=2\|\mathbf{p}\|(2\pi)^{3}\delta^{(3)}(\mathbf{p}-\mathbf{q})\delta_{\alpha\beta},$		(58)
	$\displaystyle\{b_{L\alpha}(\pm\mathbf{p},t),b_{L\beta}^{\dagger}(\pm\mathbf{q},t)\}$	$\displaystyle=2\|\mathbf{p}\|(2\pi)^{3}\delta^{(3)}(\mathbf{p}-\mathbf{q})\delta_{\alpha\beta},$		(59)

with all others being zero. Similar anti-commutation relations hold for the right-handed operators.

The time evolution of Eq.(56) and Eq.(57) starting from operators defined at $t=t_{0}$ , can be calculated using a similar method to the Majorana case, Eq.(29) and Eq.(30). We first write the massive operators $a_{i}(\mathbf{p},\lambda)$ and $b^{\dagger}_{i}(\mathbf{p},\lambda)$ in terms of the operators $a_{L\alpha}(\mathbf{p},t_{0})$ , $b_{L\alpha}(\mathbf{p},t_{0})$ , $a_{R\alpha}(\mathbf{p},t_{0})$ , and $b_{R\alpha}(\mathbf{p},t_{0})$ using the inverted relations of Eq.(56), Eq.(57), and their right-handed counterparts. Then we substitute those inverted relations at $t=t_{0}$ into the massive operators of Eq.(56) and Eq.(57) to find the time evolution to be,

	$\displaystyle\begin{split}a_{L\alpha}(\pm\mathbf{p},t)=\sum_{i=1}^{3}\sum_{\beta=e}^{\tau}&\left[V_{\alpha i}V^{\ast}_{\beta i}\left(\cos E_{i}(\mathbf{p})T-i\frac{\|\mathbf{p}\|}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\right)a_{L\beta}(\pm\mathbf{p},t_{0})\right.\\ &\left.\quad\mp V_{\alpha i}U^{\ast}_{\beta i}\frac{m_{i}}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\,b^{\dagger}_{R\beta}(\mp\mathbf{p},t_{0})\right],\end{split}$		(60)
	$\displaystyle\begin{split}a^{\dagger}_{L\alpha}(\pm\mathbf{p},t)=\sum_{i=1}^{3}\sum_{\gamma=e}^{\tau}&\left[V^{\ast}_{\alpha i}V_{\gamma i}\left(\cos E_{i}(\mathbf{p})T+i\frac{\|\mathbf{p}\|}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\right)a^{\dagger}_{L\gamma}(\pm\mathbf{p},t_{0})\right.\\ &\left.\quad\mp V^{\ast}_{\alpha i}U_{\gamma i}\frac{m_{i}}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\,b_{R\gamma}(\mp\mathbf{p},t_{0})\right],\end{split}$		(61)
	$\displaystyle\begin{split}b_{L\alpha}(\pm\mathbf{p},t)=\sum_{i=1}^{3}\sum_{\beta=e}^{\tau}&\left[V^{\ast}_{\alpha i}V_{\beta i}\left(\cos E_{i}(\mathbf{p})T-i\frac{\|\mathbf{p}\|}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\right)b_{L\beta}(\pm\mathbf{p},t_{0})\right.\\ &\left.\quad\mp V^{\ast}_{\alpha i}U_{\beta i}\frac{m_{i}}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\,a^{\dagger}_{R\beta}(\mp\mathbf{p},t_{0})\right],\end{split}$		(62)
	$\displaystyle\begin{split}b^{\dagger}_{L\alpha}(\pm\mathbf{p},t)=\sum_{i=1}^{3}\sum_{\gamma=e}^{\tau}&\left[V_{\alpha i}V^{\ast}_{\gamma i}\left(\cos E_{i}(\mathbf{p})T+i\frac{\|\mathbf{p}\|}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\right)b^{\dagger}_{L\gamma}(\pm\mathbf{p},t_{0})\right.\\ &\left.\quad\mp V_{\alpha i}U^{\ast}_{\gamma i}\frac{m_{i}}{E_{i}(\mathbf{p})}\sin E_{i}(\mathbf{p})T\,a_{R\gamma}(\mp\mathbf{p},t_{0})\right],\end{split}$		(63)

where the right-handed operators are just replacements of the PMNS matrices $U\rightarrow V$ , $V\rightarrow U$ and the handedness $a_{(L,R)}\rightarrow a_{(R,L)}$ , $b_{(L,R)}\rightarrow b_{(R,L)}$ . We define the lepton family number density, for both the left- and right-handed Dirac neutrinos,

	$\displaystyle l^{L}_{\alpha}(t,\mathbf{x})=\,:\overline{\nu_{L\alpha}}(t,\mathbf{x})\gamma^{0}\nu_{L\alpha}(t,\mathbf{x}):\,,$		(64)
	$\displaystyle l^{R}_{\alpha}(t,\mathbf{x})=\,:\overline{\nu_{R\alpha}}(t,\mathbf{x})\gamma^{0}\nu_{R\alpha}(t,\mathbf{x}):\,.$		(65)

Then, the evolution of the lepton family number density for the left-handed case Eq.(64) is obtained by substituting the time dependent form of Eq.(43),

l^{L}_{\alpha}(t,\mathbf{x})=\int^{\prime}\frac{d^{3}\mathbf{k}}{(2\pi)^{3}2|\mathbf{k}|}\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2|\mathbf{p}|}\\ \times\left[a^{\dagger}_{L\alpha}(\mathbf{k},t)a_{L\alpha}(\mathbf{p},t)\overline{u_{L}}(\mathbf{k})\gamma^{0}u_{L}(\mathbf{p})e^{-i(\mathbf{k}-\mathbf{p})\cdot\mathbf{x}}+b_{L\alpha}(\mathbf{k},t)a_{L\alpha}(\mathbf{p},t)\overline{v_{L}}(\mathbf{k})\gamma^{0}u_{L}(\mathbf{p})e^{i(\mathbf{k}+\mathbf{p})\cdot\mathbf{x}}\right.\\ +\left.a^{\dagger}_{L\alpha}(\mathbf{k},t)b^{\dagger}_{L\alpha}(\mathbf{p},t)\overline{u_{L}}(\mathbf{k})\gamma^{0}v_{L}(\mathbf{p})e^{-i(\mathbf{k}+\mathbf{p})\cdot\mathbf{x}}-b^{\dagger}_{L\alpha}(\mathbf{p},t)b_{L\alpha}(\mathbf{k},t)\overline{v_{L}}(\mathbf{k})\gamma^{0}v_{L}(\mathbf{p})e^{i(\mathbf{k}-\mathbf{p})\cdot\mathbf{x}}\right],

(66)

whereas the right-handed form is a replacement of the spinors and operators from $L$ to $R$ . Next, we write the expectation value using a similar state $|\psi_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ as Eq.(32), but with the left-handed operator replaced as $a^{\dagger}_{\sigma}(\mathbf{q},t_{0})\rightarrow a^{\dagger}_{L\sigma}(\mathbf{q},t_{0})$ . We sandwich the left-handed density for the lepton family number of Eq.(66) with the left-handed state $|\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ and integrate over $\mathbf{k},\mathbf{p}$ to result in,

\langle\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\lvert l^{L}_{\alpha}(t,\mathbf{x})\rvert\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle=\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}\frac{dq^{\prime}dq}{A}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)\mathbf{e}_{2}\cdot\mathbf{x}}\\ \times\left[\sum_{i,j}V^{\ast}_{\alpha i}V_{\sigma i}V_{\alpha j}V^{\ast}_{\sigma j}\left(\cos{E_{i}(q^{\prime})T}+i\frac{|q^{\prime}|}{E_{i}(q^{\prime})}\sin{E_{i}(q^{\prime})T}\right)\left(\cos{E_{j}(q)T}-i\frac{|q|}{E_{j}(q)}\sin{E_{j}(q)T}\right)\right].

(67)

An initial right-handed state of Eq.(67), $\langle\psi^{R}_{\sigma}(q^{0};\sigma_{q}),t_{0}\lvert l^{R}_{\alpha}(t,\mathbf{x})\rvert\psi^{R}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ , is obtained by interchanging the matrices as $U\rightarrow V,\quad V\rightarrow U$ . Lastly, we perform the integration over $q$ and $q^{\prime}$ using the approximation of Eq.(37) to write the integrand in Gaussian form. The result is the linear density expectation value of the left-handed lepton family number,

\begin{split}\lambda^{L}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})\simeq&\frac{\sigma_{q}}{(2\pi)^{1/2}}\sum_{i,j}V^{\ast}_{\alpha i}V_{\sigma i}V_{\alpha j}V^{\ast}_{\sigma j}\\ &\times\frac{1}{2}\left[(v_{i0}+v_{j0}+1+v_{i0}v_{j0})e^{i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}-v_{i0}T)^{2}+(x_{2}-v_{j0}T)^{2}]}\right.\\ &\quad\quad-(v_{i0}+v_{j0}-1-v_{i0}v_{j0})e^{-i(E_{i}(q^{0})-E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}+v_{i0}T)^{2}+(x_{2}+v_{j0}T)^{2}]}\\ &\quad\quad+(v_{i0}-v_{j0}+1-v_{i0}v_{j0})e^{i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}-v_{i0}T)^{2}+(x_{2}+v_{j0}T)^{2}]}\\ &\left.\quad\quad-(v_{i0}-v_{j0}-1+v_{i0}v_{j0})e^{-i(E_{i}(q^{0})+E_{j}(q^{0}))T}e^{-\sigma^{2}_{q}[(x_{2}+v_{i0}T)^{2}+(x_{2}-v_{j0}T)^{2}]}\right].\end{split}

(68)

Where we have defined $\lambda^{L}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})=\iint dx_{1}dx_{3}\langle\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}|l^{L}_{\alpha}(t,\mathbf{x})|\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ . The process for solving for the linear density of the right-handed lepton family number of Eq.(65) is the same. After sandwiching the right-handed density of Eq.(65) with the left-handed state $|\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ and integrating over $\mathbf{k},\mathbf{p}$ our result is,

\langle\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}|l^{R}_{\alpha}(t,\mathbf{x})|\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle=\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}\frac{dq^{\prime}dq}{A}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)\mathbf{e}_{2}\cdot\mathbf{x}}\\ \times\left[\sum_{i,j}U_{\alpha j}V^{\ast}_{\sigma j}U^{\ast}_{\alpha i}V_{\sigma i}\frac{m_{j}}{E_{j}(q^{\prime})}\sin{E_{j}(q^{\prime})T}\frac{m_{i}}{E_{i}(q)}\sin{E_{i}(q)T}\right].

(69)

Thus, we write the resulting linear density expectation value of the right-handed lepton family number,

\lambda^{R}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})\simeq\frac{\sigma_{q}}{\sqrt{2\pi}}\sum_{i,j}\sqrt{1-v_{i0}^{2}}\sqrt{1-v_{j0}^{2}}e^{-\sigma_{q}^{2}\left((v_{i0}^{2}+v_{j0}^{2})T^{2}+2x_{2}^{2}\right)}\\ \times\left\{\text{Re}\left(U^{\ast}_{\alpha i}V_{\sigma i}U_{\alpha j}V^{\ast}_{\sigma j}\right)\left[\cosh\left(2\sigma_{q}^{2}(v_{i0}+v_{j0})Tx_{2}\right)\cos\left((E_{i}(q^{0})-E_{j}(q^{0}))T\right)\right.\right.\\ \left.-\cosh\left(2\sigma_{q}^{2}(v_{i0}-v_{j0}\right)Tx_{2})\cos\left((E_{i}(q^{0})+E_{j}(q^{0}))T\right)\right]\\ -\text{Im}\left(U^{\ast}_{\alpha i}V_{\sigma i}U_{\alpha j}V^{\ast}_{\sigma j}\right)\left[\sinh\left(2\sigma_{q}^{2}(v_{i0}+v_{j0})Tx_{2}\right)\sin\left((E_{i}(q^{0})-E_{j}(q^{0}))T\right)\right.\\ \left.\left.-\sinh\left(2\sigma_{q}^{2}(v_{i0}-v_{j0})Tx_{2}\right)\sin\left((E_{i}(q^{0})+E_{j}(q^{0}))T\right)\right]\right\},

(70)

where $\lambda^{R}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})=\iint dx_{1}dx_{3}\langle\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}|l^{R}_{\alpha}(t,\mathbf{x})|\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ .

III.3 Comparison to Wave packet formulations

The results we present in Eqs.(38, 68, 70) are based on the evolution of a density operator. We will now discuss an interpretation of those results based on wave packet formulations. Specifically, wave packets have been used in connection with neutrino oscillation phenomenology (see Ref. [41] for a review) as several forms; the quantum mechanical (QM) formulation [22], and the external wave packet QFT formulation [40, 38]. We will focus on comparisons to interpretations in the quantum mechanical formulation, because it is the simplest in literature.

Often discussed in literature see [23, 33, 24, 25], neutrino oscillations occur when three types of coherence are satisfied.

1.

The different massive neutrino components are coherently produced.
2.

Coherent propagation of massive neutrino components.
3.

Coherent detection of the different massive components.

We do not consider any production or detection processes, so we can not discuss the possibility of coherent detection in our formulation. The initial density state of Eq.(32) captures the properties of coherent propagation, similar to the QM wave packet.

A prediction of the QM wave packet formulation is the coherence length $L^{\text{coh}}_{i,j}$ for oscillations [26]. The coherence length is a measure of the group velocities for the mass eigenstates that appears after integration over time. It acts as a damping factor on the oscillations in space. We find a similar damping in our density formulation. The damping appears as a real exponential component that depends quadratically on spacetime,

	$\displaystyle e^{-\sigma^{2}_{q}[(x_{2}\pm v_{i0}T)^{2}+(x_{2}\pm v_{j0}T)^{2}]}$		(71)
	$\displaystyle e^{-\sigma^{2}_{q}[(x_{2}\pm v_{i0}T)^{2}+(x_{2}\mp v_{j0}T)^{2}]}.$		(72)

These real components are a non-linear damping that is applied to the density oscillation. An important distinction between our density formulation and the QM formulation are the types of damping. The QM formulation of Eq.(3) in Ref. [22] has a single type of damping to coincide with the single type of oscillation

e^{-\sigma^{2}_{q}[(x_{2}-v_{i0}T)^{2}+(x_{2}-v_{j0}T)^{2}]},

(73)

which we denote as $(-,-)$ . Whereas, our density formulation has four different types damping, $(-,-),(-,+),(+,-),(+,+)$ . Each damping is applied to a different oscillation in Eqs.(38, 68, 70). The $(+,+)$ is the strongest damping factor and $(-,-)$ is the weakest damping, which we illustrate by minimizing the polynomials inside the exponentials assuming $x_{2}>0$ ;

$\displaystyle P_{1}\left(x_{2},T=\frac{2x_{2}}{v_{i0}+v_{j0}}\right)$	$\displaystyle=2x_{2}^{2}\frac{(v_{i0}-v_{j0})^{2}}{(v_{i0}+v_{j0})^{2}}$	$\displaystyle\text{from }(-,-),$	(74)
$\displaystyle P_{2}\left(x_{2},T=\frac{2x_{2}}{v_{i0}+v_{j0}}\right)$	$\displaystyle=2x_{2}^{2}\frac{3(v_{i0}+v_{j0})^{2}+2(v_{i0}^{2}+v_{j0}^{2})}{(v_{i0}+v_{j0})^{2}}$	$\displaystyle\text{from }(+,+),$	(75)
$\displaystyle P_{3}\left(x_{2},T=\frac{2x_{2}}{v_{i0}+v_{j0}}\right)$	$\displaystyle=2x_{2}^{2}+8x_{2}^{2}\frac{v_{j0}^{2}}{(v_{i0}+v_{j0})^{2}}$	$\displaystyle\text{from }(-,+),$	(76)
$\displaystyle P_{4}\left(x_{2},T=\frac{2x_{2}}{v_{i0}+v_{j0}}\right)$	$\displaystyle=2x_{2}^{2}+8x_{2}^{2}\frac{v_{i0}^{2}}{(v_{i0}+v_{j0})^{2}}$	$\displaystyle\text{from }(+,-).$	(77)

Notice the peak value of $(-,-)$ is damped at the peak by $P_{1}<2x_{2}^{2}$ causing it to be weak. Whereas the peak values of $(+,+)$ , $(-,+)$ , and $(+,-)$ are damped with $P_{2,3,4}>2x_{2}^{2}$ strengthening their suppression of the oscillations.

If we perturb slightly away from the peak time, the strength of the damping factors depend directly on the velocities. For example,

P_{1}\left(x_{2},T=\frac{2x_{2}}{v_{i0}+v_{j0}}+\Delta T\right)-P_{1}\left(x_{2},T=\frac{2x_{2}}{v_{i0}+v_{j0}}\right)=(v_{i0}^{2}+v_{j0}^{2})(\Delta T)^{2}+\frac{2x_{2}(v_{i0}-v_{j0})^{2}}{v_{i0}+v_{j0}}\Delta T.

(78)

So, a decrease in the velocities lead to longer damping times.

After we integrate over all space to arrive at Eq.(39), the damping factors appear as two types $\pm$ .

T^{coh}_{i,j}\simeq\frac{\sqrt{2}}{\sigma_{q}(v_{i0}\pm v_{j0})},

(79)

which we write as a coherence time $T_{i,j}^{coh}$ . In the usual QM case only a single damping factor appears for the oscillations $(v_{i0}+v_{j0})$ . However, new oscillation terms appear in our formulation corresponding to $e^{\pm i[E_{i}(q^{0})+E_{j}(q^{0})]T}$ , to which a new damping factor $(v_{i0}-v_{j0})$ is applied.

IV Comparison of the Dirac and Majorana formulations

From the expectation value equations of Eq.(38), Eq.(68), and Eq.(70) we find signatures of the Dirac and Majorana formulations. Firstly, the terms proportional to the PMNS matrix combination of $V_{\alpha i}^{\ast}V_{\sigma i}V_{\alpha j}V_{\sigma j}^{\ast}$ are common between the formulations. Those terms are the first summation from Eq.(38) of the Majorana formulation and the left-handed Dirac formulation Eq.(68). The difference between the formulations is the secondary term of the Majorana formulation and the right-handed Dirac formulation Eq.(70). The secondary term of the Majorana formulation has the PMNS combination $V^{\ast}_{\alpha i}V^{\ast}_{\sigma i}V_{\alpha j}V_{\sigma j}$ that is related to the Majorana phases and is subtracted from the common oscillation terms. Whereas, the second unitary matrix $U_{\alpha i}$ from the right-handed Dirac formulation is present in Eq.(70) and is added to Eq.(68) as $\lambda^{D}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})=\iint dx_{1}dx_{3}\langle\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}|l^{L}_{\alpha}(t,\mathbf{x})+l^{R}_{\alpha}(t,\mathbf{x})|\psi^{L}_{\sigma}(q^{0};\sigma_{q}),t_{0}\rangle$ . The sign of the different terms, subtraction for Majorana and addition for Dirac, is responsible for total lepton number violation in the Majorana case. To help emphasize the affect of the sign difference we write the expectation value of the total lepton number for the Dirac formulation case,

\sum_{\alpha}\lambda^{D}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})\simeq\frac{\sigma_{q}}{\sqrt{2\pi}}\sum_{i}|V_{\sigma i}|^{2}\left[(1+v_{i0})e^{-2\sigma_{q}^{2}(x_{2}-v_{i0}T)^{2}}+(1-v_{i0})e^{-2\sigma_{q}^{2}(x_{2}+v_{i0}T)^{2}}\right],

(80)

where integration over $x_{2}$ would result in total lepton number conservation for all times. Next, from Eq.(38), we also take the summation over the lepton family numbers $\alpha$ as follows,

\begin{split}\sum_{\alpha}\lambda^{M}_{\sigma\rightarrow\alpha}(T=t-t_{0},x_{2})\simeq&\frac{\sigma_{q}}{\sqrt{2\pi}}\sum_{i}|V_{\sigma i}|^{2}\left[v_{i0}(1+v_{i0})e^{-2\sigma_{q}^{2}(x_{2}-v_{i0}T)^{2}}-v_{i0}(1-v_{i0})e^{-2\sigma_{q}^{2}(x_{2}+v_{i0}T)^{2}}\right.\\ &\left.\hphantom{XXXXXXXXXX}+2(1-v_{i0}^{2})e^{-2\sigma_{q}^{2}(v_{i0}^{2}T^{2}+x_{2}^{2})}\cos 2E_{i}(q^{0})T\right].\end{split}

(81)

In Eq.(80) for the Dirac case, there is no oscillation term and the linear density of the total lepton number is always positive. For a nonrelativistic Majorana neutrino ( $v_{i0}\ll 1$ ) in Eq.(81), the oscillation term can dominate and the sign of the linear density can change with respect to time.

V Illustration of nonrelativistic differences

From inspection of the expectation values Eq.(38), Eq.(68), and Eq.(70); and the discussions in sections III.3 and IV we see two key results.

•

Differences between our formulation and the usual QM formulation are negligible under the ultra-relativistic assumption.
•

The neutrino masses $\frac{m_{i}m_{j}}{E_{i}E_{j}}=\sqrt{1-v_{i}^{2}}\sqrt{1-v_{j}^{2}}$ suppress any modifications by the Dirac or Majorana mass formulations.

To illustrate those results we choose $m_{\text{lightest}}=0.01\text{ eV}<q^{0}=0.2\text{ eV}$ for the linear density expectations of Eqs.(38) and (68). Even for a momentum an order of magnitude greater than the mass, differences between the Majorana and Dirac cases are not visible (Fig. 1 and Fig. 2).

Refer to caption — Figure 1: Time evolution of the linear density for the expectation value of a lepton family number with Majorana, panel (a), or Dirac, panel (b), mass. We take an arbitrary distance slice at $x_{2}=2.5$ cm. The initial momentum density is a Gaussian distribution with a width of $\sigma_{q}=0.00001$ and a mean momentum of $q^{0}=0.2$ eV. Normal mass hierarchy is considered. Oscillation parameters are the best fit values from the NuFIT 5.0 (2020) collaboration [54].

We used the best fit values of the PMNS mixing angles $\theta_{12},\theta_{23},\text{and }\theta_{13}$ , the CP violating phase $\delta_{\text{CP}}$ , and the mass squared differences $\Delta m_{21}^{2}\text{ and }\Delta m_{31}^{2}$ from the work of the NuFIT 5.0 (2020) collaboration [54]. Then, we are left to choose the absolute mass of the lightest neutrino $m_{1}\text{ or }m_{3}$ , the neutrino mass hierarchy, normal $m_{1}<m_{2}\ll m_{3}$ or inverted $m_{3}\ll m_{1}<m_{2}$ , the value of the mean momentum $q$ , the width of the initial Gaussian density $\sigma_{q}$ , the time $T=t-t_{0}$ , the distance $x_{2}$ , and the Majorana phases $\alpha_{21}\text{ and }\alpha_{31}$ . The definition of the Majorana phases is given as $\alpha_{21}=2\arg(V_{e2})$ and $\alpha_{31}=2\arg(V_{\mu 3})$ . Any changes caused by different Majorana phases are also hidden, because the differences between the Majorana and Dirac cases are not visible for $m_{\text{lightest}}=0.01\text{ eV}<q^{0}=0.2\text{ eV}$ .

The decoherence effects discussed in section III.3 from the real exponentials of Eq.(38), Eq.(68), and Eq.(70) cause the oscillations to be localized to a spacetime region. In Fig. 3, we show the 2-D spacetime contour for an electron number linear density for the $\sigma=e\rightarrow\alpha=e$ case. The region where the electron number linear density has a maximum peak of $\sim 0.4$ propagates from the spacetime origin $(t,x_{2})=(0,0)$ toward the upper right corner at a constant velocity. The initial Gaussian shape is maintained during propagation because we assumed no wave packet spreading in the approximation of Eq.(37).

We discussed in section IV how the total lepton number becomes violated in the Majorana case, and we illustrate that by decreasing the momentum to $q^{0}=0.0002<m_{lightest}=0.01$ eV in Fig.4. This leads to the distinction between the Majorana and Dirac cases, in which the Majorana case has negative expectation values. Both cases feature the small period oscillations of the inset graph because of lines 4 and 5 from Eq.(68). Those lines are not present in the usual QM formulation, so these small period oscillations can distinguish our formulation.

By comparing Fig.1 and Fig.4, the oscillations occur for longer times when the average momentum is smaller. This was discussed in section III.3, where Eq.(78) was given as an example. To recall, the strength of the damping factors depend directly on the velocities. So, a decrease in the velocities lead to longer damping times.

Next, in the distance plane of Fig.5 the smaller momenta suppress the peak value for the linear density of the lepton family numbers. This follows directly from Eq.(74), where for $q^{0}=0.2$ eV the polynomial is larger than for $q^{0}=0.0002$ eV. Thus, the propagation coherence occurs at smaller distances, which is sometimes stated as a smaller coherence length.

Even with this in mind, at smaller momenta the Dirac and Majorana cases are distinguishable. Similar to Fig.4, the Majorana case has negative expectation values and the Dirac case is always positive.

VI Conclusion

We have proposed a novel formulation for neutrino oscillations from a QFT point of view. This formulation can be applied to both relativistic and nonrelativistic energies for neutrinos. The previous formulation was limited to Majorana neutrinos with a fixed momentum. We have extended the previous formulation to the Dirac neutrino case and have applied to the case with a momentum distribution. To study the oscillation of neutrino with a momentum distribution, we have considered the linear lepton number density and the initial state with a Gaussian distribution along a one spacial direction. After taking the expectation value of the density with the state, we find a result similar to the quantum mechanical wave packet approach used to describe neutrino flavor oscillations. Our result has additional terms that result in interesting behavior. In particular, we have studied type of decoherence ascribed to wave packet separation. We have found how the peak value of lepton number density is suppressed as it propagates. Even with that decoherence our formulation can distinguish between the neutrino mass type at nonrelativistic energies. The mass types, the expectation value of the lepton number density for Majorana neutrinos can take both positive and negative value while for Dirac neutrinos, it does not change the sign.

Acknowledgements.

The work of T.M. is supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Number JP17K05418. We would like to thank Uma Sankar and Koichi Hamaguchi for useful comments and suggestions during the initial presentations of this work. We would like to thank Tomohiro Inagaki, Tomoharu Orimi and Naoki Uemura for useful discussions. Lastly, we thank the three anonymous referees for the detailed review and improvements.

Appendix A Quantization with the zero mode and the massless spinors

In this appendix, we demonstrate the quantization with the massless spinor in Eq.(8) by adding the zero mode. The zero mode corresponds a massive fermion in its rest frame. The quantization based on the Dirac bracket is carried out. For simplicity, we consider only the single flavor version of the Lagrangian in Eq.(6). We write the Lagrangian in terms of two component spinor $\eta$ including the zero mode as will be shown in Eq.(94),

\nu_{L}(x,t)=\left(\begin{matrix}0\\ \eta(x,t)\end{matrix}\right).

(82)

Then the Lagrangian is given as follows,

\mathcal{L}=\eta^{\dagger}i\overline{\sigma}^{\mu}\partial_{\mu}\eta-\frac{m}{2}\left(-\eta^{\dagger}i\sigma_{2}\eta^{\ast}+\eta^{T}i\sigma_{2}\eta\right).

(83)

By treating $\eta^{\dagger}$ as independent variables,

\displaystyle\pi_{\eta}=\frac{\partial\mathcal{L}}{\partial\dot{\eta}}=i\eta^{\dagger},

\displaystyle\pi_{\eta^{\dagger}}=\frac{\partial\mathcal{L}}{\partial{\dot{\eta}^{\dagger}}}=0,

(84)

one obtains the following two constraints,

\displaystyle\phi_{1}=\pi_{\eta}-i\eta^{\dagger}=0,

\displaystyle\phi_{2}=\pi_{\eta^{\dagger}}=0.

(85)

The Hamiltonian is given by,

H=i\int dx\,\eta^{\dagger}\bm{\sigma}\cdot\bm{\nabla}\eta+\frac{m}{2}\int dx\left(-\eta^{\dagger}i\sigma_{2}\eta^{\ast}+\eta^{T}i\sigma_{2}\eta\right).

(86)

Following Dirac, we aim to obtain the anti-commutation relation, for $\eta$ , based on Dirac bracket. For instance, the Dirac bracket for $\eta$ and $\eta^{\dagger}$ is calculated to be,

\{\eta(\mathbf{x},t),\eta^{\dagger}(\mathbf{y},t)\}_{\text{DB}}=\{\eta(\mathbf{x},t),\eta^{\dagger}(\mathbf{y},t)\}_{\text{PB}}\\ -\sum_{a=1,2}\int d^{3}x^{\prime}d^{3}y^{\prime}\{\eta(\mathbf{x},t),\phi_{a}(\mathbf{x}^{\prime},t)\}_{\text{PB}}(N^{-1})_{ab}(\mathbf{x}^{\prime},\mathbf{y}^{\prime})\{\phi_{b}(\mathbf{y}^{\prime},t),\eta^{\dagger}(\mathbf{y},t)\}_{\text{PB}},

(87)

where PB denotes the Poisson bracket and DB the Dirac bracket. $N_{ab}$ is a matrix of the Poisson bracket for the constraints,

N_{ab}(\mathbf{x},\mathbf{y})=\{\phi_{a}(\mathbf{x},t),\phi_{b}(\mathbf{y},t)\}_{\text{PB}}=\begin{pmatrix}0&-i\\ -i&0\end{pmatrix}\delta^{(3)}(\mathbf{x}-\mathbf{y})\mathbbm{1}_{2\times 2},

(88)

where $\mathbbm{1}_{2\times 2}$ represents a two by two unit matrix that acts on the two component spinors. The inverse matrix $(N^{-1})_{ab}$ is calculated to be,

(N^{-1})_{ab}=\begin{pmatrix}0&i\\ i&0\end{pmatrix}\delta^{(3)}(\mathbf{x}-\mathbf{y})\mathbbm{1}_{2\times 2}.

(89)

Then the Dirac brackets among $\eta$ and $\eta^{\dagger}$ are given as,

	$\displaystyle\{\eta(\mathbf{x},t),\eta^{\dagger}(\mathbf{y},t)\}_{\text{DB}}=-i\delta^{(3)}(\mathbf{x}-\mathbf{y})\mathbbm{1}_{2\times 2},$		(90)
	$\displaystyle\{\eta(\mathbf{x},t),\eta(\mathbf{y},t)\}_{\text{DB}}=\{\eta^{\dagger}(\mathbf{x},t),\eta^{\dagger}(\mathbf{y},t)\}_{\text{DB}}=0.$		(91)

Then the anti-commutation relations are,

	$\displaystyle\{\eta(\mathbf{x},t),\eta^{\dagger}(\mathbf{y},t)\}=\delta^{(3)}(\mathbf{x}-\mathbf{y})\mathbbm{1}_{2\times 2},$		(92)
	$\displaystyle\{\eta(\mathbf{x},t),\eta(\mathbf{y},t)\}=\{\eta^{\dagger}(\mathbf{x},t),\eta^{\dagger}(\mathbf{y},t)\}=0.$		(93)

Now, we write the two component spinor $\eta$ with the massless, non-zero mode expansion from Eq.(8) and include a zero mode contribution,

\eta(\mathbf{x},t)=\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}\sqrt{2\lvert\mathbf{p}\rvert}}\left(a(\mathbf{p},t)\phi_{-}(\mathbf{n}_{\mathbf{p}})e^{i\mathbf{p}\cdot\mathbf{x}}-b^{\dagger}(\mathbf{p},t)\phi_{-}(\mathbf{n}_{\mathbf{p}})e^{-i\mathbf{p}\cdot\mathbf{x}}\right)+\eta^{0}(t).

(94)

In the right-hand side of Eq.(94), the first term represents non-zero modes while the second term represents zero mode defined as,

\int d^{3}x\eta(\mathbf{x},t)=V\eta^{0}(t),

(95)

with $V=(2\pi)^{3}\delta^{(3)}(\mathbf{p}=0)$ . Our goal is to determine the equal-time operator relations among the creation and annihilation operators for both zero mode and non-zero modes with the knowledge of Eq.(92) and Eq.(93). Additionally, we will introduce the zero mode contribution to the lepton number density operator and its expectation value. The anti-commutation for the zero mode is extracted with Eq.(92) by integrating over $\mathbf{x}$ and $\mathbf{y}$ ,

\{\eta^{0}_{i}(t),\eta^{0\dagger}_{j}(t)\}=\frac{1}{V}\delta_{ij},

(96)

where we denote the two component spinor indices with $i,j=1,2.$ Next we focus on the non-zero modes. Recall the meaning of $\int^{\prime}$ , in Eq.(94), is from Eq.(8) and results in two momenta regions for the spinors. We use that to explicitly rewrite the integration as a single momenta region, which we label $A$ ,

\begin{split}\eta(\mathbf{x},t)=\eta^{0}(t)+\int_{\mathbf{p}\in A}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}\sqrt{2\lvert\mathbf{p}\rvert}}&\left\{[a(\mathbf{p},t)\phi_{-}(\mathbf{n}_{\mathbf{p}})-b^{\dagger}(-\mathbf{p},t)\phi_{-}(-\mathbf{n}_{\mathbf{p}})]e^{i\mathbf{p}\cdot\mathbf{x}}\right.\\ &\left.+[a(-\mathbf{p},t)\phi_{-}(-\mathbf{n}_{\mathbf{p}})-b^{\dagger}(\mathbf{p},t)\phi_{-}(\mathbf{n}_{\mathbf{p}})]e^{-i\mathbf{p}\cdot\mathbf{x}}\right\}.\end{split}

(97)

Then, the Fourier transformation is carried out,

	$\displaystyle\int d^{3}\mathbf{x}\,\eta(\mathbf{x},t)e^{-i\mathbf{q}\cdot\mathbf{x}}=\frac{1}{\sqrt{2\|\mathbf{q}\|}}\left[a(\mathbf{q},t)\phi_{-}(\mathbf{n}_{\mathbf{q}})-b^{\dagger}(-\mathbf{q},t)\phi_{-}(-\mathbf{n}_{\mathbf{q}})\right],$		(98)
	$\displaystyle\int d^{3}\mathbf{x}\,\eta(\mathbf{x},t)e^{i\mathbf{q}\cdot\mathbf{x}}=\frac{1}{\sqrt{2\|\mathbf{q}\|}}\left[a(-\mathbf{q},t)\phi_{-}(-\mathbf{n}_{\mathbf{q}})-b^{\dagger}(\mathbf{q},t)\phi_{-}(\mathbf{n}_{\mathbf{q}})\right],$		(99)

where $\mathbf{q}$ is a momentum in $A$ region. This leads to expressions for the operators of non-zero $\mathbf{q}$ ,

	$\displaystyle\int d^{3}\mathbf{x}\phi^{\dagger}_{-}(\mathbf{n}_{\mathbf{q}})\eta(\mathbf{x},t)e^{-i\mathbf{q}\cdot\mathbf{x}}=\frac{1}{\sqrt{2\|\mathbf{q}\|}}a(\mathbf{q},t),\int d^{3}\mathbf{x}\phi^{\dagger}_{-}(-\mathbf{n}_{\mathbf{q}})\eta(\mathbf{x},t)e^{i\mathbf{q}\cdot\mathbf{x}}=\frac{1}{\sqrt{2\|\mathbf{q}\|}}a(-\mathbf{q},t),$		(100)
	$\displaystyle\int d^{3}\mathbf{x}\phi^{\dagger}_{-}(\mathbf{n}_{\mathbf{q}})\eta(\mathbf{x},t)e^{i\mathbf{q}\cdot\mathbf{x}}=-\frac{1}{\sqrt{2\|\mathbf{q}\|}}b^{\dagger}(\mathbf{q},t),\int d^{3}\mathbf{x}\phi^{\dagger}_{-}(-\mathbf{n}_{\mathbf{q}})\eta(\mathbf{x},t)e^{-i\mathbf{q}\cdot\mathbf{x}}=-\frac{1}{\sqrt{2\|\mathbf{q}\|}}b^{\dagger}(-\mathbf{q},t).$		(101)

Similar expressions are found for the Hermitian conjugate $\eta^{\dagger}(t,\mathbf{x})$ ,

	$\displaystyle\int d^{3}\mathbf{x}\,\eta^{\dagger}(\mathbf{x},t)\phi_{-}(\mathbf{n}_{\mathbf{q}})e^{i\mathbf{q}\cdot\mathbf{x}}=\frac{1}{\sqrt{2\|\mathbf{q}\|}}a^{\dagger}(\mathbf{q},t),\int d^{3}\mathbf{x}\,\eta^{\dagger}(\mathbf{x},t)\phi_{-}(-\mathbf{n}_{\mathbf{q}})e^{-i\mathbf{q}\cdot\mathbf{x}}=\frac{1}{\sqrt{2\|\mathbf{q}\|}}a^{\dagger}(-\mathbf{q},t),$		(102)
	$\displaystyle\int d^{3}\mathbf{x}\,\eta^{\dagger}(\mathbf{x},t)\phi_{-}(\mathbf{n}_{\mathbf{q}})e^{-i\mathbf{q}\cdot\mathbf{x}}=-\frac{1}{\sqrt{2\|\mathbf{q}\|}}b(\mathbf{q},t),\int d^{3}\mathbf{x}\,\eta^{\dagger}(\mathbf{x},t)\phi_{-}(-\mathbf{n}_{\mathbf{q}})e^{i\mathbf{q}\cdot\mathbf{x}}=-\frac{1}{\sqrt{2\|\mathbf{q}\|}}b(-\mathbf{q},t).$		(103)

Now we can calculate the equal-time operator relations for $a(\mathbf{p},t)$ and $b(\mathbf{p},t)$ ,

	$\displaystyle\{a(\mathbf{q},t),a^{\dagger}(\mathbf{p},t)\}=(2\pi)^{3}2\lvert\mathbf{p}\rvert\delta^{(3)}(\mathbf{p}-\mathbf{q}),$		(104)
	$\displaystyle\{b(\mathbf{q},t),b^{\dagger}(\mathbf{p},t)\}=(2\pi)^{3}2\lvert\mathbf{p}\rvert\delta^{(3)}(\mathbf{p}-\mathbf{q}),$		(105)

with all other relations being zero. In addition to these anti-commutation relation for non-zero modes, the relation for zero-mode is given in Eq.(96). Notice the non-zero mode relations match Eq.(24), which were calculated using the Bogoliubov transformation. We rewrite the Hamiltonian in Eq.(86) with the creation and annihilation operators. After the normal ordering, the result is

H=H_{0}+\int^{\prime}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2|\mathbf{p}|}|\mathbf{p}|[a^{\dagger}(\mathbf{p},t)a(\mathbf{p},t)+b^{\dagger}(\mathbf{p},t)b(\mathbf{p},t)]\\ +m\int_{\mathbf{p}\in A}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}2|\mathbf{p}|}[-ia(\mathbf{p},t)a(-\mathbf{p},t)-ib(\mathbf{p},t)b(-\mathbf{p},t)+\text{h.c.}]

(106)

where $H_{0}$ is the Hamiltonian for the zero mode and it is defined as,

H_{0}=\frac{m}{2}[\eta^{0T}(t)i\sigma_{2}\eta^{0}(t)-\eta^{0\dagger}(t)i\sigma_{2}\eta^{0\ast}(t)]V.

(107)

The time evolution of the zero mode reads,

	$\displaystyle i\frac{d}{dx^{0}}\eta^{0}(t)=[\eta^{0},H]=m(-i\sigma_{2}\eta^{0\dagger}(t))$		(108)
	$\displaystyle i\frac{d}{dx^{0}}(-i\sigma_{2}\eta^{0\dagger}(t))=m\eta^{0}(t)$		(109)

The solution to both Eq.(108) and Eq.(109) is given by,

\eta^{0}(t)=\cos m(t-t_{0})\eta^{0}(t_{0})-i\sin m(t-t_{0})(-i\sigma_{2}\eta^{0\dagger}(t_{0})).

(110)

For the non-zero modes, considering $\mathbf{p}\in A$ , the equations of motion are

	$\displaystyle i\frac{d}{dt}a(\mathbf{p},t)=[a(\mathbf{p},t),H]=\|\mathbf{p}\|a(\mathbf{p},t)-ima^{\dagger}(-\mathbf{p},t),$		(111)
	$\displaystyle i\frac{d}{dt}a^{\dagger}(-\mathbf{p},t)=[a^{\dagger}(-\mathbf{p},t),H]=-\|\mathbf{p}\|a^{\dagger}(-\mathbf{p},t)+ima(\mathbf{p},t).$		(112)

These lead to solutions consistent with Eq.(29) for the multi flavor case,

	$\displaystyle a(\mathbf{p},t)=\left[\cos E(\mathbf{p})T-i\frac{\|\mathbf{p}\|}{E(\mathbf{p})}\sin E(\mathbf{p})T\right]a(\mathbf{p},t_{0})-\frac{m}{E(\mathbf{p})}\sin[E(\mathbf{p})T]a^{\dagger}(-\mathbf{p},t_{0})$		(113)
	$\displaystyle a^{\dagger}(-\mathbf{p},t)=\left[\cos E(\mathbf{p})T+i\frac{\|\mathbf{p}\|}{E(\mathbf{p})}\sin E(\mathbf{p})T\right]a^{\dagger}(-\mathbf{p},t_{0})+\frac{m}{E(\mathbf{p})}\sin[E(\mathbf{p})T]a(\mathbf{p},t_{0}).$		(114)

We used the definition $T\equiv t-t_{0}$ . The relations consistent with Eq.(30) hold for $b(\mathbf{p},t)$ and $b^{\dagger}(-\mathbf{p},t)$ . A difference between the zero-mode solution of Eq.(110) and non-zero-modes of Eq.(113) is the amplitude factors. For the zero-mode cosine and sine amplitudes are not modified by $\lvert\mathbf{p}\rvert/E(\mathbf{p})$ or $m/E(\mathbf{p})$ .

Next we study the contribution to the lepton number density including the zero-mode. The zero mode is represented with $\eta^{0}(t)$ and its charge conjugation $-i\sigma_{2}\eta^{0\dagger}(t)$ as shown in Eq.(97). We relate the zero mode to a massive Majorana field. The Fourier expansion of the massive Majorana field is

\psi_{M}(\mathbf{x},t)=\int\frac{d^{3}p}{(2\pi)^{3}2E(\mathbf{p})}\sum_{s=\pm 1}(a_{M}(\mathbf{p},s)u(\mathbf{p},s)e^{-ip\cdot x}+a_{M}^{\dagger}(\mathbf{p},s)v(\mathbf{p},s)e^{ip\cdot x}),

(115)

where $s$ is related to the spin component in an arbitrary direction at the rest frame. To extract the zero mode from that expansion, we integrate over all space

\psi_{M}(t)=\frac{\int d^{3}x\psi_{M}(\mathbf{x},t)}{V}=\frac{1}{2mV}\sum_{s=\pm 1}(a_{M}(s)u(s)e^{-imt}+a_{M}^{\dagger}(s)v(s)e^{imt})

(116)

where the operator is $a_{M}(s)\equiv a_{M}(\mathbf{p}=0,s)$ . We have defined $u(s)\equiv u(\mathbf{p}=0,s)$ and $v(s)\equiv v(\mathbf{p}=0,s)$ as the positive and negative energy solutions for the particle at rest. In spinor/chiral representation of the gamma matrices, they are given as

\displaystyle u(s)=\sqrt{m}\begin{pmatrix}\chi^{(s)}\\ \chi^{(s)}\end{pmatrix},

\displaystyle v(s)=i\gamma^{2}u(s)^{\ast}=\sqrt{m}\begin{pmatrix}-i\sigma_{2}\chi^{(s)\ast}\\ i\sigma_{2}\chi^{(s)\ast}\end{pmatrix}.

(117)

One can choose the two component spinors $\chi^{(s)}(s=\pm 1)$ as spin up or down with respect to an arbitrary direction. Thus, they satisfy

(\mathbf{n}\cdot\sigma)\chi^{(s=\pm)}=\pm\chi^{(s=\pm)}.

(118)

The zero mode of the expansion in Eq.(116) should be matched to the zero mode of Eq.(82) through the relation

\psi_{M}(\mathbf{x},t)=\nu_{L}(\mathbf{x},t)+(\nu_{L}(\mathbf{x},t))^{c}.

(119)

That relation is equivalent to Eq.(19). We compare the zero modes of both side of equation to obtain,

\psi_{M}(t)=\begin{pmatrix}-i\sigma_{2}\eta^{0\dagger}(t)\\ \eta^{0}(t)\end{pmatrix}.

(120)

To facilitate the comparison, we rewrite the right-hand side of Eq.(120) using Eq.(110),

\begin{pmatrix}-i\sigma_{2}\eta^{0\dagger}(t)\\ \eta^{0}(t)\end{pmatrix}=\frac{e^{-im(t-t_{0})}}{2}\begin{pmatrix}\eta^{0}(t_{0})-i\sigma_{2}\eta^{0\dagger}(t_{0})\\ \eta^{0}(t_{0})-i\sigma_{2}\eta^{0\dagger}(t_{0})\end{pmatrix}+\frac{e^{im(t-t_{0})}}{2}\begin{pmatrix}-(\eta^{0}(t_{0})+i\sigma_{2}\eta^{0\dagger}(t_{0}))\\ \eta^{0}(t_{0})+i\sigma_{2}\eta^{0\dagger}(t_{0})\end{pmatrix}.

(121)

The left-hand side of Eq.(120) is given by Eq.(116)

\psi_{M}(t)=\frac{1}{2\sqrt{m}V}\sum_{s=\pm 1}\left[a_{M}(s)\begin{pmatrix}\chi^{(s)}\\ \chi^{(s)}\end{pmatrix}e^{-imt}+a_{M}^{\dagger}(s)\begin{pmatrix}-i\sigma_{2}\chi^{(s)\ast}\\ i\sigma_{2}\chi^{(s)\ast}\end{pmatrix}e^{imt}\right].

(122)

Comparing both sides of Eq.(120) given by Eqs.(121) and (122), one obtains,

	$\displaystyle\eta^{0}(t_{0})$	$\displaystyle=$	$\displaystyle\frac{1}{2\sqrt{m}V}\sum_{s=\pm 1}\left[a_{M}(s)\chi^{(s)}e^{-imt_{0}}+a^{\dagger}_{M}(s)i\sigma_{2}\chi^{(s)\ast}e^{imt_{0}}\right],$		(123)
		$\displaystyle=$	$\displaystyle\frac{1}{2\sqrt{m}V}\sum_{s=\pm 1}\left[a_{M}(s)\chi^{(s)}e^{-imt_{0}}-sa^{\dagger}_{M}(s)\chi^{(-s)}e^{imt_{0}}\right].$		(123)

where we use the relation $i\sigma_{2}\chi^{(s)\ast}=-s\chi^{(-s)}$ . By substituting the two component spinors $\chi^{(s)}$ explicitly, we rewrite Eq.(123)

\eta^{0}(t_{0})=\sum_{s=\pm 1}C(s,t_{0})\chi^{(s)},

(124)

where the new operator $C(s,t)$ is given by

C(s,t)=\frac{1}{\sqrt{2m}V}\frac{a_{M}(s)e^{-imt}+sa_{M}^{\dagger}(-s)e^{imt}}{\sqrt{2}}.

(125)

That obeys the anti-commutation relation following Eq.(96),

\{C(s,t),C^{\dagger}(s^{\prime},t)\}=\frac{1}{V}\delta_{ss^{\prime}}.

(126)

Additionally, it satisfies the differential equation,

i\frac{dC(s,t)}{dt}=-smC^{\dagger}(-s,t).

(127)

The solution is given by the operator defined at $t=t_{0}$ as,

C(s,t)=\cos m(t-t_{0})C(s,t_{0})+is\sin m(t-t_{0})C^{\dagger}(-s,t_{0})

(128)

This completes the expansion of Eq.(97) in terms of the operators for the non-zero modes $a(\pm\mathbf{p},t)$ and $b(\pm\mathbf{p},t)$ , and the zero mode operator $C(s,t)$ . The result is,

\begin{split}\eta(\mathbf{x},t)=\sum_{s=\pm 1}C(s,t)\chi^{(s)}+\int_{\mathbf{p}\in A}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}\sqrt{2\lvert\mathbf{p}\rvert}}&\left\{[a(\mathbf{p},t)\phi_{-}(\mathbf{n}_{\mathbf{p}})-b^{\dagger}(-\mathbf{p},t)\phi_{-}(-\mathbf{n}_{\mathbf{p}})]e^{i\mathbf{p}\cdot\mathbf{x}}\right.\\ &\left.+[a(-\mathbf{p},t)\phi_{-}(-\mathbf{n}_{\mathbf{p}})-b^{\dagger}(\mathbf{p},t)\phi_{-}(\mathbf{n}_{\mathbf{p}})]e^{-i\mathbf{p}\cdot\mathbf{x}}\right\}.\end{split}

(129)

Next, we consider the lepton number density operator including the zero mode contribution. This can be view as an extension of Eq.(31) for a single flavor,

l^{M}(t,\mathbf{x})=:\overline{\nu_{L}}(\mathbf{x},t)\gamma^{0}{\nu_{L}}(\mathbf{x},t):=\sum_{s=\pm 1}C^{\dagger}(s,t)C(s,t)\\ +\int^{\prime}\frac{d^{3}p}{(2\pi)^{3}\sqrt{2|\mathbf{p}|}}\sum_{s=\pm 1}C^{\dagger}(s,t)(a(\mathbf{p},t)e^{i\mathbf{p}\cdot\mathbf{x}}-b^{\dagger}(\mathbf{p},t)e^{-i\mathbf{p}\cdot\mathbf{x}})(\chi^{(s)\dagger}\cdot\phi_{-}(n_{\mathbf{p}}))\\ +\int^{\prime}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}\sum_{s=\pm 1}(\phi(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s)})(a^{\dagger}(\mathbf{k},t)e^{-i\mathbf{k}\cdot\mathbf{x}}-b(\mathbf{k},t)e^{i\mathbf{k}\cdot\mathbf{x}})C(s,t)\\ +\int^{\prime}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}\int^{\prime}\frac{d^{3}p}{(2\pi)^{3}\sqrt{2|\mathbf{p}|}}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\phi_{-}(n_{\mathbf{p}}))\\ \times:(a^{\dagger}(\mathbf{k},t)e^{-i\mathbf{k}\cdot\mathbf{x}}-b(\mathbf{k},t)e^{i\mathbf{k}\cdot\mathbf{x}})(a(\mathbf{p},t)e^{i\mathbf{p}\cdot\mathbf{x}}-b^{\dagger}(\mathbf{p},t)e^{-i\mathbf{p}\cdot\mathbf{x}}):

(130)

The first term of Eq.(130) comes from the zero mode. It is given by the operator defined at $t=t_{0}$ in Eq.(128) to be,

\sum_{s=\pm 1}C^{\dagger}(s,t)C(s,t)=\frac{1}{V}(1-\cos 2m(t-t_{0}))+\cos 2m(t-t_{0})\sum_{s=\pm 1}C^{\dagger}(s,t_{0})C(s,t_{0})\\ -i\{C(-1,t_{0})C(1,t_{0})-C^{\dagger}(1,t_{0})C^{\dagger}(-1,t_{0})\}\sin 2m(t-t_{0}).

(131)

The second and third lines of Eq.(130) are more complicated, because they are a mixture of zero with non-zero modes. Line four of Eq.(130) is equivalent to a single flavor version of Eq.(31). We take the expectation value of the lepton number density with the following momentum distribution

|\Psi(q^{0},\sigma_{q},c^{0}(s_{0})),t_{0}\rangle=|\psi(q^{0};\sigma_{q}),t_{0}\rangle+c^{0}(s_{0})C^{\dagger}(s_{0},t_{0})|0(t_{0})\rangle

(132)

where $|\psi(q^{0};\sigma_{q}),t_{0}\rangle$ is similar to the initial Gaussian distribution of Eq.(32). It is given by a superposition of the non-zero modes,

\lvert\psi(q^{0};\sigma_{q}),t_{0}\rangle=\frac{1}{\sqrt{\sigma_{q}}(2\pi)^{3/4}}\int^{\prime}\frac{dq}{\sqrt{A}\sqrt{2|q|}}e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}}a^{\dagger}(\mathbf{q},t_{0})\lvert 0(t_{0})\rangle;

(133)

where $\mathbf{q}=(0,q,0)$ . In addition to non-zero mode contribution, we add the zero mode contribution with the amplitude $c^{0}(s)$ . Note that the dimension of the coefficient is $\sqrt{V}$ . The expectation value of the lepton number densities given by the sum of the following terms

\langle\Psi(q^{0},\sigma_{q},c^{0}(s_{0})),t_{0}\rvert l^{M}(t,\mathbf{x})\lvert\Psi(q^{0},\sigma_{q},c^{0}(s_{0})),t_{0}\rangle=\langle\psi(q^{0};\sigma_{q}),t_{0}\rvert l^{M}(t,\mathbf{x})\lvert\psi(q^{0};\sigma_{q}),t_{0}\rangle\\ \begin{aligned} &+c^{0}(s_{0})\langle\psi(q^{0};\sigma_{q}),t_{0}\rvert l^{M}(t,\mathbf{x})C^{\dagger}(s_{0},t_{0})\lvert 0(t_{0})\rangle+\text{h.c.}\\ &+|c^{0}(s_{0})|^{2}\langle 0(t_{0})\rvert C(s_{0},t_{0})l^{M}(t,\mathbf{x})C^{\dagger}(s_{0},t_{0})\lvert 0(t_{0})\rangle.\end{aligned}

(134)

We start with the first line by substituting Eq.(130),

\langle\psi(q^{0};\sigma_{q}),t_{0}\rvert l^{M}(t,\mathbf{x})\lvert\psi(q^{0};\sigma_{q}),t_{0}\rangle=\frac{1-\cos 2m(t-t_{0})}{V}\langle\psi(q^{0};\sigma_{q}),t_{0}\rvert\psi(q^{0};\sigma_{q}),t_{0}\rangle\\ +\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}\frac{dq^{\prime}dq}{A}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)x^{2}}\left[-\frac{m}{E(q^{\prime})}\sin{E(q^{\prime})T}\frac{m}{E(q)}\sin{E(q)T}\right.\\ \left.+\left(\cos{E(q^{\prime})T}+i\frac{|q^{\prime}|}{E(q^{\prime})}\sin{E(q^{\prime})T}\right)\left(\cos{E(q)T}-i\frac{|q|}{E(q)}\sin{E(q)T}\right)\right].

(135)

Although, this result is similar to main text of Eq.(34), the term proportional to $1-\cos 2m(t-t_{0})$ is from the zero mode. For the last line of Eq.(134), we are able to compute the matrix elements with

	$\displaystyle\langle 0(t_{0})\|C(s_{0},t_{0})\sum_{s=\pm 1}C^{\dagger}(s,t)C(s,t)C^{\dagger}(s_{0},t_{0})\|0(t_{0})\rangle=\frac{1}{V^{2}}.$		(136)
	$\displaystyle\|c^{0}(s_{0})\|^{2}\langle 0(t_{0})\rvert C(s_{0},t_{0})l^{M}(t,\mathbf{x})C^{\dagger}(s_{0},t_{0})\lvert 0(t_{0})\rangle=\frac{\|c^{0}(s_{0})\|^{2}}{V^{2}}.$		(137)

Finally, the third line of Eq.(130) is computed to be,

\begin{split}&\langle\psi(q^{0};\sigma_{q}),t_{0}\rvert\int^{\prime}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}\sum_{s=\pm 1}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s)})(a^{\dagger}(\mathbf{k},t)e^{-i\mathbf{k}\cdot\mathbf{x}}-b(\mathbf{k},t)e^{i\mathbf{k}\cdot\mathbf{x}})C(s,t)C^{\dagger}(s_{0},t_{0})|0(t_{0})\rangle,\\ &=\int^{\prime}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}\sum_{s=\pm 1}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s)})\langle\psi(q^{0},\sigma),t_{0}\rvert a^{\dagger}(\mathbf{k},t)|0(t_{0})\rangle e^{-i\mathbf{k}\cdot\mathbf{x}}\langle 0(t_{0})|C(s,t)C^{\dagger}(s_{0},t_{0})|0(t_{0})\rangle,\\ &=\int^{\prime}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s_{0})})\langle\psi(q^{0};\sigma_{q}),t_{0}\rvert a^{\dagger}(\mathbf{k},t)|0(t_{0})\rangle e^{-i\mathbf{k}\cdot\mathbf{x}}\frac{1}{V}\cos m(t-t_{0}).\end{split}

(138)

We perform the integral over $\mathbf{k}$ using the definition of Eq.(133) for the bra- vector,

\begin{split}\int^{\prime}&\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s_{0})})\langle\psi(q^{0};\sigma_{q}),t_{0}\rvert a^{\dagger}(\mathbf{k},t)|0(t_{0})\rangle e^{-i\mathbf{k}\cdot\mathbf{x}}\\ &=\int^{\prime}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s_{0})})\int^{\prime}\frac{dqe^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i\mathbf{k}\cdot\mathbf{x}}}{\sqrt{A\sigma_{q}(2\pi)^{3/2}}\sqrt{2|q|}}\langle 0(t_{0})|a(\mathbf{q},t_{0})a^{\dagger}(\mathbf{k},t)|0(t_{0})\rangle\\ &=\int_{\mathbf{k}\in A}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s_{0})})\int^{+\infty}_{+0}\frac{dqe^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-iq\bf{e}_{2}\cdot\mathbf{x}}(\cos E(q)t+i\frac{q}{E(q)}\sin E(q)t)}{\sqrt{A\sigma_{q}(2\pi)^{3/2}}\sqrt{2|q|}}\\ &\times(2\pi)^{3}2|q|\delta(k^{1})\delta(k^{2}-q)\delta(k^{3})\\ &+\int_{\mathbf{k}\in\overline{A}}\frac{d^{3}k}{(2\pi)^{3}\sqrt{2|\mathbf{k}|}}(\phi_{-}(n_{\mathbf{k}})^{\dagger}\cdot\chi^{(s_{0})})\int^{-0}_{-\infty}\frac{dqe^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-iq\bf{e}_{2}\cdot\mathbf{x}}(\cos E(q)t+i\frac{|q|}{E(q)}\sin E(q)t)}{\sqrt{A\sigma_{q}(2\pi)^{3/2}}\sqrt{2|q|}}\\ &\times(2\pi)^{3}2|q|\delta(k^{1})\delta(k^{2}-q)\delta(k^{3})\\ &=(\phi_{-}({\bf{e}}_{2})^{\dagger}\cdot\chi^{(s_{0})})\int^{+\infty}_{+0}\frac{dq}{\sqrt{A\sigma_{q}(2\pi)^{3/2}}}e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-iq\bf{e}_{2}\cdot\mathbf{x}}(\cos E(q)t+i\frac{|q|}{E(q)}\sin E(q)t)\\ &+(\phi_{-}(-{\bf{e}}_{2})^{\dagger}\cdot\chi^{(s_{0})})\int^{-0}_{-\infty}\frac{dq}{\sqrt{A\sigma_{q}(2\pi)^{3/2}}}e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-iq\bf{e}_{2}\cdot\mathbf{x}}(\cos E(q)t+i\frac{|q|}{E(q)}\sin E(q)t),\end{split}

(139)

where $\int^{+\infty}_{+0}dq$ implies integration over positive $q$ and $\int^{-0}_{-\infty}dq$ for negative $q$ . If we choose $\chi^{(s_{0})}$ such that $\mathbf{n}=\bf{e}_{2}$ in Eq.(118),

	$\displaystyle\phi_{-}({\bf{e}}_{2})^{\dagger}\cdot\chi^{(s_{0}=1)}=0,$	$\displaystyle\phi_{-}({\bf{e}}_{2})^{\dagger}\cdot\chi^{(s_{0}=-1)}=1$		(140)
	$\displaystyle\phi_{-}(-{\bf{e}}_{2})^{\dagger}\cdot\chi^{(s_{0}=1)}=-i,$	$\displaystyle\phi_{-}(-{\bf{e}}_{2})^{\dagger}\cdot\chi^{(s_{0}=-1)}=0.$		(141)

By adding all the contribution, the expectation value for the lepton number density Eq.(130) given by the matrix element with the state Eq.(132) is

\langle\Psi(q^{0},\sigma_{q},c^{0}(s_{0}=1))|l^{M}(t,\mathbf{x})|\Psi(q^{0},\sigma_{q},c^{0}(s_{0}=1))\rangle=\frac{1-\cos 2m(t-t_{0})}{V}+\frac{|c^{0}(s_{0}=1)|^{2}}{V^{2}}\\ -i\frac{c^{0}(s_{0}=1)\cos mT}{V}\int^{-0}_{-\infty}\frac{dq}{\sqrt{A\sigma_{q}(2\pi)^{3/2}2}}(\cos(E(q)T)+i\frac{q}{E(q)}\sin(E(q)T))(e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-iqx^{2}})\\ +i\frac{{c^{0}(s_{0}=1)}^{\ast}\cos mT}{V}\int^{-0}_{-\infty}\frac{dq}{\sqrt{A\sigma_{q}(2\pi)^{3/2}2}}(\cos(E(q)T)-i\frac{q}{E(q)}\sin(E(q)T))(e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}+iqx^{2}})\\ +\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}\frac{dq^{\prime}dq}{A}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)x^{2}}\\ \times\left[\left(\cos{E(q^{\prime})T}+i\frac{|q^{\prime}|}{E(q^{\prime})}\sin{E(q^{\prime})T}\right)\left(\cos{E(q)T}-i\frac{|q|}{E(q)}\sin{E(q)T}\right)\right.\\ \left.-\frac{m}{E(q^{\prime})}\sin{E(q^{\prime})T}\frac{m}{E(q)}\sin{E(q)T}\right];

(142)

For the opposite direction of the spin state of the zero mode $s_{0}=-1$ , the result is,

\langle\Psi(q^{0},\sigma_{q},c^{0}(s_{0}=-1))|l^{M}(t,\mathbf{x})|\Psi(q^{0},\sigma_{q},c^{0}(s_{0}=-1))\rangle=\frac{1-\cos 2m(t-t_{0})}{V}+\frac{|c^{0}(s_{0}=-1)|^{2}}{V^{2}}\\ +\frac{c^{0}(s_{0}=-1)\cos mT}{V}\int_{+0}^{+\infty}\frac{dq}{\sqrt{A\sigma_{q}(2\pi)^{3/2}2}}(\cos E(q)T+i\frac{q}{E(q)}\sin E(q)T)(e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-iqx^{2}})\\ +\frac{{c^{0}(s_{0}=-1)}^{\ast}\cos mT}{V}\int_{+0}^{+\infty}\frac{dq}{\sqrt{A\sigma_{q}(2\pi)^{3/2}2}}(\cos E(q)T-i\frac{q}{E(q)}\sin E(q)T)(e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}+iqx^{2}})\\ +\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}\frac{dq^{\prime}dq}{A}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)x^{2}}\\ \times\left[\left(\cos{E(q^{\prime})T}+i\frac{|q^{\prime}|}{E(q^{\prime})}\sin{E(q^{\prime})T}\right)\left(\cos{E(q)T}-i\frac{|q|}{E(q)}\sin{E(q)T}\right)\right.\\ \left.-\frac{m}{E(q^{\prime})}\sin{E(q^{\prime})T}\frac{m}{E(q)}\sin{E(q)T}\right];

(143)

Both Eq.(142) and Eq.(143) imply that even if the zero mode amplitude $c^{0}(s_{0}=\pm 1)$ vanishes in the initial state, the expectation value has the new contribution. To examine the effect, we compute the linear density as done in the text by integrating the lepton number density over the area $A=\int\int dx^{1}dx^{3}$ . The result for Eq.(143) is,

\lambda^{M}(T=t-t_{0},x^{2})=\int dx^{1}dx^{3}\langle\Psi(q^{0},\sigma_{q},c^{0}(s_{0}=-1))|l^{M}(t,\mathbf{x})|\Psi(q^{0},\sigma_{q},c^{0}(s_{0}=-1))\rangle\\ =\frac{1-\cos 2m(t-t_{0})}{L}+\frac{|c^{0}(s_{0}=-1)|^{2}}{VL}\\ +\frac{c^{0}(s_{0}=-1)\cos mT}{L\sqrt{A}}\int_{0+}^{\infty}\frac{dq}{\sqrt{\sigma_{q}(2\pi)^{3/2}2}}(\cos E(q)T+i\frac{q}{E(q)}\sin E(q)T)(e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}-iqx^{2}})\\ +\frac{{c^{0}(s_{0}=-1)}^{\ast}\cos mT}{L\sqrt{A}}\int_{0+}^{\infty}\frac{dq}{\sqrt{\sigma_{q}(2\pi)^{3/2}2}}(\cos E(q)T-i\frac{q}{E(q)}\sin E(q)T)(e^{-\frac{(q-q^{0})^{2}}{4\sigma_{q}^{2}}+iqx^{2}})\\ +\frac{1}{\sigma_{q}(2\pi)^{3/2}}\iint^{\prime}{dq^{\prime}dq}e^{-\frac{(q^{\prime}-q^{0})^{2}+(q-q^{0})^{2}}{4\sigma_{q}^{2}}-i(q^{\prime}-q)x^{2}}\\ \times\left[\left(\cos{E(q^{\prime})T}+i\frac{|q^{\prime}|}{E(q^{\prime})}\sin{E(q^{\prime})T}\right)\left(\cos{E(q)T}-i\frac{|q|}{E(q)}\sin{E(q)T}\right)\right.\\ \left.-\frac{m}{E(q^{\prime})}\sin{E(q^{\prime})T}\frac{m}{E(q)}\sin{E(q)T}\right];

(144)

where we use $V=AL$ .

References

[1] J. N. Bahcall, AAPPS Bull. 12, no.4, 12-19 (2002) [arXiv:astro-ph/0209080 [astro-ph]].
[2] K. Abe et al. (Super-Kamiokande), Phys. Rev. D 83, 052010 (2011) [arXiv:1010.0118 [hep-ex]].
[3] B. Aharmim et al. (SNO), Phys. Rev. C 81, 055504 (2010) [arXiv:0910.2984 [nucl-ex]].
[4] B. Pontecorvo, Sov. Phys. JETP 7, 172 (1958) [Zh. Eksp. Teor. Fiz. 34, 247 (1957)].
[5] Z. Maki, M. Nakagawa and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
[6] B. Abi et al. (DUNE), Eur. Phys. J. C 80, no.10, 978 (2020) [arXiv:2006.16043 [hep-ex]].
[7] K. Abe et al. (Hyper-Kamiokande Proto-), PTEP 2015, 053C02 (2015) [arXiv:1502.05199 [hep-ex]].
[8] K. Abe et al. (Hyper-Kamiokande), PTEP 2018, no.6, 063C01 (2018) [arXiv:1611.06118 [hep-ex]].
[9] E. Baussan et al. (ESSnuSB), Nucl. Phys. B 885, 127-149 (2014) [arXiv:1309.7022 [hep-ex]].
[10] S. Ahmed et al. (ICAL), Pramana 88, no.5, 79 (2017) [arXiv:1505.07380 [physics.ins-det]].
[11] F. An et al. (JUNO), J. Phys. G 43, no.3, 030401 (2016) [arXiv:1507.05613 [physics.ins-det]].
[12] M. Aker et al. (KATRIN), Eur. Phys. J. C 80, no.3, 264 (2020) [arXiv:1909.06069 [physics.ins-det]].
[13] M. J. Dolinski, A. W. P. Poon and W. Rodejohann, Ann. Rev. Nucl. Part. Sci. 69, 219-251 (2019) [arXiv:1902.04097 [nucl-ex]].
[14] W. H. Furry, Phys. Rev. 56, 1184-1193 (1939)
[15] J. N. Bahcall and H. Primakoff, Phys. Rev. D 18, 3463-3466 (1978)
[16] J. Schechter and J. W. F. Valle, Phys. Rev. D 23, 1666 (1981)
[17] Z. z. Xing, Phys. Rev. D 87, no.5, 053019 (2013) [arXiv:1301.7654 [hep-ph]].
[18] N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
[19] J. A. Formaggio and G. P. Zeller, Rev. Mod. Phys. 84, 1307-1341 (2012) [arXiv:1305.7513 [hep-ex]].
[20] M. G. Betti et al. (PTOLEMY), JCAP 07, 047 (2019) [arXiv:1902.05508 [astro-ph.CO]].
[21] A. D. Dolgov, Phys. Rept. 370, 333-535 (2002) [arXiv:hep-ph/0202122 [hep-ph]].
[22] C. Giunti, C. W. Kim and U. W. Lee, Phys. Rev. D 44, 3635-3640 (1991)
[23] E. Akhmedov, [arXiv:1901.05232 [hep-ph]].
[24] E. Akhmedov, JHEP 07, 070 (2017) [arXiv:1703.08169 [hep-ph]].
[25] E. K. Akhmedov and A. Y. Smirnov, Phys. Atom. Nucl. 72, 1363-1381 (2009) [arXiv:0905.1903 [hep-ph]].
[26] S. Nussinov, Phys. Lett. B 63, 201-203 (1976)
[27] A. S. Adam, N. J. Benoit, Y. Kawamura, Y. Matsuo, T. Morozumi, Y. Shimizu, Y. Tokunaga and N. Toyota, PTEP 2021, (2021) [arXiv:2101.07751 [hep-ph]].
[28] A. S. Adam, N. J. Benoit, Y. Kawamura, Y. Mastuo, T. Morozumi, Y. Shimizu, Y. Tokunaga and N. Toyota, [arXiv:2105.04306 [hep-ph]].
[29] N. J. Benoit, Doctor Thesis, Hiroshima University, 2022 url:https://ir.lib.hiroshima-u.ac.jp/00053253
[30] E. J. Konopinski and H. M. Mahmoud, Phys. Rev. 92, 1045-1049 (1953)
[31] B. Pontecorvo, Zh. Eksp. Teor. Fiz. 37, 1751-1757 (1959)
[32] S. M. Bilenky and C. Giunti, Int. J. Mod. Phys. A 16, 3931-3949 (2001) [arXiv:hep-ph/0102320 [hep-ph]].
[33] B. Kayser, Phys. Rev. D 24, 110 (1981)
[34] C. Giunti and C. W. Kim, Found. Phys. Lett. 14, no.3, 213-229 (2001) [arXiv:hep-ph/0011074 [hep-ph]].
[35] C. Giunti, C. W. Kim and U. W. Lee, Phys. Rev. D 45, 2414-2420 (1992)
[36] M. Blasone, L. Smaldone and G. Vitiello, J. Phys. Conf. Ser. 1275, no.1, 012023 (2019) [arXiv:1903.01401 [hep-ph]].
[37] M. Blasone, A. Capolupo, C. R. Ji and G. Vitiello, Int. J. Mod. Phys. A 25, 4179-4194 (2010) [arXiv:hep-ph/0611106 [hep-ph]].
[38] C. Giunti, JHEP 11, 017 (2002) [arXiv:hep-ph/0205014 [hep-ph]].
[39] C. Giunti, [arXiv:hep-ph/0402217 [hep-ph]].
[40] C. Giunti, C. W. Kim, J. A. Lee and U. W. Lee, Phys. Rev. D 48, 4310-4317 (1993) [arXiv:hep-ph/9305276 [hep-ph]].
[41] M. Beuthe, Phys. Rept. 375, 105-218 (2003) [arXiv:hep-ph/0109119 [hep-ph]].
[42] A. Asahara, K. Ishikawa, T. Shimomura and T. Yabuki, Prog. Theor. Phys. 113, 385-411 (2005) [arXiv:hep-ph/0406141 [hep-ph]].
[43] E. K. Akhmedov and J. Kopp, JHEP 04, 008 (2010) [erratum: JHEP 10, 052 (2013)] [arXiv:1001.4815 [hep-ph]].
[44] J. Wu, J. A. Hutasoit, D. Boyanovsky and R. Holman, Int. J. Mod. Phys. A 26, 5261-5297 (2011) [arXiv:1002.2649 [hep-ph]].
[45] M. Blasone, A. Capolupo and G. Vitiello, Phys. Rev. D 66, 025033 (2002) [arXiv:hep-th/0204184 [hep-th]].
[46] M. Blasone, A. Capolupo, O. Romei and G. Vitiello, Phys. Rev. D 63, 125015 (2001) [arXiv:hep-ph/0102048 [hep-ph]].
[47] M. Blasone, P. Jizba and G. Vitiello, Phys. Lett. B 517, 471-475 (2001) [arXiv:hep-th/0103087 [hep-th]].
[48] K. Fujii, C. Habe and T. Yabuki, Phys. Rev. D 59, 113003 (1999) [erratum: Phys. Rev. D 60, 099903(E) (1999)] [arXiv:hep-ph/9807266 [hep-ph]].
[49] K. Fujii, C. Habe and T. Yabuki, Phys. Rev. D 64, 013011 (2001) [arXiv:hep-ph/0102001 [hep-ph]].
[50] C. Giunti, Eur. Phys. J. C 39, 377-382 (2005) [arXiv:hep-ph/0312256 [hep-ph]].
[51] A. Tureanu, [arXiv:2005.02219 [hep-ph]].
[52] M. Blasone and L. Smaldone, Mod. Phys. Lett. A 35, no.38, 2050313 (2020) [arXiv:2004.04739 [hep-ph]].
[53] T. Morozumi, N. J. Benoit and Y. Kawamura, PoS CORFU2021, 063 (2022) [arXiv:2204.00971 [hep-ph]].
[54] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou (www.nu-fit.org), JHEP 09, 178 (2020) [arXiv:2007.14792 [hep-ph]].

	$\displaystyle\frac{1}{\sqrt{2\|\mathbf{p}\|}}a_{L\alpha}(\pm\mathbf{p},t)$	$\displaystyle=\sum_{i}^{3}V_{\alpha i}\frac{\sqrt{E_{i}(\mathbf{p})+\|\mathbf{p}\|}}{2E_{i}(\mathbf{p})}\left(a_{i}(\pm\mathbf{p},-)e^{-iE_{i}(\mathbf{p})t}\pm i\frac{m_{i}}{E_{i}(\mathbf{p})+\|\mathbf{p}\|}b_{i}^{\dagger}(\mp\mathbf{p},-)e^{iE_{i}(\mathbf{p})t}\right),$		(56)
	$\displaystyle\frac{1}{\sqrt{2\|\mathbf{p}\|}}b_{L\alpha}^{\dagger}(\pm\mathbf{p},t)$	$\displaystyle=\sum_{i}^{3}V_{\alpha i}\frac{\sqrt{E_{i}(\mathbf{p})+\|\mathbf{p}\|}}{2E_{i}(\mathbf{p})}\left(b_{i}^{\dagger}(\pm\mathbf{p},+)e^{iE_{i}(\mathbf{p})t}\mp i\frac{m_{i}}{E_{i}(\mathbf{p})+\|\mathbf{p}\|}a_{i}(\mp\mathbf{p},+)e^{-iE_{i}(\mathbf{p})t}\right).$		(57)