Simple and Accurate Oscillation Probabilities for Three Coupled Neutrinos Propagating in Matter

# Simple and Accurate Oscillation Probabilities for Three Coupled Neutrinos Propagating in Matter

Mikkel B. Johnson
Los Alamos National Laboratory, Los Alamos, NM 87545
Leonard S. Kisslinger
Department of Physics, Carnegie-Mellon University, Pittsburgh, PA 15213
###### Abstract

Within a conventional Hamiltonian description, we find accurate closed-form expressions for the oscillation probabilities of three coupled neutrinos propagating in matter. Subtle cancelations that occur in coefficients of our formulation are avoided for all transitions by transforming to a different set of coefficients presented in an appendix of this paper. The neutrino mass eigenvalues are easily obtained numerically as the solution of a cubic equation. Our methods are illustrated for flavor-changing transitions in the sector. The resulting analytic expressions oscillation probabilities, which are particularly simple, are also accurate to a few percent over all regions of interest at present and the envisioned future neutrino facilities. While somewhat less accurate than numerical simulations, our approximate expressions are sufficiently accuracy to obviate the need for exact computer simulations in many circumstances.

PACS Indices: 11.3.Er,14.60.Lm,13.15.+g

Keywords:

## I Introduction

In this paper, we develop methods leading to approximate analytical expressions for the oscillation probability from the exact, closed-form results in Ref. jhk0 (). These results are based on a conventional Hamiltonian and the Standard Neutrino Model ISS () and lead to expressions for the oscillation probability for all transitions .

We include the interaction of neutrinos with electrons shown to be essential by Wolfenstein in 1978 wolf () and later identified by Mikheyev and Smirnov SMAS () as a likely explanation of the deficit of solar neutrinos discovered experimentally by Davis Davis (); Clev ().

Subsequent to the work of Wolfenstein, neutrino oscillations in matter have been commonly explored using exact computer simulations. These methods have advantages for density profiles having a spatial dependence, are discussed in Ref. MS () and references therein.

Comparisons of Ref. jhk0 () to currently available analytical results f (); ahlo (); JO () are found in Ref. jhk1 (), where it is established that these results f (); ahlo (); JO () are reliable over rather narrow regions of energy, density, and baseline. Expressions found here do not suffer from these limitations.

In the present paper we use the published formulation of Ref. jhk0 (). Because we rely heavily on all results found in Ref. jhk1 (), the reader will find many of the important results of the unpublished Ref jhk1 () repeated here.

Our formulation jhk0 () is introduced in Sect. II with details relegated to appendices. The oscillation probability simplifies in part because many of the terms appearing in the exact coefficients of our approach jhk0 () are quite small [, ] and can therefore be eliminated.

The expressions we find for the oscillation probabilities are built on coefficients that are devoid of subtle cancelations occurring in the coefficients of Ref. jhk0 (). The nature of these cancelations is discussed in Sect. III, and their source is identified in Sect. III.2. The cancellations are avoided by making a transformation that leaves the partial oscillation probabilities invariant. The resulting coefficients are essential for obtaining simple expressions for the observable oscillation probabilities accurate to [about 1%]. The coefficients are tabulated in Appendix C of the present paper for all transitions .

Observable oscillation probabilities accurate to are easily found using and eigenvalues obtained by numerically by solving a cubic equation. However, the emphasis of this paper is on developing methods for obtaining simple and accurate expressions for the oscillation probabilities. These methods are based on the same coefficients .

Application of these methods is illustrated for flavor-changing transitions in the sector. The oscillation probabilities so obtained are accurate to a few percent over all regions of interest at present and envisioned future neutrino facilities. While somewhat less accurate than numerical simulations, our approximate expressions are sufficiently accurate to obviate the need for exact computer simulations in many circumstances.

Our results rely on neutrino mass eigenvalues evaluated using first-order perturbation theory in one of the small parameters of the SNM. Within the solar resonance region, the simplest expressions accurate to are based on the -expanded eigenvalues. Above the solar resonance region they are based on the -expanded eigenvalues.

