Fast Exact Algorithm for Neutrino Oscillation in Constant Matter Density

A recently published method for solving the neutrino evolution equation with constant matter density is further refined and used to lay out an exact algorithm for computing oscillation probabilities, which is moderately faster than previous methods when looping through neutrinos of different energies. In particular, the three examples of $\overset{\scriptscriptstyle{(-)}}{\nu}_e$ survival, $\overset{\scriptscriptstyle{(-)}}{\nu}_\mu$ survival and $\overset{\scriptscriptstyle{(-)}}{\nu}_e$ appearance probabilities are written in terms of mixing angles, mass differences and matter electron density. A program based on this new method is found to be roughly twice as fast as, and in agreement with, the leading GLoBES package. Furthermore, the behaviour of all relevant effective parameters is sketched out in terms of a range of neutrino energies, or matter electron densities. For instance, the $\overset{\scriptscriptstyle{(-)}}{\nu}_e$ survival probability in constant matter density is found to have no dependence on the mixing angle $\theta_{23}$ or the CP-violating phase $\delta_{13}$.


I. INTRODUCTION
A. The problem Neutrinos are produced in charged current (CC) interactions in pure flavour states |ν α ⟩, α ∈ {e, µ, τ }, which are composed of a superposition of the mass states |ν k ⟩, k ∈ {1, 2, 3} where it is assumed the mass differences have a negligible impact on kinematics -not favouring some mass states' formation over others'.U αk is the unitary PMNS matrix, and the normalisation conditions are ⟨ν k |ν j ⟩ = δ kj , so that ⟨ν α |ν β ⟩ = δ αβ .In a vacuum, it is the mass states that are eigenstates of the free Hamiltonian ( Ĥ0 ), and so whose evolution can be computed where P µ 0 is the (free) spacetime translation operator ( P 0 0 = Ĥ0 ) and p µ k is the four-momentum of mass state k.Assuming plane-wave solutions, this is solved with where |ν k ⟩ ≡ |ν k (0)⟩, which can be further simplified via the ultra-relativistic approximation where m k is the mass of ν k , E its energy, and L is the distance travelled.One can thus arrive at the familiar vacuum transition/survival probability by defining P να→ν β (L, E) ≡ * j.page@sussex.ac.uk ⟨ν β (L)|ν α ⟩ 2 , so that which is governed by the mass differences ∆m 2 kj ≡ m 2 k − m 2 j and the PMNS matrix.
When taking into account matter effects, the process is not so straightforward.Recall that previously, only the free Hamiltonian was used for spacetime translations P µ 0 (eqn 2).Strictly speaking, one should use the full Hamiltonian, which contains the CC interaction terms where Ĥeff I is the average effective interaction Hamiltonian, caused by coherent forward elastic scattering in matter.Incoherent scattering in matter is exceedingly unlikely (σ ∼ G 2 F s), so one need only consider coherent elastic scattering [3].These are described by the effective Fermi theory terms where G F is Fermi's constant, g ψ V and g ψ A are the vector and axial components of the coupling constants for the associated fermion ψ, and all the other symbols have their usual meaning.Averaging over the particles in the matter medium, one can show that [3] Ĥeff where N e and N n are the electron and neutron densities in matter respectively.The ± is + for neutrinos and − for antineutrinos [3].Now, notice that the interaction Hamiltonian acts on the flavour eigenstates, while the free spacetime translation operator acts on the mass eigenstates, creating a complicated differential equation.To deal with this, the transition amplitude ψ αβ (x) = ⟨ν β |ν α (x)⟩ is used, so that each part of the translation operator may act on its corresponding eigenstates.The neutron density term can be factored out since it affects all flavours equally, along with other common factors, and one arrives at the ODE [3] where x is now the spacial coordinate in the direction of propagation, and Notice again that it is the mass differences that come into play, and the difference with the vacuum ODE (eqn 2) is isolated to the A matrix, which vanishes for N e = 0.

