Neutrino oscillation probabilities through the looking glass

TL;DR

This paper reviews various methods for calculating neutrino oscillation probabilities in matter, assessing their accuracy and efficiency to guide future experimental and phenomenological analyses.

Contribution

It provides a comparative analysis of exact and approximate oscillation probability expressions, identifying those suitable for next-generation experiments.

Findings

Many existing expressions lack sufficient precision for upcoming experiments.

Several approximate formulas balance accuracy and simplicity effectively.

Guidance is provided for choosing appropriate oscillation expressions in analysis tools.

Abstract

In this paper we review different expansions for neutrino oscillation probabilities in matter in the context of long-baseline neutrino experiments. We examine the accuracy and computational efficiency of different exact and approximate expressions. We find that many of the expressions used in the literature are not precise enough for the next generation of long-baseline experiments, but several of them are while maintaining comparable simplicity. The results of this paper can be used as guidance to both phenomenologists and experimentalists when implementing the various oscillation expressions into their analysis tools.

Figures

Tables

Table 1. Table 1: The various terms used throughout this paper. Natural units are assumed throughout unless otherwise mentioned.

$E$	Neutrino energy
$L$	Baseline
$a$	$2E\sqrt{2}{G}_F{n}_e$
$G_F$	Fermi’s constant
$n_e$	Electron # density
$s_{ij}$	$\sin {\theta }_{ij}$
$c_{ij}$	$\cos {\theta }_{ij}$
$\Delta m_{ij}^2$	$m_i^2-m_j^2$
$\Delta m_{ee}^2$	$c_{12}^2\Delta m_{31}^2+s_{12}^2\Delta m_{32}^2$
${\Delta }_{ij}$	$\Delta m_{ij}^{2}L/\left(4E\right)$
${\Delta }_x$	$xL/\left(4E\right)$
$ϵ$	$\Delta m_{21}^2/\Delta m_{ee}^2$
$\epsilon$	$\Delta m_{21}^2/\Delta m_{31}^2$
$J_r$	$s_{23}{c}_{23}{s}_{13}{c}_{13}^2{s}_{12}{c}_{12}$
$\hat{x}$	$x$ in matter
$\tilde{x}$	Approx. $x$ in matter

Table 2. Table 2: The expansion terms of each expression. Terms that are expansion parameters in the sense of Eq.  28 are denoted with a green check ( ✓ ), while terms that are not are denoted with a red cross ( × \color[rgb]{1,0,0}\times ). Note that Madrid refers also to the expression AJLOS(31) and FL which are all generally quite similar. AM refers to both AM 2 and AM 5/2 , and DMP refers to both DMP 0 and DMP 1 .

	$ϵ$ ($\epsilon$)	$s_{13}$	$a/\Delta m_{31}^2$
Madrid(like)	$\times$	$\times$	$\times$
AKT	✓	✓	✓
MP	✓	$\times$	$\times$
DMP	✓	✓	✓
AKS	$\times$	$\times$	$\times$
MF	✓	$\times$	$\times$
AJLOS(48)	✓	$\times$	$\times$
AM	$\times$	$\times$	$\times$

Table 3. Table 3: Neutrino oscillation parameters used. All of the values are within the 1 σ 𝜎 \sigma ranges of the best fit values obtained in the global fit in Ref.  [ 1 ] .

Parameter	Value
$\Delta m_{21}^2$	$7.5\times {10}^{-5}$ eV2
$\Delta m_{ee}^2$	$2.50\times {10}^{-3}$ eV2
${\sin }^2{\theta }_{12}$	0.32
${\sin }^2{\theta }_{23}$	0.55
${\sin }^2{\theta }_{13}$	0.022
$\delta$	-0.40$\pi$

Equations

$$P_{\mu e}=$$

$$+8J_r\frac{\Delta m_{21}^2}{a}\frac{\Delta m_{31}^2}{b}\sin {\Delta }_a\sin {\Delta }_b\cos \left(\delta +{\Delta }_{31}\right),$$

$$P_{\mu e}=$$

$$+8\frac{J_r}{c_{13}^2}\frac{\Delta m_{21}^2}{a}\frac{\Delta m_{31}^2}{b}\sin {\Delta }_a\sin {\Delta }_b\cos ({\Delta }_{31}+\delta ).$$

$$P_{\mu e}={\left|{\Delta }_{32}{e}^{i\delta }{s}_{13}{c}_{13}{s}_{23}\frac{e^{2i({\Delta }_{32}-{\Delta }_a)}-1}{{\Delta }_{32}-{\Delta }_a}-{\Delta }_{21}{s}_{12}{c}_{12}{c}_{23}\frac{e^{-2i{\Delta }_a}-1}{{\Delta }_a}\right|}^2.$$

$$\tan 2{\tilde{\theta }}_{12}$$

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

$${\lambda }_{1/p}$$

$${\lambda }_{2/3}$$

$$P_{\mu e}=$$