The procedure for obtaining the expressions for the oscillation probabilities for all transitions is identified in Sect IV.1.1. As before, we illustrate the procedure for the specific case of transitions. Expressions for our simplified oscillation probabilities are given in Appendix F.

As emphasized in Ref. zx (), accurate and reliable expressions for oscillation probabilities are essential for analysis and prediction at future neutrino facilities. We calibrate these aspects of our results over all regions of energy and matter density of interest. In Sect. V, we present comparisons within the solar resonance region, and in Sect. VI we present comparisons above the solar resonance region.

## Ii Neutrino Dynamics

The dynamics of the three known neutrinos and their corresponding anti-neutrinos in matter is determined by the time-dependent Schroedinger equation,

 iddt|ν(t)> = Hν|ν(t)> , (1)

where the neutrino Hamiltonian ,

 Hν=H0v+H1 , (2)

consists of a piece describing neutrinos in the vacuum and a piece describing their interaction with matter.

Our formulation jhk0 () applies to neutrinos propagating in a uniform medium for interactions constant not only in space but also time, and it assumes that that neutrinos and anti-neutrinos represented are the structureless elementary Dirac fields of the Standard Neutrino Model ISS ().

### ii.1 The Neutrino Hamiltonian in Matter

Neutrinos are characterized by flavor and mass . These states are related by the neutrino analog of the familiar CKM matrix,

 νf = Uνm . (3)

They are produced and detected in states of a specific flavor.

The unitary matrix is often parametrized in terms of three mixing angles and a phase characterizing violation.

Using the standard abbreviations, , , etc,

 U ≡ ⎛⎜⎝c12c13s12c13s13e−iδcpU21U22s23c13U31U32c23c13⎞⎟⎠ , (4)

where,

 U21 = −s12c23−c12s23s13eiδcp U22 = c12c23−s12s23s13eiδcp U31 = s12s23−c12c23s13eiδcp U32 = −c12s23−s12c23s13eiδcp . (5)

In this representation, the index of corresponds to the electron neutrino , the index to the muon neutrino , and the index to the tau neutrino .

The perturbing Hamiltonian is determined by the interaction between electron neutrino flavor states and the electrons of the medium. Expressed as matrix,

 H1 = U−1⎛⎜⎝V00000000⎞⎟⎠U , (6)

With and the electron number density in matter.

### ii.2 The Standard Neutrino Model

For electrically neutral matter consisting of protons, neutrons, and electrons, the electron density is the same as the proton density ,

 ne = np (7) = RA ,

where is the average total nucleon number density of matter through which the neutrinos propagate and is its average proton-nucleon ratio. In the earth’s mantle, the dominant constituents of matter are the light elements so ; in the surface of a neutron star . Matrix elements of are thus

 = U∗1kVU1k′ . (8)

Using the well-known expression for , we find the corresponding the (dimensionless) interaction strength of neutrinos and anti-neutrinos with matter to be,

 ^A = ±6.50 10−2R E[GeV]ρ[gm/cm3] , (9)

with E[GeV] being the neutrino beam energy in GeV and [gm/cm] the average total density of matter through which the neutrino beam passes on its way to the detector in . For experiments close to the earth’s surface, the appropriate density is the mean density of the earth’s mantle,

 ρ[gm/cm3] = ρ0 (10) ≈ 3 .

We adopt the Standard Neutrino Model ISS (), given next, to complete our description of the neutrino Hamiltonian. Most of the parameters of the SNM are consistent with global fits to neutrino oscillation data with relatively good precision hjk1 (); khj2 (). These include the neutrino mass differences,

 m22−m21 ≡ δm221 (11) = 7.6×10−5 eV2

and

 m23−m21 ≡ δm231 (12) = 2.4×10−3 eV2 ,

which corresponds to

 α ≡ δm221δm231 (13) = 3.17×10−2 .

The mixing angles are also determined from experiment. In the SNM, the value of ,

 θ23 = π/4 , (14)

is the best-fit value from Ref. Dav (), and ,

 θ12 = π/5.4 , (15)