B. Extant Solutions and Approximations
There are many ways to approach this problem, which have already been written about extensively.A sample is shown here for context.First, the mass hierarchy can be exploited to separate out the contributions of each mass difference, depending on one's setup.This can "freeze out" one of the three transition amplitudes, and if one assumes a constant matter density, can be reduced to effective two-neutrino mixing again.These expressions are very helpful to gain a qualitative understanding of the processes, such as resonances for specific values of A CC that lead to a strong MSW effect.See [3] or [4] for more details.
However, if the background matter density varies, or one mass difference does not totally dominate, one must keep track of all the transition amplitudes.[4] and [5] use diagonalisation to compare with the vacuum case and determine the effective mass difference and mixing angles in matter.This is the standard approach, including numerical techniques: numerically diagonalise at each iteration of the evolution equation [6,7].Meanwhile, [8] uses Lagrange's formula to determine the evolution operator where E n are the three eigenvalues of H F , with rather involved expressions provided for constant matter density.The expression is written here explicitly due its similarity with what will be shown later on.For its part, [9] uses the Cayley-Hamilton theorem to decompose the evolution operator into a linear combination of second order polynomials of mixing matrices, with analytic expressions for these matrices as well as their coefficients.Recently, [10] used the eigenvectoreigenvalue identity to derive relatively simple formulae for the effective mixing angles and CP violating phase, along with perturbative approximations of these and mass differences, which can be used in the vacuum expressions as normal.A summary of many of these exact and approximate techniques can be found in [11], along with very useful accuracy and speed comparisons in the context of long baseline ν e appearance experiments.[12] provides an elegant generalisation to (3 + N ) neutrino flavours, and a generic matter potential that can include non-standard interactions (NSI).The effects of the sterile neutrinos and NSI are studied both together and independently.These last two papers also make their codes available via public GitHub repositories, referenced therein.
Lastly, two recent papers [1,2] compute a general form for the evolution operator (assuming constant matter density) in terms of a Gell-Mann basis and structure constants, using methods from [9].The first paper then derives perturbative expansions for particular electron density profiles in the Earth, while the second formulates a general method to compute the oscillation probability of any general time independent Hamiltonian for two or three active neutrino flavours.The initial approach of these will be followed in the next section, with some small tweaks, before being taken in another more specific direction than the original papers.These deviations will be highlighted throughout the derivation, and some extra details added for clarity.An efficient algorithm will then be constructed, which splits the computation up into two parts, so that applications which compute oscillation probabilities for a great many neutrinos need only perform the first part once as an initialisation, saving time.Some specific examples are the provided and compared to existing numerical computations by the GLoBES package and others.The calculation and behaviour of effective parameters will then briefly be covered.

II. SOLVING THE DIFFERENTIAL EQUATION A. The Evolution Operator
This subsection largely follows Bushra Shafaq and Faisal Akram's method in [1].First, the traceless effective Hamiltonian H is defined, since the trace acts on all flavours equally and so does not contribute to mixing The evolution equation is thus where contrary to Shafaq and Akram's paper, the 1/2E factor is kept separate from H. Assuming constant matter density, this is solved with where U (x) is the evolution operator.Now, these can be decomposed using the property that the Gell-Mann matrices (λ i , i ∈ {1, ..., 8}) and the identity matrix form a complete orthogonal basis for 3 × 3 complex matrices.H is traceless, so it does not need the identity matrix where from now on repeated dummy indices imply summation.These equations are derived from the Gell-Mann matrix identities tr λ i λ j = 2δ ij and tr λ i = 0. Now, some useful general results in linear algebra will be used: for a matrix A, with eigenvalues For u i , first note that from and so using the previous identities one can show All that is needed now are expressions for the eigenvalues E n of H.The parametric equation of a 3×3 matrix A with eigenvalues λ is so that for the traceless H, where d ijk are the symmetric structure constants of the Gell-Mann matrices The last relation between the determinant and structure constants in (eqn 23) can be derived by first multiplying the structure constant definition (eqn 24) (second equation) by h i h j h k (and summing over these indices as normal), to find tr Then from the definition of the determinant of a 3 × 3 matrix the Levi-Civita identity b3 (the index position is irrelevent here), and recalling that the Gell-Mann matrices are traceless (λ i aa = 0), one can find and thus Meanwhile, taking the derivative of the parametric equation (eqn 23) w.r.t h i gives the needed expression while different solutions to the parametric equation are used here compared to the original paper (30) These are the solutions to a depressed cubic equation for real solutions, which must be real since H is Hermitian (one can also check that a 0 , a 1 ∈ R, and 27 < 0, which imply the solutions are real).The evolution operator is then with which can be shown from the first equation of (eqn 24) and multiplied by h i h j (summing over indices).This last relation (eqn 32) was not in Bushra Shafaq and Faisal Akram's paper [1].Finally, assuming a (anti)neutrino is produced in a pure flavour state ψ αβ (0) = δ αβ and recalling P να→ν β (x) = ψ αβ (x) 2 , one therefore has the transition probability where x = L is the propagation length, as usual.This equation is of course of the same form as the vacuum case, but writing out the effective mass differences and mixing angles is saved for a later section.

