Exact and hybrid neutrino oscillation probabilities via the evolution operator

Introduction

Neutrino oscillations are a well-established, extensively studied phenomenon that conclusively demonstrates that neutrinos are mixed, massive particles [1, 2]. The current experimental and observational data can be interpreted in terms of the minimal extension of the Standard Model, where the three known flavor eigenstates $\lvert\nu_\alpha\rangle$ ($\alpha=e,\mu,\tau$) that participate in weak interactions are linear combinations of an equal number of mass eigenstates $\lvert\nu_i\rangle$ ($i=1,2,3$), via the relation $\lvert\nu_\alpha\rangle=\sum_i U_{\alpha i}^* \lvert\nu_i\rangle$. The coefficients $U_{\alpha i}$ are the elements of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix $U=O_{23}\Gamma O_{13}\Gamma^*O_{12}$, where $O_{ij}$ are orthogonal matrices representing rotations by angles $\theta_{ij}\in[0,\pi/2]$ in the respective planes, and $\Gamma=\operatorname{diag}(1,1,e^{i\delta})$, with $\delta\in[0,2\pi)$, accounts for CP violation in the leptonic sector. Additionally, the oscillations between the active flavor neutrinos are described by two squared-mass differences: $\Delta m^2_{21}\equiv m^2_2-m^2_1$ and $\Delta m^2_{31}\equiv m^2_3-m^2_1$.

Most oscillation parameters are known to better than 5% precision thanks to global analyses of data from solar, atmospheric, reactor, and accelerator neutrino experiments carried out during the last two decades [3]. Three main unknowns remain: the sign of $\Delta m^2_{31}$, the octant of $\theta_{23}$, and the value of $\delta$. The current global fit is consistent with CP conservation within $1\sigma$ for normal ordering, whereas for inverted ordering values of $\delta$ around $3\pi/2$ are favored over CP conservation at more than $3.6\sigma$ [3]. T2K was the first experiment to exclude values of $\delta$ at $3\sigma$, ruling out a large range around $\delta=+\pi/2$ [4]. A conclusive answer is expected from a new generation of experiments, some already taking data and others under construction, that are designed to address the aforementioned questions and will also contribute to other aspects of neutrino physics and astrophysics [5, 6, 7, 8]. To this end, they will perform high-precision measurements of neutrino oscillations from different sources; JUNO, for instance, began operation in 2025 and has already delivered the most precise determination of the solar oscillation parameters [9]. Proper analysis of these measurements will, in turn, require increasingly precise determinations of event rates. While the calculations can be performed numerically with the required accuracy, having better analytical expressions remains useful for extracting as much information as possible and gaining a deeper understanding of the underlying physics.

When neutrinos propagate through a medium, coherent forward scattering with background electrons differs for $\nu_e$ and $\nu_{\mu,\tau}$, introducing a density-dependent potential energy in the Hamiltonian that governs the evolution of flavor amplitudes along the neutrino path [10, 11, 12]. Consequently, flavor transformations in matter differ non-trivially from those in vacuum, and the expressions of the corresponding probabilities are more difficult to derive [12, 13]. The problem of neutrino oscillations in matter cannot be solved exactly for an arbitrary position-dependent density. In the simplest case of two-neutrino mixing, useful expressions for the crossing probability have been derived from exact solutions for linear and exponential density profiles [14, 15]. For these profiles, exact solutions have also been obtained for three neutrinos, but they are given in terms of generalized hypergeometric functions [16, 17, 18] and are not easy to handle or interpret physically.