is consistent with the recent analysis of Ref. Gon (). The mixing angle is known to be small ( at the 95% confidence level), but until recently its precise value has been quite uncertain. Results from the Daya Bay project DB () have measured its value more accurately, , which we adopt to determine our value for ,

 θ13 = 0.151 . (16)

This fixes . The CP violating phase is not known at all, and determining its value will one of the major interests at future neutrino facilities.

### ii.3 Our Hamiltonian Formulation

Our approach is based on the evaluation of the time-evolution operator using the Lagrange interpolation formula barg (). For time-independent interactions, is written,

 S(t′,t) = e−iHν(t′−t) . (17)

It depends on time only through the time difference and mat then be written in terms of the stationary state solutions of Eq. (1),

 S(t′,t) = ∑i|νmi>e−iEi(t′−t)<νmi| . (18)

Because the rest masses of neutrinos are considered to be tiny, for most cases of interest including our approach, the ultra-relativistic limit, (we take the speed of light ) is assumed. Ultra-relativistic neutrinos of energy in the laboratory frame may be expressed,

 E ≈ |→p|+m22E , (19)

where is its mass in the vacuum.

Thus, in this limit and in dimensionless variables,

 ^¯Ei → M2i−m21m23−m21 (20)

and

 ^¯H0v → ⎛⎜⎝0000α0001⎞⎟⎠ (21)

with

 α≡m22−m21m23−m21 . (22)

The distance from the source to the detector corresponding to in Eq. (17) is

 L = t′−t . (23)

The time-evolution operator, Eq. (17), expressed in dimensionless variables is then,

 S(L) = e−iHν(t′−t) (24) = e2i^¯E01ΔLe−2i^¯HνΔL ,

where is the full neutrino Hamiltonian expressed in dimensionless variables, and where

 ΔL ≡ L(m23−m21)4E . (25)

Our formulation is summarized in Appendix A.

Taking the value of from the SNM, in the high-energy limit [ Eq. (25)] becomes,

 ΔL ≈ 3.05×10−3L[Km]E[GeV] , (26)

with the baseline in . From Eqs. (26,9), we find

 ΔL^A = ±1.92 10−4R ρ[gm/cm3]L[Km] , (27)

or

 L[Km] = ±5.05 103RΔL^A ρ[gm/cm3] . (28)

Neutrinos at rest in a neutron star were recently considered in Ref. K (). It was recognized that the only modification required was to take the non-relativistic limit of the vacuum Hamiltonian. In our approach, this would mean taking

 E0i → mi+p22mi , (29)

rather than Eq. (19). In particular, the eigenvalues of the interacting Hamiltonian continue to be expressed as the solution of the same cubic equation [see Eqs. (34,35) of Ref jhk0 ()] but now with ,

 ^α≡m2−m1m3−m1 . (30)

Because it is quite straightforward to obtain the eigenvalues by solving the cubic equation numerically, if all one requires is an answer, the preferable approach would be to obtain the oscillation probability in our Hamiltonian formulation using the expressions given in Appendix A with the coefficients given in Appendix C. The advantage of the analytical result we pursue in the remainder of the paper is the insight into the underlying physics that simple analytical results provide.

## Iii Small Parameter Expansions of Observable Oscillation Probabilities

Our goal is to obtain the simplest possible expressions for the neutrino oscillation probabilities accurate to about 1%. This is accomplished using the small-parameter expansions of the exact analytical expressions of our formulation. The small-parameter expansion is facilitated by eliminating in favor of , writing , and then expanding in . Because in the SNM, a Taylor expansion in avoids having to expand separately in and to simplify expressions.

### iii.1 The Observable Oscillation Probabilities

The observable oscillation probabilities are expressed in terms of the partial oscillation probabilities defined in Ref. jhk0 () and summarized in Appendix A of this paper. Our methods are applicable to these oscillation probabilities for all transitions .

One of the observable oscillation probabilities, , is the sum over three of the four partial oscillation probabilities defined in Eq. (A.1) and is symmetric under ,

 P+ab(ΔL,^A) ≡ Pabcosδ(ΔL,^A)+Pabcos2δ(ΔL,^A) (31) + Pab0(ΔL,^A) .

The structure of is rather complicated, and simplifying this quantity is the focus of this paper.

