A Time-evolution operator S(T) in Our Hamiltonian Formulation

# Analytical Theory of Neutrino Oscillations in Matter with CP violation

## Abstract

We develop an exact analytical formulation of neutrino oscillations in matter within the framework of the Standard Neutrino Model assuming 3 Dirac Neutrinos. Our Hamiltonian formulation, which includes CP violation, leads to expressions for the partial oscillation probabilities that are linear combinations of spherical Bessel functions in the eigenvalue differences. The coefficients of these Bessel functions are polynomials in the neutrino CKM matrix elements, the neutrino mass differences squared, the strength of the neutrino interaction with matter, and the neutrino mass eigenvalues in matter. We give exact closed-form expressions for all partial oscillation probabilities in terms of these basic quantities. Adopting the Standard Neutrino Model, we then examine how the exact expressions for the partial oscillation probabilities might simplify by expanding in one of the small parameters and of this model. We show explicitly that for small and there are branch points in the analytic structure of the eigenvalues that lead to singular behavior of expansions near the solar and atmospheric resonances. We present numerical calculations that indicate how to use the small-parameter expansions in practice.

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

Keywords:

## I Introduction

In this paper, we develop an exact analytical representation of neutrino oscillations (1) in matter within the framework of the Standard Neutrino Model (SNM) (2) with 3 Dirac Neutrinos. The exact closed-form expressions we give for the time-evolution operator are obtained from using the Lagrange interpolation formula given in Ref. (3). The resulting expressions are easily evaluated without any approximations.

The paper is divided into two main parts. In the first, we summarize the main results of our theory. Details underlying the derivation are given in Appendices. We also retrieve the well-known two-neutrino flavor results as a special case of our general results.

In the second part we address other analytical formulations found in the literature. The expansion of the neutrino oscillation probability in one of the small parameters and of the standard neutrino model (SNM) for is of particular interest. The seminal work along these lines is found in Refs. (4); (5); (6). This work underlies many of the analyses and exploratory studies of experiments at present and future neutrino facilities, including our earlier work (7); (8); (9); (10).

The present paper was undertaken, and used, for the purpose of independently confirming the results of Refs. (7); (8); (9); (10). We find that the accuracy of the expanded oscillation probability is restricted by the presence of branch points in the analytic structure of the eigenvalues of neutrinos propagating in matter. We also show that the regions where the expanded results are reliable is different for expansions in  (4) and  (5); (6). We then map out regions where the expanded results are reliable by comparing numerical results to the exact results of our Hamiltonian formulation.

Another recent study (11) takes a complementary approach and finds that the predictions of Refs. (7); (8); (9); (10) can be improved in certain regions using an exact evaluation of the integral rather than the approximate one found there. It concludes that within these regions, predictions of oscillations improve for certain values of the experimental parameters.

The dimensionalities of the neutrino Hamiltonian and the parameter space characterizing the mixing of three neutrino pairs are sources of difficulty for finding a tractable representation of the oscillation probability. The Lagrange interpolation formula (3) is enormously helpful, providing an exact and formally elegant expression for the exponentiation of an matrix.

The description of two-flavor neutrino oscillations is elementary by comparison. In that case, is a matrix, and the mixing is described by a single real parameter.

## Ii Neutrino Dynamics

We will be interested in the dynamics of the three known neutrinos and their corresponding anti-neutrinos in matter. This dynamics is determined by the time-dependent Schrdinger 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.

The solutions of Eq. (1) may be expressed in terms of stationary-state solutions of the eigenvalue (EV) equation

 Hν|νmi> = Ei|νmi> , (3)

where the label “” indicates neutrino mass eigenstates, as distinguished from their flavor states sometimes denoted the label “”. In operator form, this dynamics may be expressed in terms of the time-evolution operator , which describes completely the evolution of states from time to and also satisfies the time-dependent Schrdinger equation.

We will examine neutrinos propagating in a uniform medium for interactions constant not only in space but also time. Because the Hamiltonian is then translationally invariant, attention may be restricted to states, both in the vacuum and in matter, characterized by momentum and therefore having the overall -dependence . In this case expressions may be simplified by suppressing the overall plane wave, a convention we adopt.