$$\times {\left(\frac{\Delta m_{ee}^2}{{\lambda }_{+}-{\lambda }_{-}}\right)}^2{\sin }^2\frac{({\lambda }_{+}-{\lambda }_{-})L}{4E}$$

$$+8ϵJ_r\frac{(\Delta m_{ee}^2{)}^3}{({\lambda }_{+}-{\lambda }_{-})({\lambda }_{+}-{\lambda }_0)({\lambda }_{-}-{\lambda }_0)}$$

$$\times \sin \frac{({\lambda }_{+}-{\lambda }_{-})L}{4E}\sin \frac{({\lambda }_{-}-{\lambda }_0)L}{4E}\cos \left(\delta +\frac{({\lambda }_{+}-{\lambda }_0)L}{4E}\right).$$

$${\lambda }_0$$

$$\mathrm{with}$$

$${\tilde{P}}_{\alpha \beta }(\Delta m_{ij}^2,{\theta }_{ij},\delta )=P_{\alpha \beta }(\Delta {\widetilde{m^2}}_{ij},{\tilde{\theta }}_{ij},\tilde{\delta }),$$

$${\sin }^2{\tilde{\theta }}_{13}$$

$$\Delta {\widetilde{m^2}}_{ee}$$

$${\sin }^2{\tilde{\theta }}_{12}$$

$$\mathrm{where}{a}_{12}=\frac{1}{2}(a+\Delta m_{ee}^2-\Delta {\widetilde{m^2}}_{ee}),$$

$$\Delta {\widetilde{m^2}}_{21}$$

$$\Delta {\widetilde{m^2}}_{31}$$

$${ϵ}^{\prime }\equiv ϵ\sin ({\tilde{\theta }}_{13}-{\theta }_{13})s_{12}{c}_{12}<0.015,$$

$$P_{\mu e}=4{\mathbb{C}}_{21}{\sin }^2{\Delta }_{21}+4{\mathbb{C}}_{31}{\sin }^2{\Delta }_{31}+4{\mathbb{C}}_{32}{\sin }^2{\Delta }_{32}+8\mathbb{D}\sin {\Delta }_{21}\sin {\Delta }_{31}\sin {\Delta }_{32},$$

$${\mathbb{C}}_{21}$$

$${\mathbb{C}}_{31}$$

$${\mathbb{C}}_{32}$$

$$\mathbb{D}$$

$${\mathbb{C}}_{21}^{(1)}$$

$${\mathbb{C}}_{31}^{(1)}$$

$${\mathbb{C}}_{32}^{(1)}$$

$${\mathbb{D}}^{(1)}$$

$$F_1$$

$$G_1$$

$$K_1$$

$$P_{\mu e}=$$

$$+4{\Delta }_{31}\sin (2{\Delta }_{31})c_{13}^2{s}_{13}{s}_{23}$$

$$\times \left\{-\frac{a}{\Delta m_{31}^2}{s}_{13}{s}_{23}\cos 2{\theta }_{13}+\epsilon s_{12}(-s_{13}{s}_{23}{s}_{12}+\cos \delta c_{23}{c}_{12})\right\}$$

$$-8J_r{\Delta }_{21}{\sin }^2{\Delta }_{31}\sin \delta .$$

$$P_0=$$

$$P_{\sin \delta }=$$

$$\times \left[\cos (C_{13}{\Delta }_{31})-\cos ({\Delta }_{31}+{\Delta }_a)\right],$$

$$P_{\cos \delta }=$$

$$\times \left[\sin (C_{13}{\Delta }_{31})-\sin ({\Delta }_{31}+{\Delta }_a)\right],$$

Full text

1. 1 Introduction

2. 2 Expansions under consideration

3. 2.1 The Madrid-like expressions

4. 2.1.1 The Madrid expression (2000)

5. 2.1.2 The AJLOS(31) expression (2004)

6. 2.1.3 The FL expression (2006)

7. 2.2 The AKT, MP and DMP expressions

8. 2.2.1 The AKT expression (2014)

9. 2.2.2 The MP expression (2015)

10. 2.2.3 The DMP expression (2016)

11. 2.3 Other expressions

12. 2.3.1 The AKS expression (1999)

13. 2.3.2 The MF expression (2001)

14. 2.3.3 The AJLOS(48) expression (2004)

15. 2.3.4 The AM expression (2011)

16. 2.4 Exact expressions

17. 2.4.1 The ZS expression (1988)

18. 2.4.2 Numerical diagonalization

19. 3 Comparison of expansions

20. 3.1 Expansion term

21. 3.2 Precision and speed

22. 4 Conclusions

1 Introduction