The other observable oscillation probability, , is defined in Eq. (53). It is antisymmetric under . Because is both exact and already sufficiently simple, it will not be discussed further in this paper.

### iii.2 Leading Terms of the α-Expanded w(ab)i,p

A key quantity of our formulation jhk0 () is the set of coefficients given explicitly in Ref. jhk1 (). These coefficients are polynomials in , and for the purpose of the present work, it is rather important to note that for many transitions , the terms in the polynomial expressions for with the lowest powers in (the “leading terms”) cancel against each other when or . One or both of these conditions hold for at least one or two of the eigenvalues over nearly the entire range of . The consequence is that if we express the oscillation probabilities in terms of the original set of coefficients , three powers of must be retained to achieve accuracy, whereas we might have expected to retain only two. Polynomials with two powers of are, of course, simpler than those with three.

Polynomials with accuracy with two powers of do, in fact, exist. To find them, it is important to recognize that there is nothing unique about the original set of coefficients . We could have chosen any set of coefficients say, , as long as the partial oscillation probabilities remain invariant under . An example of such a transformation is,

 w(ab)i;0− ≡ w(ab)i;0+w(ab)i;1+w(ab)i;2 w(12)i;1− ≡ w(12)i;1+w(12)i;2 w(12)i;2− ≡ w(12)i;2 . (32)

Clearly, the partial oscillation probabilities are invariant under this transformation since,

 ^¯wabi[ℓ] → (w(ab)i;0−+w(ab)i;1−(^¯Eℓ−1)+w(ab)i;2− (33) × ^¯Eℓ(^¯Eℓ−1))Δ^¯E[ℓ] = (w(ab)i;0+w(ab)i;1 ^¯Eℓ+w(ab)i;2 ^¯E2ℓ) × Δ^¯E[ℓ] .

The important points to note are, first, that the terms in the polynomial expressions for , having the lowest powers in cancel only in isolated cases. Secondly, the quantity is expressed equivalently in terms of both and . Thus, the simplest expressions for the oscillation probabilities of accuracy are found from the first two powers in in the expansion of , whereas the same accuracy using would require the first three powers of . The coefficients and are, of course, no longer equivalent when they are truncated.

Next, we compare calculations of an approximate oscillation probability using the first two leading terms of with an exact oscillation probability calculated with .

Figure 1 shows a comparison of the two calculations of over the interval for a value of . In Fig. 2 shows a comparison of the two calculations of over the interval for fixed . As expected, the calculation using the first two leading terms of agrees very well with the exact one in both cases.

We have examined the percentage error in the partial oscillation probabilities and calculated using the first two leading terms of by comparing them numerically to the exact result. The exact and approximate calculations agree within a fraction of a percent, except in the vicinity of zeros. The error unavoidably diverges sufficiently close to the zeros when they are displaced by even a small amount.

We find the same cancellation just discussed for transitions occurs for most of the other transitions as well. We have tabulated the truncated coefficients for all these transitions in Appendix C.

The coefficients , used in conjunction with a numerical evaluation of the neutrino mass eigenvalues, make it possible to find simple expressions for the observable oscillation probabilities applicable to all transitions accurate to [about 1%] throughout the entire region of interest at present and envisioned future neutrino facilities. The mass eigenvalues are easily obtained numerically as the solution of a cubic equation.

## Iv Simple Analytic Expressions for P+ab(ΔL,^A)

In this section, we develop methods that lead to simple analytic expressions for the observable oscillation probabilities of all transitions . The analytical results are simple, in part, because many of the terms appearing in the exact coefficients of our approach jhk0 () are quite small [, ] and, for this reason, can be eliminated. Application of the methods is illustrated for flavor-changing transitions in the sector.

### iv.1 Simplifying P+ab(ΔL,^A)

Consider the partial oscillation probability contributing to in Eq. (31). The idea making it possible to find remarkably simple expressions for can be explained in terms of the corresponding to as follows.

Expanding we find

 ^¯wabi[ℓ] = ∑n,iCn,i(ℓ,α)^An , (34)