For time-independent interactions, ,

 S(t′,t) = e−iHν(t′−t) , (4)

depends on time only through the time difference . Then, written in terms of the stationary state solutions of Eq. (1),

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

With the momentum dependence factored out, three basis states are then required to describe three neutrinos. The basis states correspond to a specific representation, as in descriptions of a spin-1 object. The basis should, of course, be orthogonal,

 = δij (6)

and complete,

 ∑k|M(k)>

Once the basis is chosen, wavefunctions for a neutrino state are naturally introduced as the components of this state in the chosen basis. For example, with the eigenstates of Eq. (3) expanded in the basis,

 |νmj> = ∑i|M(i)>mij , (8)

the wave functions of would be the set . With the plane wave factored out, the wave function is just a set of three numbers. Additionally, introduction of a basis makes it possible to represent neutrino states and operators such as in matrix form, with each entry in the matrix corresponding to a projection of the object being described onto the basis.

In this paper we take the Hamiltonian in Eq. (2) to be expressed in the standard representation, where the mass basis states are taken as the set of states that diagonalize the neutrino vacuum Hamiltonian , i.e. ,

 H0v|ν0mi>=E0i|ν0mi> . (9)

In matrix form

 H0v = ⎛⎜ ⎜⎝E03000E02000E01⎞⎟ ⎟⎠ , (10)

with the EV’s taken to be ordered as in the normal mass hierarchy. In the literature, the Hamiltonian is often expressed in a different basis obtained by rotating to one in which the complete neutrino Hamiltonian is diagonal as in Ref. (4).

We assume here that that neutrinos and anti-neutrinos represented by and , respectively, are the structureless elementary Dirac fields of the the Standard Neutrino Model (2). For this reason the theory is invariant under CPT, so the mass of an anti-neutrino in the vacuum is the same as that for its corresponding neutrino.

### ii.1 Flavor and Mass States

Neutrinos are produced and detected in states of good flavor, . The three flavors, electron (), muon (), and tau () correspond, respectively, to the index values . In the vacuum, each flavor state is a specific linear combination of the three mass eigenstates of the neutrino vacuum Hamiltonian . This linear combination is expressed in terms of the same set of numbers for both neutrinos and anti neutrinos

 |ν0fi> = ∑jU∗ij|M(j)> |¯ν0fi> = ∑jUij|M(j)> , (11)

where are the elements of a unitary operator , the neutrino analog of the familiar CKM matrix. It is standard to express in terms of three mixing angles and a phase characterizing violation,

 ⎛⎜⎝c12c13s12c13s13e−iδcpU21U22s23c13U31U32c23c13⎞⎟⎠ , (12)

where

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

We use here the standard abbreviation , , etc. The parameters and are determined from experiment.

Because with it follows that the relationship in Eq. (II.1) between flavor and mass states for anti-neutrinos and neutrinos in the vacuum is equivalent to .

### ii.2 Neutrino Interacting Hamiltonian

The interaction , determined by taking the electron flavor states scattering from the electrons of the medium to mediate the interaction, is then expressed as an operator in the standard representation,

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

with and the electron number density in matter.

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

 ne = np (15) = RN ,

where is the average total nucleon number density and is the 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′ . (16)

Matrix elements of are thus

 = U∗1kVU1k′ . (17)

### ii.3 Dimensionless variables

The results are most naturally expressed in dimensionless variables. We first take advantage of the global phase invariance to express all energies relative to the vacuum EV of the same momentum. We indicate that a quantity is expressed relative to by placing a bar over it, e.g.,

 ¯E0i ≡ E0i−E01 . (18)

We follow the same convention for the Hamiltonian,

 ¯Hν ≡ Hν−1E01 , (19)

so the EV equation Eq. (3) becomes

 (¯H0v+H1)|νmi> = ¯Ei|νmi>, (20)

where

 ¯H0v ≡ H0v−1E01 . (21)

Then, to express the theory in dimensionless variables we divide all energies, including the Hamiltonian, by . The stationary-states are also be determined from the dimensionless Hamiltonian ,

 ^¯Hν = ^¯H0v+^H1 , (22)