Over the last two decades neutrino oscillation measurements have become increasingly precise and are now entering the precision era. Most of the current data coming from experiments using neutrinos from the sun, reactors, the atmosphere and particle accelerators can be described in terms of three-neutrino oscillations, which depend on the six oscillation parameters: two mass splittings $\Delta m_{31}^{2}$, $\Delta m_{21}^{2}$, three mixing angles $\theta_{12}$, $\theta_{13}$ and $\theta_{23}$ and a CP-violating phase $\delta$. Many of these parameters are measured rather well as of now [1]. However, there are some remaining unknowns, as for example the value of the CP-phase $\delta$, the octant of the atmospheric angle ($\sin^{2}\theta_{23}<0.5$ or $\sin^{2}\theta_{23}>0.5$) and the neutrino mass ordering ($\Delta m_{31}^{2}>0$ or $\Delta m_{31}^{2}<0$). Note however, that from combining oscillation data with data from cosmological observations a 3.5$\sigma$ preference for normal ordering can be obtained [2, 3]. There are also some anomalies, which might suggest the existence of a fourth neutrino, see Refs. [4, 5] for the current status. However, here we will focus only on the case of standard three-neutrino oscillations.

To obtain functions for the oscillation probabilities one has to solve the corresponding Schrödinger equation. While this is easily done in the case of vacuum oscillations (i.e. a free Hamiltonian), it becomes very difficult to do it in the presence of matter due to alterations of the oscillation patterns due to the Wolfenstein matter effect [6, 7, 8, 9]. However, if one assumes a sufficiently constant matter profile in the trajectory of the neutrino analytic solutions can be found [10]. These expressions are very complicated unfortunately and do not permit for a deeper understanding of the phenomenology of three-neutrino oscillations due to the presence of the $\cos\left(\frac{1}{3}\cos^{-1}\cdots\right)$ term shown later in Eq. 22.

Therefore to obtain better insights one may try to find simpler analytic expansions around naturally appearing small parameters. Some commonly used small parameters are the matter potential $a/\Delta m_{31}^{2}$ [11], $\sin\theta_{13}$ or $\sin^{2}\theta_{13}$ [12, 13], and the ratio of mass splittings $\Delta m_{21}^{2}/\Delta m_{31}^{2}$ or $\Delta m_{21}^{2}/\Delta m_{ee}^{2}$ [11, 14, 15, 16, 17, 18, 19], where $\Delta m_{ee}^{2}\equiv\cos^{2}\theta_{12}\Delta m_{31}^{2}+\sin^{2}\theta_{12}\Delta m_{32}^{2}$ [20, 21]. In this paper we analyze, with an eye for both precision and computational speed, different expansions and show how accurate they are and how they have aged as the measurement of the oscillation parameters has evolved over the past twenty years. It is clear that simplicity is an important trait for an approximate expression. While simplicity may be somewhat in the eye of the beholder, we use computational speed as a rough proxy for simplicity.

We focus on the DUNE [22] experimental conditions of $L=1300$ km and Earth density of $\rho=3$ g/cm3. We take the matter density to be constant which is a good approximation since effects due to the variability of the density of the Earth are beyond the sensitivity of DUNE [23]. Our results will also generally apply to other long-baseline experiments such as NOVA [24], T2K/HK [25, 26]. We also discuss the second oscillation maximum which is relevant for T2HKK [27] and ESSnuSB [28].

In the next section we review the various expansions we study. In section 3 we compare how well they work comparing them among each other and also in comparison to the exact solution. We also check how fast the expansions can be computed in comparison to the exact analytic solution and to the numerical diagonalization process. Finally, section 4 contains our conclusions.

2 Expansions under consideration

In this section we present the expansions we will compare to the exact expression. We categorize them into three groups based on their forms. The first is the “Madrid-like” group named for the common city for which one expression was first written down. Other very similar expressions followed and they will be grouped together accordingly. The next group is the AKT, MP and DMP group. This is a series of works that performs two flavor rotations and then perturbation theory. The final group contains the remaining expressions. Expressions generally drop terms proportional to various smallness parameters including the ratio of $\Delta m^{2}$’s, $s_{13}$, or the matter potential.

For historical reasons (i.e. neutrino factory [29]) many of these expansions have been performed in the channel $\nu_{e}\rightarrow\nu_{\mu}$. However, in the context of long-baseline accelerators the most important channel is $\nu_{\mu}\rightarrow\nu_{e}$. Therefore we will present our results in this channel. They are related to each other through the T-relation $P(\nu_{\mu}\rightarrow\nu_{e};\delta)=P(\nu_{e}\rightarrow\nu_{\mu};-\delta)$, just switching the sign of the CP-phase555Note that sending $\delta\to-\delta$ is equivalent to sending $L\to-L$ under the assumption of CPT invariance.. If one is interested in antineutrinos, namely $P(\overline{\nu}_{\mu}\rightarrow\overline{\nu}_{e};E)=P(\nu_{\mu}\rightarrow\nu_{e};-E)$, one only has to switch the sign of the neutrino energy. We will focus on neutrinos, but our results generally apply to antineutrinos as well. Note that matter effects are not very important in the disappearance channels $\nu_{\mu}\rightarrow\nu_{\mu}$ and therefore focusing only on the appearance channel will also not affect the main message of this paper.