where, of course,

 ^¯wabi[ℓ] ≡ (w(ab)i;0+w(ab)i;1 ^¯Eℓ+w(ab)i;2 ^¯E2ℓ) (35) × Δ^¯E[ℓ] ,

and is the coefficient of in the expansion of . Then, expanding and retaining the first two terms in the power series in ,

 Cn,i(ℓ,α) ≈ αp(n)( C(ℓ)0;n,i (36) + αC(ℓ)1;n,i+ ...) .

From these expressions, we can easily see that evaluating in terms of individual accurate to about 1% fails to give the simplest expressions for when or . Said otherwise, when or , simplifying the individual simplifying does not lead to the simplest expression for .

In this situation, the approach that does work is to first “consolidate” the coefficients , which means to (1) construct the single quantities defined in Eq. (38); (2) simplify ; and, (3) construct from this simplified . The reason why consolidation works is most easily understood by considering ,

 C(ΔL,^A) ≡ −4Δ2L^¯D∑ℓ(−1)ℓ^¯wab[ℓ] (37) × j20(ΔLΔ^¯E[ℓ]) ,

where,

 ^¯wab[ℓ] ≡ cosδcp^¯w(ab)cos[ℓ]+cos2δcp^¯w(ab)cos2[ℓ] (38) + ^¯w(ab)0[ℓ] .

The quantity bears a simple relation to ,

 P+ab(ΔL,^A) = C(ΔL,^A) . (39)

The important point is that when or , the smaller coefficients are seldom among the first two leading terms of defined in Eq. (38) and, for this reason, they do not appear in the simplified . However, when the individual of Eq. (35) are simplified, these small terms are often retained and appear in when they are added together to obtain .

The quantity appearing in Eq. (37) is different for each region. In particular, within the deep solar resonance region, within the far solar resonance region, and above the solar resonance region.

In Eq. (38), the quantity is,

 ^¯w(ab)i[ℓ] = (w(ab)i;0−+w(ab)i;1−(^¯Eℓ−1)+w(ab)i;2− (40) × ^¯Eℓ(^¯Eℓ−1))Δ^¯E[ℓ] ,

where the coefficients are found in Appendix C for all transitions .

#### iv.1.1 ^Cα(^A) and CT(Rs): Independent Variables

The error in , so consolidated, is most easily seen to be of order of 1% by dropping all terms of relative size with and recognizing that [see Eq. (13)]. Said otherwise, an accuracy of about 1% is achieved by retaining just the first two terms of the -expansion of ( and ).

In this case, takes the form,

 Pabi(ΔL,^A) → −4Δ2L^¯D∑n,ℓαp(n)(C(ℓ)0;n,i+αC(ℓ)1;n,i) (41) × ^Anj20(ΔLΔ^¯E[ℓ]) .

Equation 41) is valid only when is taken to be independent of , but the argument is easily generalized to take account of the dependence of on (see Sect. IV.1.2). Similar ideas apply to expressions within the solar resonance region where .

#### iv.1.2 ^Cα(^A) and CT(Rs): Dependent Variables

Of course, and are actually independent variables. We next consider the more realistic case where depends on and depends on .

Although this discussion is somewhat technical, it explains in detail how we found the simplified expressions for given in Appendix D, on which the numerical results of Sects. V and VI documenting their accuracy is based. Because this discussion is general, it applies to all transitions . Although we assume that we are above the solar resonance region, the argument applies equally as well within the solar resonance region by replacing .

When , we first factor , Eq. (37). One of these factors is proportional to . This dependence arises from the -expanded eigenvalues appearing in Appendix B of Ref jhk1 () and associated with the term . The product of the other factors is proportional to , .

We then replace each term proportional to an even power, say , of by . Similarly, each term proportional to an odd power, say , of by .

After these replacements, a straightforward generalization of Eq. (41) gives,

 C(ΔL,^A) → −4Δ2L^¯DC3α∑n,ℓαp(n)((C(ℓ)0;n+αC(ℓ)1;n)+^Cα(C′(ℓ)0;n+αC′(ℓ)1;n))^Anj20(ΔLΔ^¯E[ℓ]) , (42)

