Deriving the Pontecorvo–Maki–Nakagawa–Sakata Matrix from Koide’s Mass Formula and Brannen’s Neutrino Mass Hypothesis: Resolving the Mystery of the θ₁₃ Rotation

1. Introduction

1.1. Koide’s Mass Formula

In 1982, Yoshio Koide first proposed Koide’s mass formula [1,2] based on Harari, Haut, and Weyers [3]. This formula is expressed as

$${\displaystyle \frac{m_e+m_{\mu }+m_{\tau }}{{\left(\sqrt{m_e}+\sqrt{m_{\mu }}+\sqrt{m_{\tau }}\right)}^2}}={\displaystyle \frac{2}{3}},$$

which elegantly describes the relationship between the masses of the three charged lepton generations.

1.2. Interpretation by Carl A. Brannen

In 2006, Carl A. Brannen provided an interpretation of Koide’s mass formula [4,5]. Let the masses of $e^{-}$,${\mu }^{-}$, and ${\tau }^{-}$ be denoted as $m_{e1}$, $m_{e2}$, and $m_{e3}$, respectively, which are experimentally determined as follows [6]:

$$\begin{array}{cc}\hfill m_{e1}& =0.510999\left[\mathrm{MeV}\right],\hfill \\ \hfill m_{e2}& =105.658376\left[\mathrm{MeV}\right],\hfill \\ \hfill m_{e3}& =1776.93\left[\mathrm{MeV}\right].\hfill \end{array}$$

According to Brannen, the square root of each mass is expressed as

$$\sqrt{m_{en}}=17.716\left(1+\sqrt{2}cos\left({\displaystyle \frac{2}{9}}+{\displaystyle \frac{2}{3}}n\pi \right)\right)\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right],$$

for $n=1,2,3.$ This provides a deeper mathematical insight into Koide’s mass formula.

As $\sqrt{m_{e1}}+\sqrt{m_{e2}}+\sqrt{m_{e3}}=53.148\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right],$ which follows Brannen’s mass formula, this is used as a normalization factor. Therefore,

$$\sqrt{m_{en}}={\displaystyle \frac{53.148}{3}}\left(1+\sqrt{2}cos\left({\displaystyle \frac{2}{9}}+{\displaystyle \frac{2}{3}}n\pi \right)\right)\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right].$$

1.3. Brannen’s Neutrino Mass Hypothesis

Brannen hypothesized that a similar relationship holds for neutrinos:

$${\displaystyle \frac{m_1+m_2+m_3}{{\left(-\sqrt{m_1}+\sqrt{m_2}+\sqrt{m_3}\right)}^2}}={\displaystyle \frac{2}{3}}.$$

Let the masses of ${\nu }_1$, ${\nu }_2$, and ${\nu }_3$ be denoted as $m_1$, $m_2$, and $m_3$, respectively. Brannen proposed the following expressions [4,5]:

For $n=1$,

$$-\sqrt{m_n}={\displaystyle \frac{53.148}{3^{12}}}\left(1+\sqrt{2}cos\left({\displaystyle \frac{2}{9}}+{\displaystyle \frac{1}{12}}\pi +{\displaystyle \frac{2}{3}}n\pi \right)\right)\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right],$$

where the minus sign appears because this expression evaluates to a negative value.

For $n=2$, 3,

$$\sqrt{m_n}={\displaystyle \frac{53.148}{3^{12}}}\left(1+\sqrt{2}cos\left({\displaystyle \frac{2}{9}}+{\displaystyle \frac{1}{12}}\pi +{\displaystyle \frac{2}{3}}n\pi \right)\right)\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right].$$

From these, the neutrino masses are calculated as

$$\begin{array}{cc}\hfill m_1& =\phantom{0}0.383463\left[\mathrm{meV}\right],\hfill \\ \hfill m_2& =\phantom{0}8.913500\left[\mathrm{meV}\right],\hfill \\ \hfill m_3& =50.711844\left[\mathrm{meV}\right].\hfill \end{array}$$

1.4. Constructing Two Three-Dimensional (3D) Mass Models

Here, a question arises: While Brannen indicates that the square root of the mass of ${\nu }_1$ is negative, what does it mean for the square root of a mass to be negative?

Could it imply that ${\nu }_1$ is antimatter? No. The observed ${\nu }_1$ should correspond to the positive square root of its mass.

Assuming Brannen’s hypothesis is valid, it is proposed that $-\sqrt{m_1}$ might be the origin of the ${\theta }_{13}$ rotation in the Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix [7]. Based on Brannen’s hypothesis, two three-dimensional (3D) mass models were constructed for neutrinos.

2. Methods