i.e., from the solutions of

 ^¯Hν|νmi> = ^¯Ei|νmi> , (23)

where the “hat” placed over a quantity indicates it is dimensionless. Thus

 ^¯Ei ≡ ¯Ei¯E03 ^¯H0v ≡ ¯H0v¯E03 , (24)

and is obtained from by replacing

 V→^A ≡ V¯E03 . (25)

The quantity is the same as that defined in Refs. (4); (7); (8); (9); (10); (11). The connection of the Hamiltonian to the full Hamiltonian is then

 Hν = 1E01+¯E03^¯Hν . (26)

### ii.4 Neutrino Vacuum Hamiltonian ^¯H0v

The case of main interest for many situations is the ultra-relativistic limit, (we take the speed of light ). For ultra-relativistic neutrinos in the laboratory frame, the energy of a neutrino in the vacuum becomes

 E0i ≈ |→p|+m2i2E , (27)

where is its mass the vacuum. Similarly, appearing in Eq. (3) may be written

 Ei ≈ |→p|+M2i2E , (28)

where is its mass in the medium. Thus, in this limit,

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

and

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

with

 α≡m22−m21m23−m21 . (31)

In this limit, the distance from the source to the detector corresponding to in Eq. (4) is

 L = t′−t . (32)

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

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

where is given in Eq. (26), and where

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

[The similar quantity as defined in Ref. (7) is exactly one-half of that appearing in Eq. (34).]

### ii.5 Neutrino Mass Eigenvalues

The neutrino mass eigenstates in a medium are solutions to the EV equation for , Eq. (23). In many familiar formulations (4); (5); (6) the full solutions, including both the eigenstates and EV’s , are required to find the neutrino oscillation probabilities.

#### Diagonalization of Neutrino Hamiltonian

The energies are solutions of the cubic equation (12)

 ^¯Ei3+a^¯Ei2+b^¯Ei+c=0 , (35)

where

 a = −(1+α+^A) b = α+^Acos2θ13+^AαC(+)2 c = −^Aαcos2θ12cos2θ13 . (36)

We have expressed in terms of a frequently recurring combination of mixing angles,

 C(±)2 ≡ cos2θ12±sin2θ12sin2θ13 . (37)

Note that the mass eigenstate energies are independent of and for both neutrinos and antineutrinos.

The solutions of Eq. (35) are expressed conveniently in terms of the quantity ,

 d = ψ+√ψ2−4γ3 γ ≡ a2−3b ψ ≡ a3−27c−3aγ . (38)

These solutions are real when

 |d1/3|2 = 22/3γ >0 , (39)

which requires

 ψ2<4γ3 , (40)

and, thus, that be complex. Because having real energies is required by Hermiticity of the neutrino Hamiltonian, Eqs. (39,40) amount to conditions on all parameter sets in terms of which is defined.

We find

 ^¯E1 = −a3−13⋅21/3(√3Im[d1/3]+Re[d1/3]) ^¯E2 = −a3+13⋅21/3(√3Im[d1/3]−Re[d1/3]) ^¯E3 = −a3+22/33Re[d1/3] . (41)

The masses are ordered so that . Because EV do not cross, for all . A simple constraint among is found from the trace of Eq. (22),

 Tr^¯Hν = ^¯E1+^¯E2+^¯E3 (42) = Tr^¯H0v+Tr^H1 = 1+α+^A≡−a .

#### Using neutrino mass eigenvalues in our Hamiltonian Formulation

In our formulation, the entire dependence of the time evolution operator on the neutrino eigenvalues occurs through three eigenvalue combinations,

 Δ^¯Eℓℓ′ ≡ ^¯Eℓ−^¯Eℓ′ Σℓℓ′ ≡ ^¯Eℓ+^¯Eℓ′ Πℓℓ′ ≡ ^¯Eℓ^¯Eℓ′ , (43)

with (and powers thereof). We denote such quantities using a bracket notation, For example,

 Δ^¯E[1] = ^¯E3−^¯E2 Δ^¯E[2] = ^¯E3−^¯E1 Δ^¯E[3] = ^¯E2−^¯E1 , (44)