B. Details and Simplifications in Vacuum
The rest of this paper departs from [1], and is entirely original work.Notice that variable quantities such as L and E only appear in the last expression (eqn 33), except for where E and N e enter into A CC at the beginning.Because of the structure of H in terms of A CC , it will turn out that most calculations can be performed with vacuum settings (A CC = 0), and small modifications added later to take into account matter effects (see the next section).Therefore, here we take a look at the details assuming a vacuum first, where all the associated quantities will be marked with a tilde for clarity H F = HF + A.
Here HF is simply HF = U MU † , and from the cyclic nature of the trace tr[ HF ] = tr[M], so that and thus the components are explicitly given by Now, ã1 and ã0 can be computed from hi and d ijk , but it is easier to use the definitions ã1 = 1 6 tr H2 and ã0 = where ∆m 2 11 2 ≡ 2∆m 2 21 ∆m 2 31 is defined for compactness.So for a vacuum, these quantities can all be substituted in to compute E n (eqn 30), X n and P να→ν β (L) (eqn 33) directly.Notice also that the eigenvalues E n only depend on the mass differences here, as one would expect.

C. Adding Matter Effects
The values calculated above must be corrected for matter effects.From H F = HF + A, the traceless matrix H can be related to the vacuum one H simply according to Corrections to Y are also easier to see in matrix notation Since only the diagonal components of H change, from a 1 = − 1 2 tr(H 2 ) one can find and using the determinant definition of a 0 , it is modified by which one can find is simply Notice that since the diagonal components of H are real, so are those of Ỹ , and therefore a 0 and a 1 and, by extension E n , are always real.X n consequently always has real diagonal components, as expected.Finally, also note that Y and H (eqn 40, 41) have higher order matter corrections in their diagonal components.This means that survival probabilities appear to be more greatly affected by constant matter density than transition probabilities are.
where Hττ = − Hee + Hµµ was used to arrive at (eqn 60).Then the eigenvalues E n are computed exactly as before (eqn 54), and all substituted into the formulae C. Electron (Anti)Neutrino Appearance Another more complex example is muon (anti)neutrino to electron (anti)neutrino transition probability.Just as before, first compute the vacuum values ã1 , ã0 , Hee and Ỹee .However, in this case since we are dealing with transition probabilities, one must additionally compute the two complex numbers Heµ =∆m With given (anti)neutrino energy and electron density, A CC is computed as in (eqn 49), then a 0 , a 1 and E n are corrected and computed respectively just as previously (eqn 51, 52, 54).These are then substituted into so that finally where the ± is positive for neutrinos and negative for antineutrinos, while IV. RESULTS

A. Speed Comparison for Example Algorithm
In order to get a sense of the speed of this algorithm, calculations of various neutrino oscillation probabilities in constant matter density were performed by the above algorithms, written in C++.The same was then done with the widely used GLoBES package [6] -also written in C++ -and the results and processing times of the two were compared.Use was made of the glbConstantDensityProbability() function, which computes the transition or survival probability between any two neutrino flavours for constant matter density, neutrino energy and baseline.It does this by diagonalising the Hamiltonian with various numerical or analytic methods, depending on which is fastest [7,9].It is a fast and reliable method, but all matrix elements must be recomputed for each change in neutrino energy E, baseline L and matter density ρ, while this paper's method must only re-perform part of the calculations.

Method
First, an initialisation step is performed, where the precomputed constants above are calculated and GLoBES is initialised.This step was not timed since it need only be performed once and so will not scale with the number of calculations.However, note that the GLoBES initialisation takes longer since the package includes many more functionalities than just constant matter neutrino oscillations.
Second, for a given flavour transition, and a range of 100 neutrino energies, 100 baselines and 100 matter densities, the oscillation probabilities were computed using three different functions: the GLoBES function, a general flavour version of this algorithm (Section II), and then the version of this algorithm tailored to the specific flavour transition (Section III).The results and total computation times (CPU time as measured by the C++ standard library std::clock() function) of these three were recorded.This process was repeated 50 times to obtain a measure of statistical uncertainty.
The ranges of neutrino energies E, baselines L and matter densities ρ were evenly spaced values in some range where the minima were always fixed (E min = 0.5 MeV, L min = 0.01 km, ρ min = 0 g/cm 3 ).Therefore, for a given set of maxima (E max , L max , ρ max ), each function was called 50 million times (100 × 100 × 100 × 50).

