Compact Perturbative Expressions for Oscillations with Sterile Neutrinos in Matter

TL;DR

This paper introduces a simplified, accurate method for calculating neutrino oscillation probabilities in matter with sterile neutrinos, applicable to current experiments, and resolves degeneracies near resonances.

Contribution

It extends a compact oscillation calculation method to include sterile neutrinos, improving accuracy and resolving degeneracies in matter effects.

Findings

Method accurately predicts oscillations with sterile neutrinos in matter.

Zeroth order expressions are exact in vacuum and highly precise with weak matter effects.

Numerical tests confirm the method's effectiveness and potential experimental applications.

Abstract

We extend a simple and compact method for calculating the three flavor neutrino oscillation probabilities in uniform matter density to schemes with sterile neutrinos, with favorable features inherited. The only constraint of the extended method is that the scale of the matter potential is not significantly larger than the atmospheric $\Delta m^2$, which is satisfied by all the running and proposed accelerator oscillation experiments. Degeneracy of the zeroth order eigensystem around solar and atmospheric resonances are resolved. Corrections to the zeroth order results are restricted to no larger than the ratio of the solar to the atmospheric $\Delta m^2$. The zeroth order expressions are exact in vacuum because all the higher order corrections vanish when the matter potential is set zero. Also because all the corrections are continuous functions of matter potential, the zeroth order precision is much better than $\Delta m^2_\odot/\Delta m^2_\text{atm}$ for weak matter effect. Numerical tests are presented to verify the theoretical predictions of the exceptional features. Moreover, possible applications of the method in experiments to check the existence of sterile neutrinos are discussed.

Tables

Table 1. Table 1: Mixing parameters and vacuum eigenvalues used for the numerical calculations Fernandez-Martinez  et al. ( 2016 ); Kopp  et al. ( 2013 ); de Gouvêa and Kobach ( 2016 ); Aartsen  et al. ( 2016 ) . In different conventions to define the PMNS matrix (orders of U 23 subscript 𝑈 23 U_{23} and U sterile subscript 𝑈 sterile U_{\text{sterile}} , where U sterile = U 34 ​ U 24 ​ U 14 subscript 𝑈 sterile subscript 𝑈 34 subscript 𝑈 24 subscript 𝑈 14 U_{\text{sterile}}=U_{34}\,U_{24}\,U_{14} , see Eq.  5 ), some of the parameters are different, formulas to relate the parameters in both conventions are in Appendix A . In both conventions the energy eigenvalues in vacuum are Δ ​ m 21 2 = 7.5 × 10 − 5 ​ eV 2 Δ subscript superscript 𝑚 2 21 7.5 superscript 10 5 superscript eV 2 \Delta m^{2}_{21}=7.5\times 10^{-5}\,\text{eV}^{2} , Δ ​ m 31 2 = 2.5 × 10 − 3 ​ eV 2 Δ subscript superscript 𝑚 2 31 2.5 superscript 10 3 superscript eV 2 \Delta m^{2}_{31}=2.5\times 10^{-3}\,\text{eV}^{2} and Δ ​ m 41 2 = 0.1 ​ eV 2 Δ subscript superscript 𝑚 2 41 0.1 superscript eV 2 \Delta m^{2}_{41}=0.1\,\text{eV}^{2} .

$U_{\text{PMNS}}\equiv$	$s_{12}^2$	$s_{13}^2$	$s_{23}^2$	${\delta }_{23}/\pi$	$s_{14}^2$	$s_{24}^2$	${\delta }_{24}/\pi$	$s_{34}^2$	${\delta }_{34}/\pi$
$U_{\text{sterile}}{U}_{23}{U}_{13}{U}_{12}$	0.3	0.02	0.44	-0.40	0.02	0.01	0.10	0.1	0
$U_{23}{U}_{\text{sterile}}{U}_{13}{U}_{12}$			0.49	-0.39		0.02	0.50	0.09	0.08

Table 2. Table 2: Computation time consumed by different methods. Since a real time will depend on a specific computer’s performance, 0th order PZ’s (rotation method developed in this paper) time is set to be one unit time. FMN is from Fong, Minakata and Nunokawa Fong  et al. ( 2019 ) ; KP is from Klop and Palazzo Klop and Palazzo ( 2015 ) ; the analytical solution is given by Li  et al. ( 2018 ) ; the numerical method is referred to Eigen 3.3.7, http://eigen.tuxfamily.org .

Method	PZ(0)	PZ(1)	FMN(0)	KP	Analytical Solution	Numerical Method
Time Units	1.0	1.8	2.1	0.16	2.2	5.7

Equations

$$ϵ$$

$$\Delta m_{ee}^2$$