in the case of . We will generally use this bracket notation also for other quantities in our formulation that depend on two indices , such as and .

An expression for ,

 Σ[ℓ] = −a−^¯Eℓ , (45)

follows from Eq. (42). An equivalent expression for in terms of is found by subtracting Eq. (35) for and that for and dividing through by . We find

 0 = (^¯Eℓ2+^¯Eℓ^¯Eℓ′+^¯Eℓ′2) (46) + a(^¯Eℓ+^¯Eℓ′)+b = (^¯Eℓ+^¯Eℓ′)2−^¯Eℓ^¯Eℓ′ + a(^¯Eℓ+^¯Eℓ′)+b ,

giving

 Σ[ℓ]2−Π[ℓ]+aΣ[ℓ]+b = 0 . (47)

Then, using Eq. (45),

 Π[ℓ] = Σℓ(Σℓ+a)+b (48) = b+a^¯Eℓ+^¯Eℓ2 .

Finally, having observed that powers of the quantities given in Eq. (II.5.2) will appear in various expressions, we note that and with and involve linear combinations of eigenvalues with powers . Such terms are equivalently represented by a linear combination of three terms, one proportional to , one proportional to , and one independent of , obtained by using the equation of motion repeatedly. We later make use of this fact to simplify various expressions.

## Iii The S-Matrix in Our Hamiltonian Formulation

The probability for neutrinos to oscillate from the initial state of flavor to a final state of flavor is found from the time-evolution operator as

 P(νa→νb) = |Sab(t′,t)|2 (49) ≡ Pab(t′−t) ,

where

 |Sab(t′,t)|2 = |<ν0fb|S(t′,t)|ν0fa>|2 . (50)

We accordingly determine here from defined in Eq. (4).

In this section we review the formulation of neutrino oscillations based on the Lagrange interpolation formula as used in Ref. (3). This formulation leads to exact, closed-form expressions for the time-evolution operator and the partial oscillation probabilities that are linear combinations of spherical Bessel functions in the eigenvalue differences whose coefficients are polynomials in the neutrino CKM matrix elements, the neutrino mass differences squared, the strength of the neutrino interaction with matter, and the neutrino mass eigenvalues in matter. We are led quite naturally to such expressions for all the partial oscillation probabilities in terms of these basic quantities. The numerical results given later in this paper are based on this formulation.

### iii.1 Time-Evolution Operator

The overall phase in Eq. (33) does not contribute to , so for the purpose of calculating the oscillation probability, we may take

 S(L) → e−i^¯HνΔL . (51)

Then, with neutrinos created and detected in flavor states, which are coherent linear combinations of the neutrino vacuum mass eigenstates given in Eq. (9),

 |ν0fa> = ∑jU∗aj|M(j)> , (52)

we see that the mass eigenstate components of the flavor states contribute coherently to the time-evolution operator. Thus,

 <ν0fb|e−i^¯HνΔL|ν0fa> (53) =  .

This coherence leads to the oscillation phenomenon.

The elegant formulae for are obtained from the Lagrange interpolation formula, Eqs. (9,11) of Ref. (3),

 Ue−i^¯HνΔLU−1 = ∑ℓFℓexp−i^¯EℓΔL , (54)

where and

 Fℓ ≡ Πj≠ℓU^¯HνU−1−1^¯Ej^¯Eℓ−^¯Ej . (55)

For three neutrinos, the sum in Eq. (54) runs over three values of and the product in Eq. (55) over two values of .

Using the convention that , written without parentheses enclosing , denotes the matrix elements of the operator ,

 Oab ≡  , (56)

the matrix elements of given in Eq. (55) may be compactly written

 Fabℓ = ^¯D[ℓ] , (57)

where (3),

 ^¯W[1] ≡ (U^¯HνU−1−1^¯E3)(U^¯HνU−1−1^¯E2) ^¯W[2] ≡ (U^¯HνU−1−1^¯E3)(U^¯HνU†−1^¯E1) ^¯W[3] ≡ (U^¯HνU−1−1^¯E2) (58) × (U^¯HνU†−1^¯E1)