Results
As alluded to, the whole process was performed for various maximum values, to discern any E, L or ρ dependence on the results.Ten different values for each were used, according to so that the whole method above was carried out one thousand times.
The computed probabilities were always exactly the same (having copied any unit conversion factors from the GLoBES code), so are not shown here.However, the computation times are presented in figure 1 for three example oscillations.Dependence on the three E max , L max and ρ max parameters is shown separately.The plot showing E max dependence averages over all L max and ρ max dependence, and likewise for the other two cases.Statistical error was propagated throughout, and added quadratically to a systematic error or σ sys = 0.01 s, being the resolution limit of the CPU time measuring function.
For zero matter density, the GLoBES calculation appears slightly faster, though is still within one standard deviation from this paper's algorithms.Every other configuration shows these to be almost or around twice as fast as GLoBES.To be specific, if T GLoBES is the CPU time taken by the GLoBES function averaged over all the data above, and likewise T general and T specific are the CPU times taken by the general and specific flavour algorithms from this paper respectively, ν e appearance algorithm was compared to other calculations described in [11].Use was made of Peter Denton's code on github, where a fork was made (https://github.com/Jamicus96/Nu-Pert-Compare) to add this paper's algorithm in two separate cases: • The first step, or initialisation of quantities in vacuum, is performed separately, before any speed comparison (branch "compare_JP_precomp").
• The initialisation is included in the speed comparison (branch "compare_JP").
It was found that the former was marginally faster than the fastest exact calculation in Peter Denton's code package ("ZS") -on the order of 6% faster -while the latter slightly slower -around 17% slower.It was estimated that accounting for the initialisation time, this paper's algorithm becomes faster than the "ZS" algorithm after 3 loops (3 probability calculations for different neutrino energies).The differences here are small, and one must bear in mind that some approximate solutions described in [11], and included in the code, are significantly faster.

B. Effective Parameters
From the correspondence between the PMNS matrix and X n matrices found earlier (eqn 33), as well as that between the eigenvalue differences and the mass differences, one can find the effective parameters (including matter effects) for some relationship between (k, j) and (m, n) indices.To deduce this relationship, note that (71) must hold for the vac-uum case, so the relationship need only be shown for that simplified case.Now, the eigenvalues of the traceless matrix H (14) in the vacuum case are clearly so the vacuum eigenvalues Ẽn must be assigned to these in some order.Next, looking at the definition of E n (30), the 1 3 cos −1 (...) term is always between 0 and π 3 .This means that and E 0 > 0 always hold.Therefore, the ordering of λ 1 , λ 2 and λ 3 determines their relationship.This ordering depends on the mass ordering itself: for Normal Ordering (NO) λ 3 > 0 > λ 2 > λ 1 , and for Inverted Ordering (IO) λ 2 > λ 1 > 0 > λ 3 .Therefore, the relationship between the indices in (71) is shown in table 1.For example in normal ordering, Notice that any dependence on energy or electron density in these parameters comes only from factors of A CC .Their values can thus be plotted on a simple graph against A CC ∝ EN e , without having to vary E and N e independently.Additionally, negative values of A CC can be used to plot the behaviour of antineutrinos, since flipping the sign of A CC is effectively the only difference.See figure 2a for the fractional scaling of these parameters in the MeV neutrino energy scale (for lithospheric electron density).To see large scale absolute changes, such as mass differences changing ordering, one must go to the GeV scale, as shown in figure 2b.As one might expect, "regime changes" (sudden changes in evolution of effective parameters) appear when A CC reaches the same scale as the mass differences.Notice also that ∆m 2 21 and s 12 2 are the most strongly affected by matter in the MeV scale.
which follows the full matter effect oscillation formula more closely.

V. CONCLUSIONS
Clearly, the derived algorithm for constant matter density transition and survival probabilities is relatively simple and efficient compared with both numerical and previous analytic solutions [4][5][6][7][8][9].
First the fact that this is an exact solution allows one to draw useful information such as anti electron neutrino survival probability being independent of the CP-phase for constant matter densities, and the effective parameters being only dependent on A CC (E, N e ) -not independently on E and N e .
Second, the performance of the example algorithms is noteworthy, computing roughly two times faster than the GLoBES package for non-zero matter density.Meanwhile, GLoBES is roughly as fast for zero matter density, being less than one standard deviation away.Ignoring the pre-computation of terms independent of energy and baseline, this paper's algorithm is also 6% faster than the fast "ZS" algorithm described above.Including this pre-computing as an initialisation, it becomes faster than "ZS" after 3 calculations.The algorithm presented here is thus optimal for applications involving a large number of oscillation probability calculations at different energies or baselines, such as oscillation analyses.
At the very least, this presents a relatively easy to implement and fast tool to compute oscillation probabilities for approximately constant matter density profiles, which scales up well for large numbers of calculations.It does not require the use and implementation of an entire package such as GLoBES, when only a simple oscillation probability is desired.

FIG. 1 :
FIG. 1: Comparison of the CPU time taken by new neutrino oscillation algorithms against the same GLoBES calculation.

TABLE 1 :
Index correspondence for effective parameters from equation 71, for Normal Ordering (NO) and Inverted Ordering (IO).