While the choice of notation does not affect the precision of these formulas, it does affect their general clarity and overall usefulness. To this end, we have chosen to use uniform notation throughout this article as much as possible with the various terms defined in Tab. 1. In doing so we have made several simplifying manipulations, in each case maintaining the exact same mathematical expression. The relationship between the notation used here and the original notation used is mentioned below whenever applicable. While these definitions represent fairly commonly used definitions in the literature, some differ by factors of two or other slight changes, so care is required when making comparisons.

2.1 The Madrid-like expressions

In this subsection we list a few nearly identical expressions and discuss their similarities and differences.

2.1.1 The Madrid expression (2000)

This expression was derived in Ref. [12] (by Cervera, Donini, Gavela, Gomez Cádenas, Hernández, Mena, and Rigolin, Madrid hereafter) and can be written666In Ref. [12] $\tilde{J}\equiv 8J_{r}$ is used instead, the definition of $\Delta_{ij}$ differs by a factor of $L/2$, and terms $A\equiv a/(2E)$ and term $B\equiv|b|/(2E)$ are used. While Ref. [12] defines their $b$-like parameter with absolute value signs, we note that they are not necessary since the probability is even in $b$.

[TABLE]

where $b\equiv a-\Delta m^{2}_{31}$.

The form of this expressions suggests the square of two summed amplitudes. It is not exactly such a sum due to an extra factor of $c_{13}$ in the interference term which provides the correct CP-violating term in vacuum. Writing Eq. 1 as the sum of two amplitudes has been examined in various forms in Refs. [30, 16, 31, 32].

There are two ways to correct this. The first is to drop one of the $c_{13}$’s in the Jarlskog invariant ($J_{r}$). Alternatively, if we add in a factor of $c_{13}$ to one of the amplitudes, we recover the Jarlskog in vacuum correctly while still writing the expression as the sum of two amplitudes. The natural place to put it is on the $a$ (21) term [30, 31, 32, 33] as this reproduces the vacuum expression exactly. We note that this provides a negligible change to the precision of the equation as we are correcting an already subleading term (the solar term) by a small amount. On the other hand we can add the $c_{13}$ term to the $b$ (31) term. Doing so improves the precision of the Madrid expression for neutrinos by about an order of magnitude at the first oscillation maximum, although this effect is due to a lucky cancellation for the parameters used. The improvement is more modest elsewhere, and for antineutrinos the precision is a bit worse than the Madrid expression. As such we do not include such an expression in our subsequent analyses.

In addition, considering the previously identified importance of $\Delta m^{2}_{ee}$ [21] and the fact that there is no reason to use $\Delta m^{2}_{31}$ or $\Delta m^{2}_{32}$ unless both are treated separately, we have also examined how Eq. 1 performs with $\Delta m^{2}_{31}\to\Delta m^{2}_{ee}$. We find that this change results in somewhat better performance in some cases (modest improvement at $E\gtrsim$ few GeV and considerable improvement at and below the second maximum), in the region of interest for DUNE around a few GeV the performance is essentially the same as the Madrid expression.

2.1.2 The AJLOS(31) expression (2004)

In Ref. [15] (by Akhmedov, Johansson, Lindner, Ohlsson, and Schwetz, AJLOS hereafter) several expressions are introduced each with different expansion parameters. We label them by the equation numbers in the original paper777In Ref. [15] $\alpha=\varepsilon$, $A=a/\Delta m^{2}_{31}$, and $\Delta=\Delta_{31}$ are used. We refer to the equation numbers in the version on the arXiv, not in the journal..