$$\textbf{H}=\frac{1}{2E}\Big{[}\textbf{U}_{\text{PMNS}}\,\text{diag}(0,\,\Delta m^{2}_{21},\,\Delta m^{2}_{31},\,\Delta m^{2}_{41},\,...,\,\Delta m^{2}_{N1})\,\textbf{U}^{\dagger}_{\text{PMNS}}+\text{diag}(a(x),\,0,\,0,\,b(x),\,...,\,b(x))\Big{]},$$

$$a$$

$$b$$

$${\text{U}}_{\text{PMNS}}\equiv {\text{U}}_{23}({\theta }_{23},{\delta }_{23}){\text{U}}_{\text{sterile}}{\text{U}}_{\text{13}}({\theta }_{13}){\text{U}}_{12}({\theta }_{12}),$$

$${\text{U}}_{\text{sterile}}^{\text{3+1}}\equiv {\text{U}}_{34}({\theta }_{34},{\delta }_{34}){\text{U}}_{24}({\theta }_{24},{\delta }_{24}){\text{U}}_{14}({\theta }_{14}).$$

$$|\tilde{\nu }⟩$$

$$=$$

$$\tilde{\text{H}}$$

$$\displaystyle=\left(\begin{array}[]{cc}\tilde{H}&\\ &\frac{M^{2}}{2E}\end{array}\right)+\tilde{\textbf{H}}_{M}.$$

$$\tilde{H}={\tilde{H}}_0+{\tilde{H}}_1.$$

$$\tilde{H}_{0}=\frac{1}{2E}\left(\begin{array}[]{ccc}\lambda_{a}&&s_{13}c_{13}\Delta m^{2}_{ee}+\epsilon\,bk_{13}c_{24}c_{34}e^{-i\delta_{34}}\\ &\lambda_{b}&\\ s_{13}c_{13}\Delta m^{2}_{ee}+\epsilon\,bk_{13}c_{24}c_{34}e^{i\delta_{34}}&&\lambda_{c}\end{array}\right),$$

$$k_{ij}\equiv \frac{s_{i4}{s}_{j4}}{ϵ}\sim \mathcal{O}(1),i,j\in \{1,2,3\}$$

$${\lambda }_a$$

$${\lambda }_b$$

$${\lambda }_c$$

$$({\tilde{H}}_1{)}_{12}$$

$$({\tilde{H}}_1{)}_{23}$$

$$({\tilde{H}}_1{)}_{13}$$

$$({\tilde{\text{H}}}_M{)}_{14}$$

$$({\tilde{\text{H}}}_M{)}_{24}$$

$$({\tilde{\text{H}}}_M{)}_{34}$$

$$({\tilde{\text{H}}}_M{)}_{44}$$

$$|\hat{\nu }⟩$$

$$={\text{U}}_{13}^{†}({\tilde{\theta }}_{13},{\alpha }_{13}){\text{U}}_{\text{sterile}}^{†}{\text{U}}_{23}^{†}({\theta }_{23},{\delta }_{23})|\nu {⟩}_{\text{f}},$$

$$\hat{\text{H}}$$

$$\textbf{U}_{13}=\left(\begin{array}[]{cc}U_{13}&\\ &1\end{array}\right).$$

$$\hat{H}\equiv U_{13}^{†}({\tilde{\theta }}_{13},{\alpha }_{13})\tilde{H}{U}_{13}({\tilde{\theta }}_{13},{\alpha }_{13}).$$

$${\hat{H}}_0$$

$$\displaystyle=\frac{1}{2E}\left(\begin{array}[]{ccc}\lambda_{-}&&\\ &\lambda_{0}&\\ &&\lambda_{+}\end{array}\right),$$

$${\hat{H}}_1\equiv U_{13}^{†}({\tilde{\theta }}_{13},{\alpha }_{13}){\tilde{H}}_1{U}_{13}({\tilde{\theta }}_{13},{\alpha }_{13}).$$

$${\lambda }_{\mp }$$

$${\lambda }_0$$

$$\cos 2{\tilde{\theta }}_{13}$$

$${\alpha }_{13}$$

$$({\hat{H}}_1{)}_{12}$$