In this study, most of the computations were performed using Microsoft Excel 2021 with its default double-precision arithmetic (approximately 15 significant digits). For plotting and numerical integration, Python was used. In those cases, the Excel outputs were rounded to 12 decimal places and hard-coded into the scripts. Herein, all reported values are rounded to six decimal places for readability; this rounding does not affect the accuracy of the results within reasonable numerical precision.

2.1. Construction of the Neutrino 3D Mass Models

Let the masses of ${\nu }_1$, ${\nu }_2$, and ${\nu }_3$ be denoted as $m_1$, $m_2$, and $m_3$, respectively, satisfying

$$m_1+m_2+m_3=60.008807\left[\mathrm{meV}\right]={7.746535}^2\left[\mathrm{meV}\right].$$

In 3D space, let the origin be $\mathrm{O}$; the radius is defined as

$$r=7.746535\left[{\mathrm{meV}}^{{\displaystyle \frac{1}{2}}}\right],$$

where r represents the radius of a sphere described by $x^2+y^2+z^2=r^2$.

Points ${\mathrm{P}}_{-}$ and ${\mathrm{P}}_{+}$ are defined as

$$\begin{array}{cc}\hfill {\mathrm{P}}_{-}& =\left(-\sqrt{m_1},\sqrt{m_2},\sqrt{m_3}\right),\hfill \\ \hfill {\mathrm{P}}_{+}& =\left(\sqrt{m_1},\sqrt{m_2},\sqrt{m_3}\right),\hfill \end{array}$$

where both lie on the aforementioned sphere. Further, three additional points on the sphere are defined as

$$\begin{array}{cc}\hfill {\mathrm{N}}_1& =(r,0,0),\hfill \\ \hfill {\mathrm{N}}_2& =(0,r,0),\hfill \\ \hfill {\mathrm{N}}_3& =(0,0,r).\hfill \end{array}$$

The models are constructed in two patterns, based on the square roots of the neutrino masses and the following combinations:

(i) $\left(-\sqrt{m_1},\sqrt{m_2},\sqrt{m_3}\right)$;

(ii) $\left(\sqrt{m_1},\sqrt{m_2},\sqrt{m_3}\right)$.

2.1.1. Combination $\left(-\sqrt{m_1},\sqrt{m_2},\sqrt{m_3}\right)$

Vectors and Dot Products

The unit vector (hereafter referred to as the $\sqrt{m}$ unit vector) is defined as

$$\overrightarrow{{\mathrm{OM}}_{-}}={\displaystyle \frac{\overrightarrow{{\mathrm{OP}}_{-}}}{∥\overrightarrow{{\mathrm{OP}}_{-}}∥}}.$$

The following vectors originating from $\mathrm{O}$ are defined as:

$$\begin{array}{cc}\hfill {\mathit{M}}_{-}& =\overrightarrow{{\mathrm{OM}}_{-}}=(-0.079938,0.385404,0.919279),\hfill \\ \hfill {\mathit{N}}_{1-}& =\overrightarrow{{\mathrm{ON}}_1}=(r,0,0),\hfill \\ \hfill {\mathit{N}}_{2-}& =\overrightarrow{{\mathrm{ON}}_2}=(0,r,0),\hfill \\ \hfill {\mathit{N}}_{3-}& =\overrightarrow{{\mathrm{ON}}_3}=(0,0,r).\hfill \end{array}$$

The dot products are calculated as follows:

$$\begin{array}{cc}\hfill {\mathit{N}}_{1-}·{\mathit{M}}_{-}& =-0.079938r=-0.619244\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right]=-\sqrt{m_1},\hfill \\ \hfill {\mathit{N}}_{2-}·{\mathit{M}}_{-}& =0.385404r=2.985548\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right]=\sqrt{m_2},\hfill \\ \hfill {\mathit{N}}_{3-}·{\mathit{M}}_{-}& =0.919279r=7.121225\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right]=\sqrt{m_3}.\hfill \end{array}$$

To align the direction of ${\mathit{M}}_{-}$ with the x-axis, ${\mathit{M}}_{-}$, ${\mathit{N}}_{1-}$, ${\mathit{N}}_{2-}$, and ${\mathit{N}}_{3-}$ are rotated about the origin in 3D space.

Initial Coordinates

The initial coordinates are expressed as