The first expression (#31) drops higher order terms proportional to $\varepsilon$ and $s_{13}$.

[TABLE]

We note that up to a factor of $c_{13}^{2}$ in each the second and third term this expression is otherwise identical to the Madrid expression in Eq. 1.

2.1.3 The FL expression (2006)

Also the authors of Ref. [16] (by Friedland and Lunardini, FL hereafter) write the probability as the sum of two amplitudes. Using our notation888In Ref. [16] $\Delta_{1}=2(\Delta_{32}-\Delta_{a})$, $\Delta_{2}=-2\Delta_{a}$, $G_{1}=\Delta_{32}\sin 2\theta_{13}e^{i\delta}$, and $G_{2}=-\Delta_{21}\sin 2\theta_{12}$ were used. they obtain

[TABLE]

Note that this expression was derived actually in the context of Non-Standard neutrino Interactions (NSI) and that it reduces to eq. 3 once all the NSI parameters are set to zero. This expression is identical to Eq. 1 up to using $\Delta m^{2}_{32}$ instead of $\Delta m^{2}_{31}$ and factors of $c_{13}$.

2.2 The AKT, MP and DMP expressions

In [34, 17, 18] a different technique was used. Two-flavor rotations were performed to simply diagonalize the Hamiltonian by focusing on the largest off-diagonal terms first. This means that all channels ($\nu_{\alpha}\to\nu_{\beta}$ for $\alpha,\beta\in\{e,\mu,\tau\}$) are handled simultaneously. In AKT the focus is on the vacuum mass eigenstate basis whereas MP and DMP focus on the flavor basis. This choice affects the order of the two, two-flavor rotations and hence the precision. Roughly speaking, AKT performs a 12 rotation followed by 13 rotation whereas MP and DMP first perform a 13 rotation followed by a 12 rotation as sketched below.

2.2.1 The AKT expression (2014)

The authors (Agarwalla, Kao, Takeuchi) of Ref. [34] (AKT hereafter) perform the rotations going from largest contribution in the Hamiltonian to smallest in the mass basis. They begin with the 12 rotation, followed by the 23 rotation and then after commuting with the 12 rotation, they absorb this 23 rotation into a 13 rotation.

Using this approach the effective mixing angles can be written as

[TABLE]

Note that $\delta$ and $\theta_{23}$ are treated as constant in matter. The eigenvalues of the Hamiltonian are now given by

[TABLE]

From this expressions we obtain the mass splittings in matter simply via the relation $\Delta\widetilde{m^{2}}_{ij}=\lambda_{i}-\lambda_{j}$. The oscillation probability can now be obtained by replacing the vacuum parameters with the matter parameters in the vacuum oscillation probability, see Eqs. 11 and 12 below. While $\Delta m^{2}_{ee}$ was not explicitly used in Ref. [34], it appears in several places nonetheless and we made the substitution here for simplicity.

2.2.2 The MP expression (2015)

After one rotation in the 13 sector, we have an expression that is an expansion in $\epsilon c_{12}s_{12}$. This is the expression in Ref. [17] (by Minakata and Parke, MP hereafter). The expression is

[TABLE]

where999In Ref. [17] $\Delta m^{2}_{\rm ren}=\Delta m^{2}_{ee}$ is used and the definition of the eigenvalues is shifted by $\epsilon\Delta m^{2}_{ee}s_{12}^{2}$.

[TABLE]

This expression is very accurate except near the solar resonance when ${\Delta_{21}>1}$.

2.2.3 The DMP expression (2016)

To address the solar resonance, an additional rotation was performed in Ref. [18] (by Denton, Minakata, and Parke, DMP hereafter). The order of the rotations (12 then 13 after an initial constant 23 rotation) is chosen to diagonalize the largest remaining off-diagonal term at each step. This procedure also removes both level crossings with the minimal number of new angles. After these two rotations, perturbation theory is now possible everywhere, although the zeroth order expression (DMP0) is sufficiently precise for future long-baseline experiments.

Here, as in AKT in section 2.2.1 above, the authors do not derive formulas for the oscillation probabilities directly, but rather for the oscillation parameters in matter and then write the probability

[TABLE]

where $\widetilde{x}$ refers to the approximate expression for the quantity $x$ evaluated in matter. That is, to an excellent approximation the oscillation probability in matter has the same form as the expression in vacuum with only $\Delta m^{2}_{21}$, $\Delta m^{2}_{31}$, $\theta_{12}$, and $\theta_{13}$ replaced by their approximate matter equivalents ($\widetilde{\theta}_{23}$ and $\widetilde{\delta}$ are roughly constant in matter). The authors obtain to zeroth order the following expressions101010In Ref. [18] $\phi=\widetilde{\theta}_{13}$, ${\psi}=\widetilde{\theta}_{12}$, $\lambda_{i}=\widetilde{m^{2}}_{i}$ were used. The current notation is consistent with Ref. [35]. The following expression from [18] is also useful: $\cos^{2}(\theta_{13}-\widetilde{\theta}_{13})=c_{\widetilde{13}}^{2}c_{13}^{2}+s_{\widetilde{13}}^{2}s_{13}^{2}+\sin 2\widetilde{\theta}_{13}c_{13}s_{13}=(\Delta\widetilde{m^{2}}_{ee}+\Delta m^{2}_{ee}-a\cos 2\theta_{13})/(2\Delta\widetilde{m^{2}}_{ee})$.,

[TABLE]

Note that we have $\tilde{\delta}=\delta$ and $\tilde{\theta}_{23}=\theta_{23}$ at zeroth order. The expansion parameter is

[TABLE]

and is zero in vacuum confirming that this expression returns the exact expression in vacuum.

We recall the vacuum expression here, which we write in the following form,

[TABLE]

where the coefficients are,

[TABLE]

DMP0 is then Eqs. 11, 12 where the vacuum parameters are replaced with the approximate matter ones given in Eq. 9.

Also note that $\Delta\widetilde{m^{2}}_{ee}$ (which was further explored in [35]) is the same as $\lambda_{+}-\lambda_{-}$ in Eq. 7 above from the MP formula in 2.2.2. Successive orders of precision can also be calculated by following perturbation theory in a straightforward fashion as is done through second order in [18] with compact expressions provided through first order or by correcting the mixing angles directly [36]. Here we focus on zeroth and first orders only (DMP0 and DMP1 hereafter respectively) as they are already extremely precise. Successive orders add $\sim 2.5$ additional orders of magnitude of precision if desired and expressions through second order exist in [18].

The first order corrections to the coefficients from Eq. 12 are given by the following expressions,

[TABLE]

where $\epsilon^{\prime}\equiv\epsilon\sin(\widetilde{\theta}_{13}-\theta_{13})s_{12}c_{12}$ as shown in Eq. 10 above, and

[TABLE]

and the $F_{2},G_{2},K_{2}$ expressions are related to the above by making the transformation $c_{\widetilde{12}}^{2}\leftrightarrow s_{\widetilde{12}}^{2}$, $c_{\widetilde{12}}s_{\widetilde{12}}\to-c_{\widetilde{12}}s_{\widetilde{12}}$, and $m_{1}\leftrightarrow m_{2}$. This correction is DMP1. Note that the expressions of Eq. 12 are also invariant under this transformation [18].

2.3 Other expressions

Here we list other expressions in the literature that do not fall into the above two categories.

2.3.1 The AKS expression (1999)

The oldest expression we consider in this paper is the one derived in Ref. [11] (by Arafune, Koike, and Sato, AKS hereafter). Here the authors obtain

[TABLE]

2.3.2 The MF expression (2001)

In Ref. [14] (by Freund, MF hereafter) the author separates the oscillation probability111111We use Eq. 36 of Ref. [14] since Eq. 38 behaves poorly for large values of $\theta_{13}$. in sub-terms, $P_{\mu e}=P_{0}+P_{\sin\delta}+P_{\cos\delta}+P_{1}+P_{2}+P_{3}$. These terms are given by121212In Ref. [14] $\hat{C}=C_{13}$ is used and there is a typo wherein $\cos\theta_{13}^{2}$ is written instead of $\cos^{2}\theta_{13}$.

[TABLE]

The expression in Eq. 16 is similar to the AJLOS(48) expression in Eq. 17 below.

2.3.3 The AJLOS(48) expression (2004)

The second AJLOS (Ref. [15]) expression (#48) focuses only on $\varepsilon$ as a smallness parameter. The first two orders, $P_{\mu e}=P_{\mu e}^{(0)}+\varepsilon P_{\mu e}^{(1)}$, are

[TABLE]

where $C_{13}\equiv\sqrt{\sin^{2}2\theta_{13}+(a/\Delta m^{2}_{31}-\cos 2\theta_{13})^{2}}$ (this factor is also used in MF in the previous subsection, 2.3.2). We note, $C_{13}$ is equivalent to $\Delta\widetilde{m^{2}}_{ee}/\Delta m^{2}_{ee}$ from DMP [18, 35] after the change $\Delta m^{2}_{31}\to\Delta m^{2}_{ee}$.

There is also a third expression in the AJLOS paper, Eq. 66, but this is designed for the solar sector and does quite poorly for long-baseline oscillations which are dominated by the atmospheric term. For this reason we do not consider it here.

2.3.4 The AM expression (2011)

We study the expression obtained by the authors of Ref. [13] (by Asano and Minakata, AM hereafter). Here the authors use as expansion parameters $s_{13}\simeq\sqrt{\Delta m_{21}^{2}/\Delta m_{31}^{2}}$. The authors divide their expressions in powers of $s_{13}$, $P_{\mu e}^{(0)}+P_{\mu e}^{(1)}+P_{\mu e}^{(3/2)}+P_{\mu e}^{(2)}+P_{\mu e}^{(5/2)}+\ldots$ Each superscript refers to the order of $s_{13}^{2}$ and $\Delta m^{2}_{21}/\Delta m^{2}_{31}$. For our channel of interest they obtain131313In Ref. [13] $r_{\Delta}=\varepsilon$ and $\Delta=2\Delta_{31}/L$ are used.

[TABLE]

where $r_{A}\equiv a/\Delta m_{31}^{2}$. The $5/2$ order term is given by

[TABLE]

We will consider the precision and speed both of the expression through second order (AM2 hereafter) and through $5/2$ order (AM5/2 hereafter).

2.4 Exact expressions

For completeness, we discuss two different means of exactly calculating the oscillation probabilities. The first uses the analytic solution to a cubic equation and the second involves numerically diagonalizing the Hamiltonian. We have verified that these are equivalent up to numerical precision, $\sim 10^{-13}$.

Several pieces of software designed to solve neutrino oscillations in matter also exist in the literature [7, 37, 38]. While these may offer modest improvements in speed over off-the-shelf linear algebra packages, we have verified that they still do not compete with analytic expressions in terms of speed. In fact, several of them are the same, in part or in full, as the ZS expression in the next section. NuSquids [38] is a bit different than the other options in that once it has solved the differential equation for a given matter profile and baseline, extracting the probability (or the flux more generally) for a given energy is fairly efficient. It starts to become computational efficient for long-baseline at $\gtrsim 1000$ energy computations [39].

2.4.1 The ZS expression (1988)

It is possible to solve the characteristic polynomial of the Hamiltonian in matter directly. Using this approach one can express the eigenvalues in matter in terms of the vacuum parameters which involves solving a completely general cubic equation141414The original solution of the cubic equation was from [40] based on work by Scipione del Ferro and Niccolò Fontana Tartaglia in the sixteenth century.. This was done in Ref. [10] (by Zaglauer and Schwarzer, ZS hereafter), where the authors obtain

[TABLE]

The mass splittings in matter are then given by

[TABLE]

where $\widehat{x}$ refers to the given quantity in matter and151515Note, that in Ref. [10] $D=a$ and there are two typos: the root in the denominator of $S$ should be to the $3/2$ power not $1/3$ and the numerator of $e^{-i\widehat{\delta}}$ should have a factor of $s_{23}c_{23}$ instead of $s_{23}c_{23}^{2}$.

[TABLE]

The parameter $S$ is the one mentioned in the introduction that has $\cos\left(\frac{1}{3}\cos^{-1}\cdots\right)$ which ends up in all the final parameters: mass squared differences, mixing angles, and the CP-phase. The mixing angles and the CP-phase are

[TABLE]

where

[TABLE]

As in the case of DMP the oscillation probabilities in matter can now be obtained by simply replacing the vacuum parameters with the matter parameters in the vacuum oscillation probability from Eq. 11.

2.4.2 Numerical diagonalization

It is also possible to diagonalize the Hamiltonian numerically. The Schrödinger equation in matter can be written in matrix form as

[TABLE]

where

[TABLE]

Here $\psi_{\alpha\beta}(t)=\langle\nu_{\beta}|\nu_{\alpha}(t)\rangle$ is the oscillation amplitude, $U=U_{23}(\theta_{23})\times U_{13}(\theta_{13},\delta)$ $\times U_{12}(\theta_{12})$ is the PMNS-matrix and $\mathbb{M}^{2}=\operatorname{diag}(0,\Delta m_{21}^{2},\Delta m_{31}^{2})$. The matter potential is given by $\mathbb{A}=\operatorname{diag}(a,0,0)$, where as always in this work $a$ is constant. In principle this equation can be solved easily by writing $\mathcal{H}=RDR^{\dagger}$, where $D=\operatorname{diag}(d_{1},d_{2},d_{3})$ is the diagonal matrix containing the eigenvalues $d_{i}$ of $\mathcal{H}$, and $R$ is the diagonalization matrix of eigenvectors. Then the solution is easily found to be

[TABLE]

This diaginalization can be performed numerically using for example the packages Eigen [41] or HEigensystem [42, 43], among others. We will refer to this method as Diag from here on.

3 Comparison of expansions

We now compare the usefulness of each expression. First we describe the behavior of the various expressions under several useful limits. Next, we define two metrics: precision and speed/simplicity. Precision can be quantified as either the error or fractional error between a given expression and the exact expression. We focus on fractional error since that is more relevant for experiments although it can become misleading when the probability is small or goes to zero such as in the high energy limit. While the simplicity of an expression is a somewhat subjective metric, the computational efficiency is somewhat more scientific and quantitative.

3.1 Expansion term

Each approximate expression is an expansion in one or more parameters. In order to clearly show how each expression behaves, Tab. 2 shows which parameters each formula is expanded in. In order to qualify as an expansion parameter we require that the probability recovers the exact (to all orders) expression as that parameter goes to zero. That is, $x$ is an expansion parameter if and only if

[TABLE]

We note that as many expressions drop higher order terms of more than one parameter at a time, it is quite common for expressions to not be true expansions in the sense of Eq. 28 in that all of the parameters that were treated as small numbers simultaneously.

We find that DMP (at any order) as well as AKT are the only expressions that are an expansion in $s_{13}$ and the matter potential. Also, while several expressions are expansions in $\epsilon$ or $\varepsilon$ (including DMP and AKT), some are not an expansion in $\epsilon$ or $\varepsilon$ either despite treating $\varepsilon$ as a smallness parameter, such as the Madrid-like expressions. In addition to the parameters listed in the table we note that none of the expression are exact as $L\to 0$, the so-called vacuum mimicking regime [44].

3.2 Precision and speed

In this section we compare the different expressions. As a benchmark point we use the standard oscillation parameters in Tab. 3. We take $\delta$ to be slightly off-maximal to avoid any unintentional cancellations and to require both $\sin\delta$ and $\cos\delta$ terms to be correct. We choose as benchmark baseline $L=1300$ km and density $\rho=3$ g/cm3, the configuration for the DUNE experiment [45, 46, 47] although our results are applicable to any current or future long-baseline experiment, including those focusing on the second oscillation maximum.

The probabilities for selected expressions are shown in Fig. 1 for both neutrinos (left) and antineutrinos (right). They should be compared to the exact curves given by ZS. To compare the precision of the formulas in a more quantitative way we show $\frac{|P_{\text{test}}-P|}{P}=\frac{|\Delta P|}{P}$, where $P=P_{\text{ZS}}$ is the exact formula. The result is shown in Fig. 2. As one can see for low energies DMP gives the best results, while for $E>2$ GeV the most precise result is AM5/2. The precision of DMP0 and AKT is the worst at the atmospheric resonance ($\sim 11$ GeV) before leveling off, although the probability is approaching zero thus this region is less relevant experimentally. We also show the precision of the several Madrid-like expressions in Fig. 3. The precision of the remaining expressions can be found in Fig. 4. Note that the sharp dips are not representative of improved precision, rather they represent a crossing between the exact and approximate expressions. In addition, the peaks in the errors are when the oscillation probability, and thus the denominator, goes to zero.

In order to more clearly compare the precision of each expression, we show in Fig. 5 the precision with which each expression reconstructs the first and second oscillation maxima. We focus on the oscillation maxima because the heights (probability) of the maxima are an important test for CP violation [48] and the locations (energy) of the maxima are an important test for the atmospheric mass splitting [49]. The horizontal line at 1% is to guide the eye. Since DUNE and other next generation long-baseline experiments are aiming to reach near the percent level in precision, we cannot introduce theoretical errors larger than 1%. We have also verified that these results are generally robust under changes to the oscillation parameters, although for certain specific values of say the CP-violating phase some of the fairly precise expressions may appear to perform much better if there is a crossing between the approximate and exact expressions at the oscillation maximum.

We also measured the computational time on a single core i7 processor using c++ with the gnu v7.3.0 compiler161616We have performed the same test using Fortran with the gfortran compiler, obtaining qualitatively the same results.. Our result can be found in Fig. 6 along with the precision at the first oscillation maximum. The green (red) dots are those for which the probability at the first maximum is better (worse) than 1%. We also measure the computational time for the exact vacuum expression from Eq. 11 for comparison. All of the approximate expressions and ZS are faster than the diagonalization by about an order of magnitude or more, depending on the expansion of interest. The DMP expressions have the best precision among the approximate expressions by a considerable amount.

When computational speed is the primary, if 1% precision is required then MP is the best expression, although it is right at $\sim 1\%$ and has much larger errors at $\sim 1$ GeV. After MP, the next simplest expression that is also better than 1% is DMP0 which is precise at the per mille level or better. The software used to perform these tests Nu-Pert-Compare is publicly available at github.com/PeterDenton/Nu-Pert-Compare [50].

4 Conclusions

With the advent of next generation long-baseline neutrino experiments, neutrino oscillation physics will truly enter the precision era. Over the years, various approximate expressions describing oscillation probabilities have been written down. In this paper we have normalized the notation as much as possible and categorized similar expressions. We then directly compared them with the exact solution with an eye for precision and speed where the latter provides a rough proxy for simplicity.

We have found that the DMP expressions from Ref. [18] (Denton, Minakata, Parke) are the most precise expressions available in most cases, even at zeroth order, with the option for considerably improved precision at higher orders. In terms of simplicity they are comparable to any of the other expressions available as shown in Fig. 6. It is interesting to see within the expression from Ref. [13] (Asano and Minakata) that while the addition of the $5/2$ term adds considerably to the precision making it about as precise as any expression, leads to a considerable loss in simplicity and speed as can be seen from Eq. 19 and Fig. 6. In the same vein we can see in Fig. 6 that DMP1 adds much more precision (close to two orders of magnitude) than the already quite precise DMP0 with only a modest cost in complexity due to the compact form of the first order correction.

Finally, we note that most of the most precise expressions, AKT and the DMP expressions, all naturally use $\Delta m^{2}_{ee}$ which is the $\nu_{e}$ average of $\Delta m^{2}_{31}$ and $\Delta m^{2}_{32}$. This was not “forced” or put in by hand, rather it naturally appears. This suggests that whenever one atmospheric mass splitting is required, the correct expression to use is $\Delta m^{2}_{ee}$ since there is no reason to prefer one of $\Delta m^{2}_{31}$ or $\Delta m^{2}_{32}$ over the other.

Acknowledgements

We thank Chris Weaver and Carlos Argüelles for helpful discussions. GB acknowledges support from the MEC and FEDER (EC) Grant SEV-2014-0398, FIS2015-72245-EXP, and FPA2014-54459 and the Generalitat Valenciana under grant PROMETEOII/2017/033, also partial support along with SJP from the European Union FP7 ITN INVISIBLES MSCA PITN-GA-2011-289442 and InvisiblesPlus (RISE) H2020-MSCA-RISE-2015-690575. PBD acknowledges the United States Department of Energy under Grant Contract desc0012704. This manuscript has been authored (SJP) 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. SJP thanks IFT of Madrid for wonderful hospitality during part of this work. This project has received funding/support from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement Nos 690575 and 674896. CAT is supported by the Spanish grants FPA2017-90566-REDC (Red Consolider MultiDark), FPA2017-85216-P and SEV-2014-0398 (MINECO/AEI/FEDER, UE), as well as PROMETEO/2018/165 (Generalitat Valenciana) and the FPI fellowship BES-2015-073593. CAT acknowledges the hospitality of the Fermilab Theoretical Physics Department, where part of this work was done.