and

 ^¯D[1] = (^¯E3−^¯E1)(^¯E2−^¯E1) ^¯D[2] = (^¯E1−^¯E2)(^¯E3−^¯E2) ^¯D[3] = (^¯E1−^¯E3)(^¯E2−^¯E3) . (59)

Equations (III.1,III.1) use the same bracket notation introduced in Eq. (34). The result in Eqs. (53,54,55) is immediately verified to be correct by inserting a complete set of intermediate eigenstates of in Eq. (III.1).

It follows from the unitarity of that is Hermitian,

 ^¯W[ℓ]† = ^¯W[ℓ] (60)

and that the two factors in Eqs. (III.1) commute with each other. We find from Eq. (III.1) that

 Tr^¯W[ℓ] = ^¯D[ℓ] . (61)

Equation (60) establishes the reflection symmetry

 Fab∗ℓ = Fbaℓ . (62)

Explicit expressions for are easily found in terms of . The entire dependence of on occurs through three operators independent of , , and ,

 ^¯W[ℓ] = ^¯W0[ℓ]+cosδcp^¯Wcos[ℓ] (63) + isinδcp^¯Wsin[ℓ] ,

with , and real and independent of . Details are given in Appendix A.

### iii.2 Total Oscillation Probability

Expressions for may be obtained directly from ,

 P(νa→νb) = |Sab(t′,t)|2 (64) = Re[Sab(L)]2+Im[Sab(L)]2 .

Convenient expressions for and defined by Eq. (54) are presented in Appendix A. In our Hamiltonian formulation, the dependence of on the CP violating phase is very simple and follows from Eqs. (54,63) noting that , Eq. (57). We thus find,

 Re[Sab(t′,t)] = δab−2∑ℓ^¯Wab0[ℓ]^¯D[ℓ]sin2^¯EℓΔL (65) − 2cosδcp∑ℓ^¯Wabcos[ℓ]^¯D[ℓ]sin2^¯EℓΔL + sinδcp∑ℓ^¯Wabsin[ℓ]^¯D[ℓ]sin2^¯EℓΔL ,

and

 Im[Sab(t′,t)] = −2sinδcp∑ℓ^¯Wabsin[ℓ]^¯D[ℓ]sin2^¯EℓΔL (66) − ∑ℓ^¯Wab0[ℓ]^¯D[ℓ]sin2^¯EℓΔL − cosδcp∑ℓ^¯Wabcos[ℓ]^¯D[ℓ]sin2^¯EℓΔL ,

where is defined in Eq. (34). Approximate expressions for in terms of the parameters of were obtained from in Refs. (5); (6) by an expansion in .

### iii.3 Partial Oscillation Probabilities

Using somewhat different techniques, the oscillation probability may be expressed through a set of functions that express how depends on the CP violating phase  (4). In our Hamiltonian formulation there are four such terms,

 Pab = δ(a,b)+Pab0+Pabsinδ+Pabcosδ (67) + Pabcos2δ ,

with linear in , linear in , quadratic in , and independent of . Although only the overall oscillation probability is a true probability, guaranteed to be strictly positive everywhere, we find it convenient to refer to these four terms as “partial oscillation probabilities.” Approximate expressions for the partial oscillation probabilities expanded in the small parameter of the SNM in Ref. (4).

Obtaining expressions for the partial oscillation probabilities from the time-evolution operator requires additional analysis, given in Appendix B. In terms of spherical Bessel functions, we find there,

 Pabsinδ(ΔL,^A) = sinδcp4ΔL^D∑ℓ(−1)ℓwabsin[ℓ] (68) × j0(2^Δ[ℓ]) ,

where , and therefore , is anti-symmetric under . The other three partial oscillation probabilities are individually symmetric under . We find

 Pabcosδ(ΔL,^A) = −cosδcp4Δ2L^¯D∑ℓ(−1)ℓwabcos[ℓ] × Δ^¯E[ℓ]j20(^Δ[ℓ]) Pabcos2δ(ΔL,^A) = −cos2δcp4Δ2