An important case that admits exact analytical solutions is the problem of neutrino oscillations in a medium that, to a good approximation, can be considered uniform. The complete solution for three flavors was presented in 1988 by Zaglauer and Schwarzer [19] in terms of the energy eigenvalues and mixing parameters (angles and CP phase) in matter. Subsequently, a formulation independent of the specific mixing parametrization was given by also diagonalizing the Hamiltonian in matter [20]. One drawback of such solutions is that the mixing quantities in matter have a rather complicated dependence on the vacuum parameters, making it difficult, in particular, to extract information on CP-violation effects. In 2002, Kimura, Takamura, and Yokomakura [21, 22] presented simpler exact analytical solutions for all oscillation channels, directly considering the products of the elements of the mixing matrix in matter that appear in the flavor transition amplitudes. More recently, Denton, Parke, and Zhang [23] gave very neat formulas for these products in terms of the eigenvalues of the Hamiltonian and its principal minors, with the benefit that the oscillation probabilities in matter follow from those in vacuum.

An alternative approach that does not involve the mixing parameters in matter is to avoid the diagonalization of the Hamiltonian by directly using the evolution operator in the flavor basis. In a uniform medium, this operator is given by the exponential of a constant $3 \times3$ matrix, which, thanks to the Cayley-Hamilton theorem, can be expressed as a second-order polynomial in that matrix. There exist different methods to determine the coefficients of the polynomial, such as the Lagrange interpolation formula referred to in the pioneering 1980 work by Barger et al. [12] and the method implemented by Ohlsson and Snellman [24, 25] for real $U$.

Due to the complexity of exact solutions, many studies over the years have addressed the problem by means of perturbative calculations, employing series expansions in the matter potential or in small parameters, such as $s_{13} = \sin\theta_{13}$ and/or the mass-hierarchy parameter $\Delta m^2_{21}/\Delta m^2_{31}$ [26, 27, 28, 29, 30, 31, 32, 33, 34]. In Ref. [35], the authors examined how well the various expansions work by comparing them with each other and with the exact solution. Unitarity is not usually satisfied automatically by oscillation formulas derived from standard perturbation theory; rather, it must be imposed order by order in the expansion parameter. While approximate expressions are useful for a better understanding of the problem, they are not without limitations: they have a restricted range of applicability, in the limit of zero potential they do not always reproduce the complete formulas for neutrino oscillations in vacuum, and sometimes they lack the desired simplicity.

In this work, we re-examine the problem of three-neutrino oscillations in matter with constant density in terms of the evolution operator in the flavor basis, thus avoiding the introduction of the complicated mixing angles in matter. As in Ref. [12], our starting point is the Lagrange-Sylvester formula for the exponential evolution matrix, but instead of approximate solutions, we seek the exact analytical solution. As we will show, Vieta’s formulas, expressed in terms of the traces of the Hamiltonian, are a very convenient way to simplify the calculations. Since only two of the oscillation probabilities, say $P_{\mu e}$ and $P_{\mu\tau}$, are independent [29], we derive explicit formulas for both and use them to determine the remaining probabilities via the unitarity condition. Our expressions possess an elegant, systematic, vacuum-like structure, which facilitates their application to both exact and approximate analyses. At $V = 0$, they reduce to the formulas for flavor transitions in vacuum. Furthermore, as a way of corroborating their correctness, we have also verified (not shown) that when expanding them to first order in $s_{13}$ or in the mass hierarchy, the results coincide exactly with those of Ref. [29], which were derived via series expansions of the evolution operator.

As is well known, for three neutrinos propagating in a medium, the exact energy eigenvalues are the real roots of a cubic equation, given by unenlightening trigonometric formulas. Our probabilities, however, depend on the matter effects only through these eigenvalues and the entries of the Hamiltonian, so their structure is preserved when the eigenvalues are replaced by approximate ones. This leads to a hybrid approach: exact formulas for the probabilities evaluated with approximate eigenvalues. We illustrate it with approximate eigenvalues obtained from two successive rotations of the Hamiltonian, following Ref. [34], and show that, for a DUNE-like configuration, the resulting expressions reproduce the exact probabilities with absolute errors of at most about $2 \times10^{-6}$ in all channels, for both neutrinos and antineutrinos, corresponding to relative errors below $10^{-2}\%$ in $P_{ee}$, $P_{\mu e}$, and $P_{\mu\tau}$. In these channels, the absolute errors are at least two orders of magnitude smaller than those of the zeroth-order two-rotation probabilities, which use the same eigenvalues. Furthermore, since these approximate eigenvalues reduce to the vacuum energies of the mass eigenstates for $V = 0$, our flavor transition probabilities reproduce the complete vacuum expressions in that limit, even when written in terms of the approximate eigenvalues.