where and are the two sets of first two leading terms of .

#### iv.1.3 Summary and Discussion

To summarize, by expressing in terms of the single coefficient in Eq. (40), the entire set of terms and are collected together with the consequence that the first two leading terms of is .

At the same time, some of the higher-order terms that would have been retained by considering and independently are naturally discarded, leading to the result of . It is the result of because the combination has fewer terms when one of the two is smaller than the other. Inspection of the coefficients given in Appendix D indicates that this is often the case.

Since the first two leading terms constitute the simplest and most accurate expressions for the oscillation probability, this method is guaranteed to lead uniquely to simplest and most accurate expression for to about 1%, which is what we set out to find. We note in passing that we have managed to cast “simplicity” into mathematical language. Some might find this interesting, since “simplicity is normally considered to be a subjective concept.

The resulting oscillation probabilities were seen to be accurate to a few percent over all regions of interest at present and envisioned future neutrino facilities. These analytic oscillation probabilities are vastly more accurate than the familiar analytic expressions on the interval and for values of extending from well into the asymptotic region, . While somewhat less accurate than numerical simulations, our approximate expressions are sufficiently accuracy to obviate the need for exact computer simulations in many circumstances. The accuracy of the approximate found in this way is sufficient for analysis and prediction of experiments at present and future neutrino facilities. Furthermore, the accuracy is easy improved applying a well-defined correction procedure.

Note the following caveats. Depending on the values of the parameters of the SNM, for specific regions it may happen that (1) some pieces of the first two leading terms are small enough to drop to maintain a specific accuracy goal; (2) the terms of relative order may happen to be anomaly large and thus retained to maintain an accuracy goal.

These results may be used to define “effective” partial oscillation probabilities. The “effective” partial oscillation probability would represent the dependence of on . Similarly, an “effective” partial oscillation probability would represent the dependence of on , and an “effective” would represent terms independent of . These effective partial oscillation determine, in turn, effective . Note, however, that these effective partial oscillation probabilities are not guaranteed to be the simplest to .

## V Simplified Partial Oscillation Probabilities, ^A<0.1

In this section, we consider the simplified oscillation probability over the interval , where, in this paper, . This interval encompasses the solar resonance, which is located at . For this reason, the interval is referred to in this paper as the solar resonance region. Note that this definition of the solar resonance region differs from that adopted in Ref. jhk0 (), where .

Here, we split up the solar resonance region into two sub-regions. The sub-region is referred to as the deep solar region, and the sub-region as the far solar region. For these sub-regions, we find it sufficient to retain only the first leading term of eigenvalue difference, , , where is the first term of the Taylor series expansion of in appearing in Appendix D.2.

The simplified within the solar resonance region is given explicitly in Appendix F by applying the procedure of Sect. IV. We examine this numerically in the present section. By comparing to its exact counterpart, we will establish that describes the - and -dependence of the exact oscillation probability within the solar resonance region much more reliably than those presently found in the literature.

### v.1 General Considerations

From the simplified expressions for given in Appendix F, where our results are expressed our in terms of effective partial oscillation probabilities and effective , we will see that the effective is generally non-vanishing. For this reason, is sensitive to the CP violating phase, throughout this region. We will also see there that the situation is different above the solar resonance region.

Within the solar resonance region, the expansion in is the appropriate one to use to obtain the eigenvalues. We use the representation of the -expanded eigenvalues given in Appendix of Ref. jhk1 (). This representation motivates splitting the solar resonance region into the deep solar resonance region, and the far solar resonance region, .

The advantage of this representation is that by making the replacement,

 ^A → αRs ^Cθ(^A) → αCT(Rs) , (43)

in the eigenvalues of the deep solar resonance region, and,

 ^A → Rs/α ^Cθ(^A) → ^ACT(Rs) , (44)

in those of the far solar resonance region, the function of ,

 CT(Rs) ≡ √1+R2s−2Rscos2θ12 , (45)

appears in both sets of eigenvalues.

Expressions for the effective partial oscillation probabilities given in Appendix F have been simplified using the properties that for all , is positive and less than or equal to one and that

 CT(Rs