$$\left[\begin{array}{c}{\mathit{M}}_{-}\\ {\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ {\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}-0.079938& 0.385404& 0.919279\\ r& 0& 0\\ 0& r& 0\\ 0& 0& r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

Rotation in the $xy$-Plane

Using $cos{\alpha }_{-}=-0.203091$ and $sin{\alpha }_{-}=-0.979160$, corresponding to ${\alpha }_{-}\approx -1.775310\left[\mathrm{rad}\right]$, rotation in the $xy$-plane is applied as follows:

$$\left[\begin{array}{c}{\mathit{M}}_{-}\\ {\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ {\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}0.393607& 0& 0.919279\\ -0.203091r& -0.979160r& 0\\ \phantom{-}0.979160r& -0.203091r& 0\\ 0& 0& r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

Rotation in the $xz$-Plane

Using $cos{\beta }_{-}=0.393607$ and $sin{\beta }_{-}=-0.919279$, corresponding to ${\beta }_{-}\approx -1.166244\left[\mathrm{rad}\right]$, rotation in the $xz$-plane is applied as follows:

$$\left[\begin{array}{c}{\mathit{M}}_{-}\\ {\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ {\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ -0.079938r& -0.979160r& \phantom{-}0.186698r\\ \phantom{-}0.385404r& -0.203091r& -0.900121r\\ \phantom{-}0.919279r& 0& \phantom{-}0.393607r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

Rotation in the $yz$-Plane

Similarly, using $cos{\gamma }_{-}=0.220094$ and $sin{\gamma }_{-}=0.975479$, corresponding to $\gamma \approx -1.348885\left[\mathrm{rad}\right]$, rotation in the $yz$-plane is applied as follows:

$$\left[\begin{array}{c}{\mathit{M}}_{-}\\ {\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ {\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ -0.079938r& -0.033388r& \phantom{-}0.996241r\\ \phantom{-}0.385404r& -0.922748r& 0\\ \phantom{-}0.919279r& \phantom{-}0.383955r& \phantom{-}0.086631r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

For visualization purposes, the z-component of ${\mathit{N}}_{2-}$ can be set to 0. The x-components of ${\mathit{N}}_{1-}$, ${\mathit{N}}_{2-}$, and ${\mathit{N}}_{3-}$ denote $-\sqrt{m_1}$, $\sqrt{m_2}$, and $\sqrt{m_3}$, respectively, thus relating each vector with the respective neutrino.

2.1.2. Combination $\left(\sqrt{m_1},\sqrt{m_2},\sqrt{m_3}\right)$

Vectors and Dot Products

The unit vector is defined as

$$\overrightarrow{{\mathrm{OM}}_{+}}={\displaystyle \frac{\overrightarrow{{\mathrm{OP}}_{+}}}{∥\overrightarrow{{\mathrm{OP}}_{+}}∥}}.$$

The vectors originating from $\mathrm{O}$ are defined as follows:

$$\begin{array}{cc}\hfill {\mathit{M}}_{+}& =\overrightarrow{{\mathrm{OM}}_{+}}=(0.079938,0.385404,0.919279),\hfill \\ \hfill {\mathit{N}}_{1+}& =\overrightarrow{{\mathrm{ON}}_1}=(r,0,0),\hfill \\ \hfill {\mathit{N}}_{2+}& =\overrightarrow{{\mathrm{ON}}_2}=(0,r,0),\hfill \\ \hfill {\mathit{N}}_{3+}& =\overrightarrow{{\mathrm{ON}}_3}=(0,0,r).\hfill \end{array}$$

The dot products are calculated as follows:

$$\begin{array}{cc}\hfill {\mathit{N}}_{1+}·{\mathit{M}}_{+}& =0.079938r=0.619244\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right]=\sqrt{m_1},\hfill \\ \hfill {\mathit{N}}_{2+}·{\mathit{M}}_{+}& =0.385404r=2.985548\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right]=\sqrt{m_2},\hfill \\ \hfill {\mathit{N}}_{3+}·{\mathit{M}}_{+}& =0.919279r=7.121225\left[{\mathrm{MeV}}^{{\displaystyle \frac{1}{2}}}\right]=\sqrt{m_3}.\hfill \end{array}$$

To align the direction of ${\mathit{M}}_{+}$ with the x-axis, ${\mathit{M}}_{+}$, ${\mathit{N}}_{1+}$, ${\mathit{N}}_{2+}$, and ${\mathit{N}}_{3+}$ are rotated about the origin in 3D space.

Initial Coordinates

The initial coordinates are expressed as

$$\left[\begin{array}{c}{\mathit{M}}_{+}\\ {\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ {\mathit{N}}_{3+}\end{array}\right]=\left[\begin{array}{ccc}0.079938& 0.385404& 0.919279\\ r& 0& 0\\ 0& r& 0\\ 0& 0& r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

Rotation in the $xy$-Plane

Using $cos{\alpha }_{+}=0.203091$ and $sin{\alpha }_{+}=-0.979160$, corresponding to ${\alpha }_{+}\approx -1.366282\left[\mathrm{rad}\right]$, rotation in the $xy$-plane is applied as follows:

$$\left[\begin{array}{c}{\mathit{M}}_{+}\\ {\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ {\mathit{N}}_{3+}\end{array}\right]=\left[\begin{array}{ccc}0.393607\phantom{r}& 0& 0.919279\\ 0.203091r& -0.979160r& 0\\ 0.979160r& \phantom{-}0.203091r& 0\\ 0& 0& r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

Rotation in the $xz$-Plane

Using $cos{\beta }_{+}=0.393607$ and $sin{\beta }_{+}=-0.919279$, corresponding to ${\beta }_{+}\approx -1.166244\left[\mathrm{rad}\right]$, rotation in the $xz$-plane is applied as follows:

$$\left[\begin{array}{c}{\mathit{M}}_{+}\\ {\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ {\mathit{N}}_{3+}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0.079938r& -0.979160r& -0.186698r\\ 0.385404r& \phantom{-}0.203091r& -0.900121r\\ 0.919279r& 0& \phantom{-}0.393607r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

Rotation in the $yz$-Plane

Finally, using $cos{\gamma }_{+}=-0.220094$ and $sin{\gamma }_{+}=-0.975479$, corresponding to $\gamma \approx -1.792707\left[\mathrm{rad}\right]$, rotation in the $yz$-plane is applied as follows:

$$\left[\begin{array}{c}{\mathit{M}}_{+}\\ {\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ {\mathit{N}}_{3+}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0.079938r& \phantom{-}0.033388r& \phantom{-}0.996241r\\ 0.385404r& -0.922748r& 0\\ 0.919279r& \phantom{-}0.383955r& -0.086631r\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

The x-components of ${\mathit{N}}_{1+}$, ${\mathit{N}}_{2+}$, and ${\mathit{N}}_{3+}$ denote $\sqrt{m_1}$, $\sqrt{m_2}$, and $\sqrt{m_3}$, respectively, thus relating each vector with the respective neutrino.

3. Results

Figure 1 shows the results of the two 3D neutrino mass models.

The relationship between the two models can be expressed as follows:

$$\left[\begin{array}{c}{\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ {\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}0.984990& 0.000000& -0.172610\\ 0.000000& 1.000000& \phantom{-}0.000000\\ 0.172610& 0.000000& \phantom{-}0.984990\end{array}\right]\left[\begin{array}{c}{\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ {\mathit{N}}_{3+}\end{array}\right].$$

This simplifies to

$$\left[\begin{array}{c}{\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ {\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}cos{\theta }_{13}& 0& sin{\theta }_{13}\\ 0& 1& 0\\ -sin{\theta }_{13}& 0& cos{\theta }_{13}\end{array}\right]\left[\begin{array}{c}{\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ {\mathit{N}}_{3+}\end{array}\right].$$

Here, ${\theta }_{13}\approx -0.173479\left[\mathrm{rad}\right]\left(\approx -{9.939597}^{\circ }\right)$.

If the cyclic permutations of the 3–1, 1–2, and 2–3 rotations are considered, the relationship can also be written as

$$\left[\begin{array}{c}{\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ {\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}cos{\theta }_{31}& 0& -sin{\theta }_{31}\\ 0& 1& 0\\ sin{\theta }_{31}& 0& cos{\theta }_{31}\end{array}\right]\left[\begin{array}{c}{\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ {\mathit{N}}_{3+}\end{array}\right].$$

Here, ${\theta }_{31}\approx 0.173479\left[\mathrm{rad}\right]\left(\approx {9.939597}^{\circ }\right)$.

4. Discussion

4.1. Correspondence to the Conjugation and Parity (CP) Violation

Each component is extended to a complex value to account for the charge conjugation and parity (CP) violation [8,9]. The extended vectors are given as follows:

$$\begin{array}{cc}\hfill {\mathit{N}}_{3-}⟹-i{\mathit{N}}_{3-}& =(-0.919279ri,-0.38395ri,-0.086631ri),\hfill \\ \hfill {\mathit{N}}_{3+}⟹-i{\mathit{N}}_{3+}& =(-0.919279ri,-0.38395ri,\phantom{-}0.086631ri),\hfill \end{array}$$

which relates each vector with the respective neutrino.

Then, the relationship between the two matrices can be expressed as

$$\left[\begin{array}{c}{\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ -i{\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}cos{\theta }_{31}& 0& -isin{\theta }_{31}\\ 0& 1& 0\\ -isin{\theta }_{31}& 0& cos{\theta }_{31}\end{array}\right]\left[\begin{array}{c}{\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ -i{\mathit{N}}_{3+}\end{array}\right].$$

This relationship can also be written as

$$\left[\begin{array}{c}{\mathit{N}}_{1-}\\ {\mathit{N}}_{2-}\\ e^{-i{\displaystyle \frac{1}{2}}\pi }{\mathit{N}}_{3-}\end{array}\right]=\left[\begin{array}{ccc}cos{\theta }_{31}& 0& -e^{i{\delta }_{31}}sin{\theta }_{31}\\ 0& 1& 0\\ e^{-i{\delta }_{31}}sin{\theta }_{31}& 0& cos{\theta }_{31}\end{array}\right]\left[\begin{array}{c}{\mathit{N}}_{1+}\\ {\mathit{N}}_{2+}\\ e^{-i{\displaystyle \frac{1}{2}}\pi }{\mathit{N}}_{3+}\end{array}\right],$$

where ${\delta }_{31}={\displaystyle \frac{1}{2}}\pi \left[\mathrm{rad}\right]$.

The two states are distinguished as the states in which the square root of the mass of ${\nu }_1$ is negative, referred to as mass negative eigenstates $\left({\nu }_{1-},{\nu }_{2-},{\nu }_{3-}\right)$, and the states in which it is positive, referred to as mass eigenstates $\left({\nu }_1,{\nu }_2,{\nu }_3\right)$.

By relating each vector with the respective neutrino, the following relation can be written:

$$\left(\begin{array}{c}{\nu }_{1-}\\ {\nu }_{2-}\\ {\nu }_{3-}\end{array}\right)=\left(\begin{array}{ccc}cos{\theta }_{31}& 0& -e^{i{\delta }_{31}}sin{\theta }_{31}\\ 0& 1& 0\\ e^{-i{\delta }_{31}}sin{\theta }_{31}& 0& cos{\theta }_{31}\end{array}\right)\left(\begin{array}{c}{\nu }_1\\ {\nu }_2\\ {\nu }_3\end{array}\right).$$

Hereafter, the above matrix is denoted as

$${\mathit{U}}_{31}=\left(\begin{array}{ccc}cos{\theta }_{31}& 0& -e^{i{\delta }_{31}}sin{\theta }_{31}\\ 0& 1& 0\\ e^{-i{\delta }_{31}}sin{\theta }_{31}& 0& cos{\theta }_{31}\end{array}\right).$$

4.2. Product with the Tribimaximal Mixing Matrix

The tribimaximal mixing matrix [10] ${\mathit{U}}_{\mathrm{TBM}}$ is defined by the product of two unitary matrices, as follows:

$${\mathit{U}}_{\mathrm{TBM}}=\left(\begin{array}{ccc}{\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{3}}}\\ {\displaystyle \frac{\omega }{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{{\omega }^2}{\sqrt{3}}}\\ {\displaystyle \frac{{\omega }^2}{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{\omega }{\sqrt{3}}}\end{array}\right)\left(\begin{array}{ccc}{\displaystyle \frac{1}{\sqrt{2}}}& 0& -{\displaystyle \frac{1}{\sqrt{2}}}i\\ 0& 1& 0\\ {\displaystyle \frac{1}{\sqrt{2}}}& 0& {\displaystyle \frac{1}{\sqrt{2}}}i\end{array}\right)=\left(\begin{array}{ccc}{\displaystyle \frac{2}{\sqrt{6}}}& {\displaystyle \frac{1}{\sqrt{3}}}& 0\\ -{\displaystyle \frac{1}{\sqrt{6}}}& {\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{2}}}\\ -{\displaystyle \frac{1}{\sqrt{6}}}& {\displaystyle \frac{1}{\sqrt{3}}}& -{\displaystyle \frac{1}{\sqrt{2}}}\end{array}\right),$$

where $\omega =e^{i{\displaystyle \frac{2}{3}}\pi }$, i.e., the complex cube root of unity.

However, ${\mathit{U}}_{\mathrm{TBM}}$ does not allow extraction of rotation angles ${\theta }_{12}$ and ${\theta }_{23}$. Therefore, in this study, the following form of the tribimaximal mixing matrix is adopted:

$${\mathit{U}}_{\mathrm{TBM}}=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos{\theta }_{23}& sin{\theta }_{23}\\ 0& -sin{\theta }_{23}& cos{\theta }_{23}\end{array}\right)\left(\begin{array}{ccc}cos{\theta }_{12}& sin{\theta }_{12}& 0\\ -sin{\theta }_{12}& cos{\theta }_{12}& 0\\ 0& 0& 1\end{array}\right)=\left(\begin{array}{ccc}{\displaystyle \frac{2}{\sqrt{6}}}& {\displaystyle \frac{1}{\sqrt{3}}}& 0\\ -{\displaystyle \frac{1}{\sqrt{6}}}& {\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{2}}}\\ {\displaystyle \frac{1}{\sqrt{6}}}& -{\displaystyle \frac{1}{\sqrt{3}}}& {\displaystyle \frac{1}{\sqrt{2}}}\end{array}\right),$$

where

$$\begin{array}{cc}\hfill cos{\theta }_{12}& ={\displaystyle \frac{2}{\sqrt{6}}},sin{\theta }_{12}={\displaystyle \frac{1}{\sqrt{3}}},{\theta }_{12}\approx 0.615480\left[\mathrm{rad}\right]\left(\approx {35.264390}^{\circ }\right),\hfill \\ \hfill cos{\theta }_{23}& ={\displaystyle \frac{1}{\sqrt{2}}},sin{\theta }_{23}={\displaystyle \frac{1}{\sqrt{2}}},{\theta }_{23}={\displaystyle \frac{1}{4}}\pi \left[\mathrm{rad}\right]\left(={45.000000}^{\circ }\right).\hfill \end{array}$$

Hereafter, the above ${\mathit{U}}_{\mathrm{TBM}}$ matrix is denoted as

$${\mathit{U}}_{12}=\left(\begin{array}{ccc}cos{\theta }_{12}& sin{\theta }_{12}& 0\\ -sin{\theta }_{12}& cos{\theta }_{12}& 0\\ 0& 0& 1\end{array}\right),$$

$${\mathit{U}}_{23}=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos{\theta }_{23}& sin{\theta }_{23}\\ 0& -sin{\theta }_{23}& cos{\theta }_{23}\end{array}\right).$$

This tribimaximal mixing matrix is regarded as a transformation matrix between the mass negative eigenstates and the flavor eigenstates of neutrinos, as follows:

$$\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\\ {\nu }_{\tau }\end{array}\right)={\mathit{U}}_{23}{\mathit{U}}_{12}\left(\begin{array}{c}{\nu }_{1-}\\ {\nu }_{2-}\\ {\nu }_{3-}\end{array}\right).$$

Accordingly, the relationship between the mass and flavor eigenstates can be expressed as

$$\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\\ {\nu }_{\tau }\end{array}\right)={\mathit{U}}_{23}{\mathit{U}}_{12}{\mathit{U}}_{31}\left(\begin{array}{c}{\nu }_1\\ {\nu }_2\\ {\nu }_3\end{array}\right).$$

The product of matrices ${\mathit{U}}_{31}$, ${\mathit{U}}_{12}$, and ${\mathit{U}}_{23}$ can be approximated as follows:

$$\left(\begin{array}{ccc}\phantom{-}0.804241& 0.577350& -0.140935i\\ -0.402120-0.122054i& 0.577350& \phantom{-}0.696493+0.070468i\\ -0.402120+0.122054i& 0.577350& -0.696493+0.070468i\end{array}\right).$$

Could this be interpreted as the PMNS matrix?

The absolute values of each component are as follows:

$$\left(\begin{array}{ccc}0.804241& 0.577350& 0.140935\\ 0.420235& 0.577350& 0.700049\\ 0.420235& 0.577350& 0.700049\end{array}\right).$$

These values closely match the “Leptonic Mixing Matrix” provided by NuFIT 6.0 [11].

4.3. Neutrino Oscillation

The validity of the derived PMNS matrix depends on whether the neutrino oscillation [12] predictions calculated with it agree with the experimental data.

4.3.1. Probability Calculation

The oscillation probability formula for each neutrino can be derived as follows. Let the flavor states before and after oscillation be ${\nu }_{\alpha }$ and ${\nu }_{\beta }$, respectively. Then, the calculation involves the following steps:

(i) Decomposing flavor eigenstate ${\nu }_{\alpha }$ into mass eigenstates $\left({\nu }_1,{\nu }_2,{\nu }_3\right)$ using the PMNS matrix;

(ii) Applying phase shifts due to the time evolution of each mass eigenstate; and

(iii) Reconstructing flavor eigenstate ${\nu }_{\beta }$ from mass eigenstates $\left({\nu }_1,{\nu }_2,{\nu }_3\right)$ using the inverse PMNS matrix.

In Step (ii), the phase shift of each ${\nu }_n$ ($n=1$, 2, 3) is shifted by ${\theta }_n$, where ${\theta }_n$ depends on the neutrino mass, propagation distance, and energy. In this paper, as the estimated neutrino masses are known, the calculation proceeds in a straightforward manner.

The unit conversion constant for the phase shift is given by

$${\displaystyle \frac{1}{2\hslash c\left[\mathrm{eV}·\mathrm{m}\right]}}\approx 2.533865\times {10}^6\left[{\mathrm{eV}}^{-1}·{\mathrm{m}}^{-1}\right].$$

Thus,

$${\theta }_1\approx {\displaystyle \frac{-2.533865\times {10}^6\left[{\mathrm{eV}}^{-1}·{\mathrm{m}}^{-1}\right]·{m_1}^2·L}{E}}\left[\mathrm{rad}\right]\approx {\displaystyle \frac{-0.000373\times {10}^{-3}\left[\mathrm{GeV}·{\mathrm{km}}^{-1}\right]·L}{E}}\left[\mathrm{rad}\right],$$

$${\theta }_2\approx {\displaystyle \frac{-2.533865\times {10}^6\left[{\mathrm{eV}}^{-1}·{\mathrm{m}}^{-1}\right]·{m_2}^2·L}{E}}\left[\mathrm{rad}\right]\approx {\displaystyle \frac{-0.201317\times {10}^{-3}\left[\mathrm{GeV}·{\mathrm{km}}^{-1}\right]·L}{E}}\left[\mathrm{rad}\right],$$

$${\theta }_3\approx {\displaystyle \frac{-2.533865\times {10}^6\left[{\mathrm{eV}}^{-1}·{\mathrm{m}}^{-1}\right]·{m_3}^2·L}{E}}\left[\mathrm{rad}\right]\approx {\displaystyle \frac{-6.516319\times {10}^{-3}\left[\mathrm{GeV}·{\mathrm{km}}^{-1}\right]·L}{E}}\left[\mathrm{rad}\right],$$

where L is the propagation distance and E is the neutrino energy.

In Step (iii), the inverse of the PMNS matrix is required.

Representing the PMNS matrix as

$$\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\\ {\nu }_{\tau }\end{array}\right)=\left(\begin{array}{ccc}{U}_{e1}& U_{e2}& U_{e3}\\ U_{\mu 1}& U_{\mu 2}& U_{\mu 3}\\ U_{\tau 1}& U_{\tau 2}& U_{\tau 3}\end{array}\right)\left(\begin{array}{c}{\nu }_1\\ {\nu }_2\\ {\nu }_3\end{array}\right),$$

and noting that the PMNS matrix is unitary, its inverse is simply its Hermitian adjoint matrix. Therefore,

$$\left(\begin{array}{c}{\nu }_1\\ {\nu }_2\\ {\nu }_3\end{array}\right)=\left(\begin{array}{ccc}\overline{U_{e1}}& \overline{U_{\mu 1}}& \overline{U_{\tau 1}}\\ \overline{U_{e2}}& \overline{U_{\mu 2}}& \overline{U_{\tau 2}}\\ \overline{U_{e3}}& \overline{U_{\mu 3}}& \overline{U_{\tau 3}}\end{array}\right)\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\\ {\nu }_{\tau }\end{array}\right).$$

The oscillation probability from ${\nu }_{\alpha }$ to ${\nu }_{\beta }$, denoted as $P\left({\nu }_{\alpha }\to {\nu }_{\beta }\right)$, is given by

$$P\left({\nu }_{\alpha }\to {\nu }_{\beta }\right)={\left|\sum _{j=1}^3\left(U_{\alpha j}{e}^{i{\theta }_j}\overline{U_{\beta j}}\right)\right|}^2={\left|U_{\alpha 1}{e}^{i{\theta }_1}\overline{U_{\beta 1}}+U_{\alpha 2}{e}^{i{\theta }_2}\overline{U_{\beta 2}}+U_{\alpha 3}{e}^{i{\theta }_3}\overline{U_{\beta 3}}\right|}^2.$$

For antineutrinos, the corresponding probability is

$$P\left(\overline{{\nu }_{\alpha }}\to \overline{{\nu }_{\beta }}\right)={\left|\sum _{j=1}^3\left(\overline{U_{\alpha j}}{e}^{i{\theta }_j}{U}_{\beta j}\right)\right|}^2,$$

or alternatively:

$$P\left(\overline{{\nu }_{\alpha }}\to \overline{{\nu }_{\beta }}\right)={\left|\sum _{j=1}^3\left(U_{\alpha j}{e}^{-i{\theta }_j}\overline{U_{\beta j}}\right)\right|}^2.$$

Following the Tokai-to-Kamioka (T2K) [13,14] experimental setup, $P\left({\nu }_{\alpha }\to {\nu }_{\beta }\right)$ and $P\left(\overline{{\nu }_{\alpha }}\to \overline{{\nu }_{\beta }}\right)$ were calculated numerically using Python (see Appendix A.1 for more details). When the propagation distance is fixed at $L=295\left[\mathrm{km}\right]$, the oscillation probabilities depend on the neutrino energy (Figure 2).

4.3.2. Energy Distribution of the Muon (Anti-)Neutrino Beam

Based on this experimental setup, the emitted ${\nu }_{\mu }$ (or $\overline{{\nu }_{\mu }}$) beam energy is not precisely $0.6\left[\mathrm{GeV}\right]$; instead, it exhibits a spread in its distribution. Although the exact form of the beam energy distribution is unknown, it is assumed to be represented by a function such as

$$y=2446.048140x^2{e}^{-{\displaystyle \frac{x}{0.3}}},$$

where $x\left[\mathrm{GeV}\right]$ represents the beam energy and $y\left[{\mathrm{GeV}}^{-1}\right]$ denotes the number density of the emitted ${\nu }_{\mu }$ (or $\overline{{\nu }_{\mu }}$). This distribution is shown in Figure 3, where the function peaks at $x=0.6$. This distribution was evaluated numerically using Python (see Appendix A.2 for more details).

The expected number of emitted ${\nu }_{\mu }$ (or $\overline{{\nu }_{\mu }}$) in the energy range of 0.1–0.2 $\left[\mathrm{GeV}\right]$ is given by

$${\int }_{0.1}^{1.2}y\mathrm{d}x\approx 100.000000.$$

4.3.3. Probability Density and Expected Number of (Anti-)Neutrinos

By combining the two functions, the probability density for ${\nu }_e$ (or $\overline{{\nu }_e}$) was calculated numerically using Python (see Appendix A.3 for more details), as shown in Figure 4. The resulting distribution closely resembles that reported in Figure 3 of Ref. [14].

Integration over the 0.1–0.2 $\left[\mathrm{GeV}\right]$ interval yields the following expectations:

For 100 neutrinos, the expected number of ${\nu }_e$ events is approximately $3.921579\left[\mathrm{Events}\right]$; for 100 antineutrinos, the expected number of $\overline{{\nu }_e}$ events is approximately $2.023569\left[\mathrm{Events}\right]$.

4.4. Comparison with Experimental Data and Model Limitations

As can be seen, the expectation for ${\nu }_e$ is less than $4\left[\mathrm{Events}\right]$, which is smaller than anticipated.

Therefore, a comparison with experimental data is warranted.

Let the rotation angles in the reordered 1–2, 1–3, and 2–3 rotations be denoted as ${\theta }_{12}^{\prime }$, ${\theta }_{13}^{\prime }$, ${\theta }_{23}^{\prime }$, and ${\delta }_{\mathrm{CP}}^{\prime }$, and assume

$${\mathit{U}}_{23}^{\prime }{\mathit{U}}_{13}^{\prime }{\mathit{U}}_{12}^{\prime }={\mathit{U}}_{23}{\mathit{U}}_{12}{\mathit{U}}_{31}.$$

From the (3,1) element of the matrices, the following are obtained:

$$e^{-i{\delta }_{\mathrm{CP}}^{\prime }}sin{\theta }_{13}^{\prime }=-e^{i{\delta }_{31}}cos{\theta }_{12}sin{\theta }_{31},$$

which yields

$$sin{\theta }_{13}^{\prime }=cos{\theta }_{12}sin{\theta }_{31}\approx 0.140935,{\theta }_{13}^{\prime }\approx 0.141406\left[\mathrm{rad}\right]\left(\approx {8.101975}^{\circ }\right),$$

$${\delta }_{\mathrm{CP}}^{\prime }=-{\delta }_{31}=-{\displaystyle \frac{1}{2}}\pi \left[\mathrm{rad}\right].$$

The Daya Bay experiment, as reported by Chen [15], measured the following:

$${{sin}^2\left(2{\theta }_{13}^{\prime }\right)}_{\exp }=0.0851\pm 0.0024.$$

By contrast, the present model yields ${sin}^2\left(2{\theta }_{13}^{\prime }\right)\approx 0.077873,$ a deviation of approximately $-3\sigma$.

This discrepancy indicates that the model, in its current form, is insufficient to fully reproduce the experimental results. Furthermore, this deviation may indicate missing elements in the proposed framework, such as additional interactions, higher-order effects, or new physics beyond the present assumptions. This limitation highlights the need for theoretical refinement. Future experimental and theoretical studies will be essential for clarifying the mechanisms underlying the observed values.

5. Conclusions

In this study, the implications of Koide’s mass formula and Brannen’s neutrino mass hypothesis were examined by constructing two 3D neutrino mass models. The PMNS matrix was shown to be derived by introducing mass negative eigenstates as intermediate states linking the mass and flavor eigenstates. This approach successfully reproduces the tribimaximal mixing structure and yields a PMNS matrix with elements that closely resemble those obtained in global analyses.

Despite some encouraging results, when applied to oscillation observables, the proposed model predicts ${sin}^2\left(2{\theta }_{13}^{\prime }\right)\approx 0.077873$, which deviates from the Daya Bay measurement by approximately $-3\sigma$. This discrepancy highlights the limitations of the proposed framework and points to missing ingredients—additional interactions, higher-order effects, or new physics beyond the assumptions considered in this study.

Appendix A.1. Code for Figure 2

The following Python script was used to compute the oscillation probabilities shown in Figure 2.

Appendix A.2. Code for Figure 3

The following Python script was used to compute the beam–energy distribution shown in Figure 3.

Appendix A.3. Code for Figure 4

The following Python script was used to compute the ${\nu }_e$ (or $\overline{{\nu }_e}$) probability density shown in Figure 4 by combining the two functions.