$$\displaystyle\times\bigg{\{}$$

Full text

††preprint: FERMILAB-PUB-19-042-T

I Introduction

Since the discovery of neutrino oscillations, Fukuda et al. (1998), which determined that neutrinos are massive particles, many studies of neutrino scenarios beyond the three-flavor Standard Model have been performed.

One promising solution to the origin of the neutrino masses is a theoretical scheme with additional sterile neutrinos. In such a scheme neutrino oscillations will be modified because of the additional mixing with sterile neutrinos. In matter, calculations of neutrino propagation will be significantly more complicated since the sterile neutrinos also change the Wolfenstein matter effect term Wolfenstein (1978) in the Hamiltonian. There have been some analytical derivations of the matter effect in a 3+1 scenario, i.e. one sterile neutrino Li et al. (2018) in addition to the 3 active ones. However, the exact analytical solutions are impossible for more than one sterile neutrino because a quintic or even higher order equation will be encountered. Consequently, alternative perturbation approaches should be considered.

A satisfying perturbative framework, regardless of the existence of sterile neutrinos, is expected to possess the following properties: the expansion parameter is small; crossings of zeroth order eigenvalues are avoided anywhere; the approximated values go to the exact ones in vacuum. Recently, a compact perturbative framework achieving all the objectives above was developed by Denton, Minakata and Parke (DMP) to calculate the propagation of neutrinos in matter under the assumption of the standard three-flavor scheme Minakata and Parke (2016); Denton et al. (2016, 2018).

The main focus in this paper is to extend the principle and method of the DMP framework to schemes with sterile neutrinos when the scale of matter potential $a$ is smaller or comparable to $\Delta m^{2}_{\text{atm}}$, which is the case of all running and proposed accelerator neutrino oscillation experiments. The expansion parameter Denton et al. (2016); Nunokawa et al. (2005); Parke (2016), which will be retained by the extension, is

[TABLE]

The perturbative Hamiltonian will have no diagonal elements, and all its off-diagonal elements are proportional to $\epsilon$ and vanish in vacuum. Crossings of the zeroth order active eigenvalues will be resolved by a series of real or complex rotations; whereas crossings of the large sterile eigenvalues will not be considered since this will happen only if the matter effect is extremely large.

The structure of this paper is listed as follows: in Section II, we derive details of the rotations. This gives the zeroth order PMNS matrix and eigenvalues. The perturbative Hamiltonian is also determined by the rotations. In Section III, we discuss the higher order corrections by perturbative expansions after the rotations. A numerical test will also be presented to verify the predicted precision. In Section IV we use these perturbative expressions to calculate the oscillation probabilities of different channels and baselines. Moreover, potential applications of the method are discussed. We compare our method to some former works in Section V. Section VI is the conclusion. All other remarks and supplementary materials that are useful can be found in the Appendices.

II Rotations to derive zeroth order approximations and perturbative Hamiltonian

The principle of the method in Minakata and Parke (2016); Denton et al. (2016, 2018) is that by implementing a series of rotations of the Hamiltonian, one can disentangle the crossings of the diagonal elements and diminish the off-diagonal elements to arbitrary scales. In particular:

In the given Hamiltonian in flavor basis, find the sector with leading order (largest absolute value) off-diagonal element, then perform a rotation to diagonalize this sector. 2. 2.

Use the rotated Hamiltonian to replace the initial one and repeat the process until all the off-diagonal elements are smaller than the expected scale and the diagonal element crossings are eliminated.

In principle, the above process is not designated to any specific dynamical system and is also applicable to the schemes with sterile neutrinos.

However, this scheme must be implemented with considerable care otherwise the resulting analytic expressions becoming extremely long and complicated. First, one has to carefully choose the extension to the PMNS matrix to include sterile neutrinos as the standard choice here is far from optimal. Second, one has to decide whether or not one deals with all level crossings of the diagonal elements of the Hamiltonian or restrict the range of applicability of the result. We address these issues in depth in the following sub-sections.

II.1 PMNS matrix in vacuum

If we assume a $3+N$ scheme, i.e. there are $N$ sterile neutrinos in the scheme, the Hamiltonian in the flavor basis will be

[TABLE]

where $a$ and $b$ are Wolfenstein’s matter potentials Wolfenstein (1978):

[TABLE]

For earth matter the neutron number density $N_{n}$ is approximately equal to the electron number density $N_{e}$ so that $b\approx a/2$.

The PMNS matrix $\textbf{U}_{\text{PMNS}}$ in vacuum, which relates the flavor basis and the mass basis, is the product of a series of (complex) rotations Maki et al. (1962); Pontecorvo (1968). In the Standard Model, the convention is chosen to be $\textbf{U}^{\text{SM}}_{\text{PMNS}}\equiv\textbf{U}_{23}\,\textbf{U}_{\text{13}}\,\textbf{U}_{12}$. In the $3+N$ scheme there will be extra rotations mixing with sterile neutrinos. It is natural to require that the convention is equivalent to that of the $3\nu$SM case if all the extra rotations are trivial. Therefore, we will keep the relative positions of the three rotation matrices in the active sector when defining the PMNS matrix with sterile neutrinos.

Also it is observed that both the second and the third row vanish in the matter potential term in Eq. 2, thus we will keep $\textbf{U}_{23}$ as the first rotation in the PMNS matrix so the R.H.S of Eq. 2 will be independent of the 2-3 mixing parameters if we perform the $\textbf{U}_{23}$ rotation. The last step to determine the convention of the PMNS matrix is finding places after the $\textbf{U}_{23}$ for the rotations mixing with the sterile neutrinos. By trying different choices to simplify the calculation processes, we adopt the following convention of the PMNS matrix:

[TABLE]

where $\textbf{U}_{\text{sterile}}$ is the product of all the rotations mixing with sterile neutrinos. This choice leads to significant reductions in the complexity of the calculations and the resulting expressions. Physics, of course, is independent of this choice.

In the following sections, we will use the 3+1 scheme as an example to develop the expressions for the schemes with sterile neutrinos. In particular, we choose 111Convention of the CP phases is chosen to simplify the calculation process. Different conventions can be related by pure phase transformations.

[TABLE]

Current global fits Dentler et al. (2018); Dev et al. (2019); Diaz (2017) suggest $|U_{i4}|\sim 0.1$, so in this paper we assume that $\textbf{U}_{\text{sterile}}\simeq\mathbb{1}+\mathcal{O}(\sqrt{\epsilon})$, which means that $s_{i4}\sim\mathcal{O}(\sqrt{\epsilon})$ for $i=1,\,2,\,3$. The small parameter $\epsilon$ is defined in Eq. 1.

The convention in Eq. 4 is different from the usual one used by many papers in which $\textbf{U}_{\text{sterile}}$ comes before (i.e. on the left side of) all the three rotations in the active sector (see e.g., Blennow et al. (2018)). We will derive the relations of the mixing angles and phases connecting both conventions in Appendix A.

II.2 $\textbf{U}_{23}$ and $\textbf{U}_{\text{sterile}}$ rotations

We first define a rotated basis $|\tilde{\nu}\rangle$ by

[TABLE]

$|\nu\rangle_{\text{f}}$ is the flavor basis. After the rotations, the Hamiltonian becomes

[TABLE]

In the above equation $M^{2}(b)\equiv\Delta m^{2}_{41}+b\,c^{2}_{14}c^{2}_{24}c^{2}_{34}$, $\tilde{H}$ is a $3\times 3$ submatrix in the active sector and in $\tilde{\textbf{H}}_{M}$ all the elements not in the 4th column or row vanish.

Based on the scales we can distribute the elements of $\tilde{H}$ into two parts, i.e.

[TABLE]

The leading order term is

[TABLE]

where

[TABLE]

and the diagonal elements, which can be approximations to the eigenvalues are

[TABLE]

In the first order term $\tilde{H}_{1}$, all the diagonal elements vanish, and the off-diagonal elements are

[TABLE]

None-zero elements of $\tilde{\textbf{H}}_{M}$ are listed below (the Hamiltonian is a Hermitian matrix)

[TABLE]

Since $s_{i4}\sim\mathcal{O}(\sqrt{\epsilon})$, it is easy to see that $\tilde{\textbf{H}}_{M}\sim\mathcal{O}(\sqrt{\epsilon})$. Although $\tilde{\textbf{H}}_{M}$ is not as small as $\mathcal{O}(\epsilon)$, it will be a part of the perturbative Hamiltonian. However, this does not mean that the first order corrections must be as large as $\mathcal{O}(\sqrt{\epsilon})$. The mass of the heavy sterile neutrino will be an alternative parameter which controls scales of the correction terms. More specifically, in a perturbative expression, all non-zero elements of $\tilde{\textbf{H}}_{M}$ are divided by $M^{2}$. For large $M^{2}$ the quotient gives a small term in the perturbation expansion. Another condition that is necessary for $\tilde{\textbf{H}}_{M}$ being a perturbative Hamiltonian is that it consists of terms proportional to $a$ and $b$, which means that it vanishes in vacuum. This is crucial because we require the perturbative expressions to be exact in vacuum.

II.3 $\textbf{U}_{13}$ rotation

Now the dominating off-diagonal term (except the ones in $\tilde{\textbf{H}}_{M}$) comes from the (1-3) sector of $\tilde{H}_{0}$. Because of the complex phase $\delta_{34}$, the rotation will not be real. Let us assume that the rotation is $\textbf{U}_{13}(\tilde{\theta}_{13},\alpha_{13})$, where $\tilde{\theta}_{13}$ is the real rotation angle and $\alpha_{13}$222Here we are not using the usual phase symbol $\delta$ since $\alpha_{13}$ is not an effective physical phase in matter. In Appendix B it can be eliminated by implementing a pure phase transformation of the neutrino basis. is the complex phase. After this rotation the neutrino basis becomes

[TABLE]

where $\textbf{U}^{\dagger}_{\text{sterile}}=\textbf{U}^{\dagger}_{14}(\theta_{14})\textbf{U}^{\dagger}_{24}(\theta_{24},\delta_{24})\textbf{U}^{\dagger}_{34}(\theta_{34},\delta_{34})$. The Hamiltonian becomes

[TABLE]

Since the 4th index is not engaged in the rotation, we can just focus on the first three indices and define a $3\times 3$ submatrix $U_{13}$ to be the active sectors of $\textbf{U}_{13}$, i.e.

[TABLE]

After the rotation, the sub-Hamiltonian in the active sector $\tilde{H}$ becomes

[TABLE]

We require the (1-3) sector of $\tilde{H}$ to be diagonalized by $U_{13}(\tilde{\theta}_{13},\alpha_{13})$. Since the (1-3) sector of $\tilde{H}_{1}$ vanishes, it is equivalent to diagonalizing this sector of $\tilde{H}_{0}$, i.e.

[TABLE]

with $\lambda_{\pm}$ and $\lambda_{0}$ to be determined. Simultaneously $\tilde{H}_{1}$ becomes

[TABLE]

It can be shown that

[TABLE]

The real rotation angle and the complex phase can be determined by

[TABLE]

The elements of $\hat{H}_{1}$ are

[TABLE]

The Hamiltonian in the sterile sector becomes

[TABLE]

At the end of this subsection we define a real parameter $\epsilon^{\prime}$ and a phase $\alpha_{\epsilon}$

[TABLE]

Obviously $\epsilon^{\prime}\sim\epsilon$ and $(\hat{H}_{1})_{23}=e^{i\alpha_{\epsilon}}\epsilon^{\prime}\Delta m^{2}_{ee}/2E$. It is not hard to see that in the Standard Model $\epsilon^{\prime}=|\epsilon\,\sin{(\tilde{\theta}_{13}-\theta_{13})}\,s_{12}c_{12}|$, which reconciles with the one defined in Denton et al. (2016). The two new defined parameters will frequently emerge in the following sections. Since in vacuum, $a,\,b=0$, $\tilde{\theta}_{13}=\theta_{13}$ and $\alpha_{13}=0$, $\epsilon^{\prime}$ must be zero then, as shown in Fig. 1. This guarantees that the perturbative expressions will be exact in vacuum.

II.4 $\textbf{U}_{12}$ rotation

As pointed out in Denton et al. (2016), to resolve the $\lambda_{1}$ and $\lambda_{0}$ crossing at the solar resonance, one more rotation that diagonalizes the (1-2) sector is necessary. Again, since $(\hat{H}_{1})_{12}$ is complex, the rotation cannot be real in general. We assume that the rotation in (1-2) sector is $\textbf{U}_{12}(\tilde{\theta}_{12},\alpha_{12})$, and after this rotation, the neutrino basis becomes

[TABLE]

where $\textbf{U}_{\text{sterile}}=\textbf{U}^{\dagger}_{14}(\theta_{14})\textbf{U}^{\dagger}_{24}(\theta_{24},\delta_{24})\textbf{U}^{\dagger}_{34}(\theta_{34},\delta_{34})$ The Hamiltonian becomes

[TABLE]

Similar to the case of the (1-3) rotation, we can again define a $3\times 3$ submatrix $U_{12}$ by

[TABLE]

Now we require the $U_{12}(\tilde{\theta}_{12},\alpha_{12})$ to diagonalize the (1-2) sector of $\hat{H}$. After the rotation the sub-Hamiltonian is

[TABLE]

where $\check{H}_{0}$ and $\check{H}_{1}$ are in zeroth and first order respectively, i.e.

[TABLE]

The diagonal elements of $\check{H}_{0}$ are

[TABLE]

The real rotation angle and the complex phase can be determined by

[TABLE]

Values of $\sin^{2}\tilde{\theta}_{13}$ and $\sin^{2}\tilde{\theta}_{12}$ are plotted in Fig. 2. After this (1-2) rotation, crossings of the first two diagonal elements $\lambda_{1,2}$ have been resolved, as shown in the top panels of Fig. 3. They will be the zeroth order eigenvalues in the following perturbation expansions in the next section. The difference between 3+1 and $3\nu$SM is small in both panels of Fig. 2 and the bottom panels of Fig. 3 but not insignificant.

The Hamiltonian in the sterile sector now is

[TABLE]

From $\tilde{\textbf{H}}_{M}$ to $\check{\textbf{H}}_{M}$, we implemented two rotations in (1-3) and (1-2) sectors. Because the active and sterile sectors were not mixed by the two rotations, the elements are still combinations of the terms proportional to $s_{i4}\sim\mathcal{O}(\sqrt{\epsilon})$. Elements of $\check{\textbf{H}}_{M}$ can be found in Appendix C.

II.5 Crossings of $M^{2}$

In principle, there are still some possible crossings of the diagonal elements, namely the crossings to the fourth diagonal element. Since both the (1-3) and the (1-2) rotations are in the active space (first three rows and columns), the fourth element is still

[TABLE]

since $\Delta m^{2}_{41}$ is much larger than the active eigenvalues in vacuum. Thus, the crossings to $M^{2}$ can only happen with very high neutrino energy, as shown in the top panels of Fig. 3. From the figure we can see that if $Y_{e}\rho=1.4\,\text{g}\,\cdot\text{cm}^{-3}$, for the earth’s crust, the neutrino energy must be ${\cal O}(1)\,\text{TeV}$. Considering the energy scales of the current and future accelerator based oscillation experiments, we are therefore not considering the energy region of these additional crossings, so they will not effect our result. For much higher energy experiments these additional level crossings would have to be dealt with using matter additional rotations.

II.6 Summary of the rotations

Now $\check{H}_{0}$’s diagonal elements, $\lambda_{1,2,3}$, do not cross (crossings to $M^{2}$ will not happen in the energy region of interest). All the off-diagonal elements in the active sectors are of scale $\epsilon^{\prime}$. We will distribute all the diagonal elements to the zeroth order Hamiltonian and all the off-diagonal elements to the perturbative Hamiltonian, i.e.

[TABLE]

The zeroth order effective PMNS matrix in matter is

[TABLE]

Since all possible degeneracies have been removed in the energy scale which we are interested in, we are free to implement a perturbation expansion to achieve even better accuracy. The process of reducing errors by performing rotations and perturbative expansions is summarized in Fig. 4.

For the scenario with more than one sterile neutrino, although it is more complicated, the rotation method developed here is still applicable. Using the same convention of Eq. 4 to define the PMNS matrix and implement the rotations in the sequence of $U_{23}\rightarrow U_{\text{sterile}}\rightarrow U_{13}\rightarrow U_{12}$ as above.

III Perturbative expressions

Since all the crossings of the zeroth order eigenvalues have been resolved (except for the crossings with $M^{2}$, which are not in the energy region of interest) by the rotations and all the off-diagonal elements are small, we can now calculate the higher order corrections to the eigenvalues and eigenvectors by perturbation methods .

We define V to be the exact PMNS matrix in matter. It can be related to the zeroth order $\textbf{U}^{m}_{\text{PMNS}}$ by

[TABLE]

where $\textbf{W}_{n}$ is $n$th order correction. The exact eigenvalues are

[TABLE]

where $\lambda_{1,2,3}$ are defined in Eq. 42 and $\lambda_{4}=M^{2}$, $\lambda^{(n)}_{i}$ is the $n$th order correction.

First order corrections to the eigenvalues are

[TABLE]

First order corrections to the eigenstates are determined by $\textbf{W}_{i}$ defined in Eq. 48, which are

[TABLE]

The detailed first and second order formulas of the perturbation expansions can be found in Appendix D. In general, with crossings of the zeroth order eigenvalues ruled out, perturbative expansions can go to arbitrary precision. However, numerical tests will suggest that it is sufficient to terminate the approach at second order.

III.1 Numerical precision test

We now test the accuracy of our perturbative expressions. We choose the $\nu_{\mu}\rightarrow\nu_{e}$ channel and 1300km baseline of DUNE to do the numerical test. The density of the earth crust is chosen to be $Y_{e}\rho=1.4\,\text{g}\cdot\text{cm}^{-3}$, $b=a/2$ and all the mixing parameters are listed in Table 1. The exact oscillation probabilities can be figured out by Li et al. (2018) or given by a computer algebra system333Only considering the 3+1 scheme, an analytical solution is still possible since one just need solve a quartic equation; but it is not the case for schemes with more sterile neutrinos. The results are shown in Fig. 5. The error in the zeroth order expression is expected to be no more than $\epsilon\sim 10^{-2}$, which is confirmed by the red curve in the plot; the green curve depicts the error of the first order perturbative expansion, which is under $\epsilon^{2}\sim 10^{-4}$; to second order, the error further declines to $\epsilon^{3}\sim 10^{-6}$, which also coincides with the prediction. In Fig. 5 the expectation values are obtained by averaging over the fast oscillation terms, i.e. the terms with angular velocities proportional to $(\lambda_{4}-\lambda_{i})$. More specifically,

[TABLE]

Based on the numerical results, we confirm that at least the second order perturabtive expansion is significantly more accurate than any experimental results Acciarri et al. (2015); Patterson (2012); Abe et al. (2011, 2015); Kelly and Parke (2018).

IV Oscillation probabilities and detecting sterile neutrinos

In this section we will discuss a possible application of the perturbative expressions above for detecting sterile neutrinos. The principle of the approach is that one can calculate the theoretical predictions of the oscillation probabilities in different schemes and compare them with the experimental results. Usually for a given baseline and neutrino energy the predictions from different schemes are close, therefore it is essential to figure out sufficiently accurate expressions for the oscillation probabilities. A similar discussion can be found in Fong et al. (2019).

In a scheme with $N$ sterile neutrinos, the neutrino oscillation probabilities for $\nu_{\alpha}\rightarrow\nu_{\beta}$ ($\alpha,\beta\in\{e,\mu,\tau\}$) are

[TABLE]

where $\lambda^{(ex)}_{i}$ are exact eigenvalues. We can chose the zeroth order results as an approximation, i.e. we adopt

[TABLE]

where $U^{m}_{\text{PMNS}}$ is defined in Eq. 47 and

[TABLE]

where $\lambda_{1,2,3}$ are defined in Eq. 42 and $\lambda_{4}=M^{2}(b)$. For the mass of the sterile neutrino, since it is significantly larger than the active ones, the oscillations related to it will be averaged out.

Former and running experimental facilities have provided parameter fitting results of neutrino oscillations for different schemes. With these parameters, for future baselines, one can predict the probabilities in different schemes and this is a potential approach to determine the existence of sterile neutrinos Fong et al. (2019). We present the probabilities given by the 3+1 scheme and the differences of the probabilities |\big{<}P_{3+1}\big{>}-P_{3\nu\text{SM}}|, in different channels, in Fig. 6, 7 and 8. The probabilities in the Standard Model are given by Zaglauer and Schwarzer (1988); Barger et al. (1980); the 3+1 scheme is calculated by the 0th order rotation method developed in this paper. All the parameters are given in Table 1.

In the figures we can identify several regions in which the differences are significantly larger than errors of the perturbation expansions. For example, in the $\nu_{\mu}\rightarrow\nu_{e}$ channel, around the band of $L/E\simeq\,700\,(\text{km/GeV})$, |\big{<}P_{\text{3+1}}\big{>}-P_{3\nu\text{SM}}| may be larger than 0.02, the differences will be even larger than 0.05 if $L/E\gtrsim 1500\,\text{km/GeV}$ and the baseline is longer than 500 km. In this channel baselines of T2K/HyperK, NOVA and DUNE (estimated) are marked Abe et al. (2017); Acero et al. (2019); Papadimitriou et al. (2018). For the channel of $\nu_{\mu}\rightarrow\nu_{\mu}$, shifts from the $3\nu$SM will be more than 0.05 with $L/E\simeq\,1000(\text{km/GeV})$ and the baseline is longer than $1000\,\text{km}$. For the $\nu_{\mu}\rightarrow\nu_{\tau}$ channel, the scale of the greatest difference is larger than 0.16 if $L/E\simeq\,500(\text{km/GeV})$ or $\simeq\,1500(\text{km/GeV})$. Future experiments may measure the oscillation probabilities with baselines and neutrino energies in the region of interest predicted above and compare the results with the numerical outcomes.

V Compare to existing approximation formulas

Approximation methods to calculate neutrino oscillations in matter in 3+1 scheme have been studied by many researchers, for example see Klop and Palazzo (2015); Agarwalla et al. (2016a, b); Fong et al. (2017, 2019). All these works chose to ignore the Hamiltonian’s fourth row/column (except for the fourth diagonal element) in zeroth order approximations; thus the problem was reduced to $3\times 3$ case. However, solving a three dimensional eigensystem is still not simple (many approximation methods developed for the $3\nu SM$ scheme). Fong, Minakata and Numokawa (FMN) adopted the exact 3-dim solutions Fong et al. (2019) which was complicated, see Zaglauer and Schwarzer (1988); Kimura et al. (2002). Klop and Palazzo (KP) used one more approximation method for the $3\times 3$ submatrix Klop and Palazzo (2015) which introduced extra errors. Fig. 9 compares fractional errors of the rotations method developed in this paper (PZ) with 0th order FMN and KP, assuming baselines of T2K/Hyper-K (T2K/HK), NOVA, HyperK-Korea (T2HKK) and DUNE. Compared with KP, just to 0th order PZ is significantly more precise for almost all baselines and energy ranges. Based on Fig. 9, FMN’s precision is similar to PZ’s for most baselines and energy regions, however, we can still identify PZ’s advantage for T2K’s baseline or low energy ($\lesssim 1\text{GeV}$) regions.

We compare computing time of the different methods in Table 2, since any specific computing time heavily depends on performance of the computer, we provide a list of relative computing time, i.e 0th order PZ’s computing time is set to one unit time. The speed of a numerical method (using Eigen 3.3.7, http://eigen.tuxfamily.org) is also included in the comparison. A similar comparison of the methods for the $3\nu$SM scheme can be found in Barenboim et al. (2019). Table 2 shows that compared with the rotation method developed in this paper (PZ), only the KP method is faster, however, its advantage in speed will be offset by the relatively poor precision. The FMN method is comparable in time consumed to the exact analytical solution. For experimentalists, the speed of a numerical method for evaluating the oscillation probability is relevant because it effects the time and computing resources consumed by large multi-dimensional parameter scans.

Besides simplicity and better precision, the rotation method of this paper also gives explicit expressions of zeroth order eigenvalues and mixing angles and CP phases with high precision which are not covered by any former references.

VI Conclusions

A compact and simple technique for calculating neutrino oscillation probabilities in matter for schemes with sterile neutrino has been developed from the extension of an analogous method for the $3\nu$SM model Denton et al. (2016). The extended method is appropriate to conditions in which the Wolfenstein matter potentials defined in Eq. 3 are not significantly larger than $\Delta m^{2}_{\text{atm}}$, meaning that it may be applied to all the current and proposed accelerator neutrino oscillation experiments. The zeroth order eigensystem of the Hamiltonian in the active space (i.e. the three dimensions included in $3\nu$SM) derived by the method is non-degenerate. Meanwhile, numerical study shows that crossings of the zeroth order eigenvalues involving the sterile one only happen with large matter potential (high neutrino energy for the earth’s crust), which is out of the paper’s scope of discussion. An additional crucial advantage of the method developed in this paper is that errors of the zeroth order results are small when the matter potentials are $\leq\Delta m^{2}_{ee}$ and vanish in vacuum because the matter potential terms are factors of all the perturbative terms.

We implement a series of complex or real rotations to kill the leading order off-diagonal elements and resolve crossings of the diagonal elements of the Hamiltonian. The rotation angles and phases are the zeroth order mixing parameters of the effective PMNS matrix in matter. In the rotated Hamiltonian, the diagonal elements are the zeroth order eigenvalues, whereas the off-diagonal elements are the perturbing Hamiltonian. Based on this arrangement, perturbation expansions are performed after the rotations to achieve better accuracy. When the matter effect is comparable to the vacuum mixing effect, i.e. matter potentials defined in Eq. 3 are comparable to $\Delta m^{2}_{31}$, the expansion parameter is no larger than $\Delta m^{2}_{21}/\Delta m^{2}_{31}\simeq 0.03$; when the matter effect itself is weak, the perturbative Hamiltonian will be higher order because it consists of terms proportional to $a\epsilon$ or $b\epsilon$.

Finally numerical tests show that the rotation method developed in this paper balances precision with computing speed nicely; comparisons with a numerical method and some previous approximate methods have been applied to demonstrate the rotation method’s advantages as presented in this paper. For the first order perturbation expansion of this paper, absolute errors of the oscillation probabilities are shown to be no more than $10^{-4}$. This precision is sufficient to distinguish the schemes with sterile neutrinos from the $3\nu$SM model, which makes the method developed by this paper suitable to explore the existence of sterile neutrinos.

Acknowledgements

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.

This project has received funding/support from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 690575 and No 674896.

We thank Peter Denton for discussions and comments on this paper and Heath O’Connell for a careful reading of this manuscript. X.Z. thanks Melissa Clegg for discussing the English expressions in this paper.

Appendix A Mixing angles and phases in the new convention of the PMNS matrix

The PMNS matrix in the new and the usual convention can be expressed as

[TABLE]

We will express the parameters of the new convention (without the prime) in formulas of the parameters of the usual convention (with prime). We notice that $\textbf{U}_{14}\textbf{U}_{23}=\textbf{U}_{23}\textbf{U}_{14}$ then

[TABLE]

and

[TABLE]

where A is a traceless real diagonal matrix. Solving Eq. 58 we get the following relations

[TABLE]

and

[TABLE]

The approximated formulas, with $\mathcal{O}(\epsilon)$ corrections, are also listed below

[TABLE]

Appendix B Complex phases convention

In Section II.1 we chose $\textbf{U}_{12}$ and $\textbf{U}_{13}$ to be real; however, now $\alpha_{12}$ and $\alpha_{13}$ are non-zero. To recover the initial convention of the complex phases we need to implement a phase transformation. Firstly, we multiply the 1st row by $e^{-i\,\alpha_{12}}$ and the 1st column by $e^{i\,\alpha_{12}}$; then the 3rd row is multiplied by $e^{i(\alpha_{13}-\alpha_{12})}$ and the 3rd column is multiplied by $e^{i(-\alpha_{13}+\alpha_{12})}$. Finally all the complex phases are absorbed into $\textbf{U}_{23}$, $\textbf{U}_{24}$ and $\textbf{U}_{34}$. The zeroth order phases are

[TABLE]

Appendix C Elements of $\check{\textbf{H}}_{M}$

Since the Hamiltonian must be Hermitian, we will just present the 4th column.

[TABLE]

Appendix D Perturbation expansions

D.1 First order corrections

Since all diagonal elements have been absorbed into the zeroth order Hamiltonian, by Eq. 50, the first order corrections to the eigenvalues, which are the diagonal elements of the perturbative Hamiltonian, are zero, i.e.

[TABLE]

As for the eigenvectors, to first order, are derived from Eq. 51 as follows:

[TABLE]

D.2 Second order corrections

The second order corrections to eigenvalues are

[TABLE]

or explicitly as

[TABLE]

where $\Delta\lambda_{ij}\equiv\lambda_{i}-\lambda_{j}$.

The second corrections to eigenvectors can be calculated by the corrections to the PMNS matrix:

[TABLE]

We list the elements of $\textbf{W}_{2}$ below:

[TABLE]