The article is organized as follows. In Section 2, we introduce the evolution-operator formalism for a potential symmetric about the midpoint of the neutrino path and, as a first application, derive $P_{\mu e}$ and $P_{\mu\tau}$ in vacuum, which also allows us to introduce our notation. Having the vacuum expressions at hand, we can check the matter results by setting the charged-current potential to zero; although straightforward, this check is rarely shown explicitly in the literature. In § 3, we solve the problem of neutrino oscillations in matter of constant density, explicitly deriving each of the six independent components of the evolution operator in simplified form. These results are used in § 4 to obtain exact expressions for $P_{ee}$, $P_{\mu e}$, $P_{\mu\tau}$, and the remaining channels. In Section 5, we illustrate the hybrid approach with approximate eigenvalues and assess its performance for a DUNE-like configuration. § 6 contains a discussion of our results and the conclusions. Technical details are collected in Appendix § 8, the explicit verification of the $V \to0$ limit is given in Appendix § 9, and the derivation of the approximate eigenvalues is presented in Appendix § 10.

Flavor evolution in a symmetric potential

Let us consider neutrinos of a given flavor produced at time $t_0$. In terms of the evolution operator $\hat{U}(t,t_0)$ the state of the system $|\psi(t)\rangle$ at time $t \ge t_0$ can be expressed as $|\psi(t)\rangle= \hat{U}(t,t_0)|\psi(t_0)\rangle$, where $|\psi(t_0)\rangle= |\nu_\alpha\rangle$, with $\alpha= e,\mu,\tau$. Expanding $|\psi(t)\rangle$ in terms of the flavor basis, we have¹

with $\psi_\beta^{(\alpha)}(t) = \langle\beta|\hat{U}(t,t_0)|\alpha\rangle= U_{\beta\alpha}(t,t_0)$ and $\psi_\beta^{(\alpha)}(t_0) = \delta_{\beta\alpha}$.

For typical energies, neutrinos are ultra-relativistic particles and, following the standard approach, we introduce the position $r=t$ ($\hbar=c=1$) along the neutrino path. Then, the probability $P_{\alpha\beta}(r)=|\psi_\beta^{(\alpha)}(r)|^2$ to have a state of flavor $\beta$ at a distance $L=r-r_0$ from the production point $r_0$ can be expressed as

The quantities $U_{\beta\alpha}(r,r_0)$ are elements of the $3 \times3$ matrix $\mathcal{U}(r,r_0)$ representing the evolution operator in the flavor basis:

where for notational brevity we have omitted the arguments of the matrix elements. This matrix obeys the differential equation

where $H(r)$ is the effective Hamiltonian that governs the evolution of the flavor amplitudes. For neutrinos of energy $E$ propagating in normal matter $(e,p,n)$, after subtracting a global phase, it takes the form

where $H_0 = \operatorname{diag}(0, \Delta_{21}, \Delta_{31})$, $Y = \operatorname{diag}(1,0,0)$, and $V(r) = \sqrt{2}G_{\mathrm{F}}n_e(r)$. The first term in the previous equation determines the flavor evolution in the vacuum, while the potential $V(r)$ represents the net effect of the coherent interactions of neutrinos with the background particles. Here, $G_{\mathrm{F}}$ is the Fermi constant and $n_e(r)$ is the number density of electrons along the neutrino path. For the antineutrinos one has to replace $V(r)$ with $-V(r)$ and $U$ with $U^*$ (i.e., $\delta$ with $-\delta$) [36].

Since the matrices $\mathcal{O}_{23}$ and $\Gamma$ commute with $Y$, we can write

where $\mathcal{O}_{23}^{\mathrm{t}}$ is the transpose of $\mathcal{O}_{23}$ and