Describing neutrino oscillations in matter with Magnus expansion

# Describing neutrino oscillations in matter with Magnus expansion

A. N. Ioannisian, A. Yu. Smirnov Yerevan Physics Institute, Alikhanian Br. 2, 375036 Yerevan, Armenia
Institute for Theoretical Physics and Modeling, 375036 Yerevan, Armenia
International Centre for Theoretical Physics, Strada Costiera 11, 34014 Trieste, Italy
Institute for Nuclear Research, Russian Academy of Sciences, Moscow, Russia
###### Abstract

We present new formalism for description of the neutrino oscillations in matter with varying density. The formalism is based on the Magnus expansion and has a virtue that the unitarity of the S-matrix is maintained in each order of perturbation theory. We show that the Magnus expansion provides better convergence of series: the restoration of unitarity leads to smaller deviations from the exact results especially in the regions of large transition probabilities. Various expansions are obtained depending on a basis of neutrino states and a way one splits the Hamiltonian into the self-commuting and non-commuting parts. In particular, we develop the Magnus expansion for the adiabatic perturbation theory which gives the best approximation. We apply the formalism to the neutrino oscillations in matter of the Earth and show that for the solar oscillation parameters the second order Magnus adiabatic expansion has better than accuracy for all energies and trajectories. For the atmospheric and small 1-3 mixing the approximation works well ( accuracy for ) outside the resonance region (2.7 - 8) GeV.

###### pacs:
14.60.Pq, 95.85.Ry, 14.60.Lm, 26.65.+t

## I Introduction

Neutrino physics enters the era of precision measurements, studies of the sub-leading oscillation effects and searches for new physics beyond the standard neutrino scenario. The neutrino flavor conversions become a tool of exploration of other particles and objects such as interiors of the Earth and stars. One of the key elements of these studies is neutrino oscillations in matter with varying density, and in particular, the oscillations inside the Earth. The latter is relevant for the solar, supernova and atmospheric neutrinos, as well as for the cosmic and accelerator neutrinos. In this connection it is important to have precise analytical or semi-analytical expressions for oscillation probabilities valid in wide energy ranges. These expressions allow us to simplify numerical computations but also to gain a deeper insight into physics involved. The results can be of special interest in view of discussions of future experiments with the megaton-scale fine structured underwater/underice detectors.

Several analytic and semi-analytic approaches to computing probabilities in matter with non-constant density have been developed recently which use various perturbation theories deHolanda:2004fd (), Ioannisian:2004jk (), Ioannisian:2004vv (), Akhmedov:2004rq (), Blennow:2004qd (), Akhmedov:2005yj (), LM (), Akhmedov:2006hb (), deAquino:2007sx (), Supanitsky:2007ks (), Liao:2007re (). In the previous publications Ioannisian:2004jk (), Ioannisian:2004vv (), we have proposed a formalism which describes the neutrino oscillations in matter with low density. It make use of smallness of the matter potential in comparison with the kinetic term: , where is the mass squared difference and is the energy of neutrino. Essentially, the expansion parameter is given by the integral along the trajectory

 I=∫dx V(x)cosϕ(x),

where is the adiabatic phase. The first approximation works very well at low energies MeV Ioannisian:2004jk (). Validity of the results can be extended to higher energies if the second order term, , is taken into account Ioannisian:2004vv (). It can be further improved in certain energy ranges if expansion is performed with respect to the deviation of the potential from some average value.

The problem of this and some other similar approaches is that the unitarity of oscillation amplitudes is not guaranteed, and in fact, is violated at high energies Ioannisian:2004vv (). This violation, in turn, can produce certain problems in numerical computations. In this paper we propose the new type of perturbation theories which maintain the unitarity explicitly in each order of expansion, and therefore at any truncation of the series. The approach is based on the Magnus expansion magnus (), burum () which was previously used for description of the nonadiabatic neutrino conversion in medium with monotonously varying density D'Olivo:1990xs (), D'Olivo:1992mw () AguilarArevalo:2003an (). Recently the first order Magnus expansion has been applied to the low energy neutrino oscillations in matter of the Earth Supanitsky:2007ks (). The formula for the regeneration factor in the Earth has been obtained which generalizes our result in Ioannisian:2004vv (). In this paper we develop various perturbation theories using explicitly two orders of the Magnus expansion. Since the Magnus expansion is an expansion in power of commutators, it is the second order that provides non-trivial new results. As a part of the present study we reproduce the formula from Supanitsky:2007ks ().

Essentially, the restoration of unitarity in the Magnus expansion is achieved by an effective re-summation of certain contributions to oscillation amplitudes. This leads to higher accuracy of the semi-analytic results and allows us to further extend the range of applications of the approach. Furthermore, it gives better understanding of the previously obtained results and their limits of validity.

We illustrate an accuracy of the approximations computing the transition probabilities for neutrinos crossing the core of the Earth. We find that for the solar oscillation parameters the second order Magnus adiabatic expansion has better than accuracy for all energies and all trajectories. For the atmospheric and small 1-3 mixing the approximation works very well ( accuracy for ) below 2.7 GeV and above 8 GeV for . In the region, (2.7 - 8) GeV, where the MSW resonances in the core and in the mantle as well as the parametric resonances take place, a special consideration is required.

The paper is organized as follows. In sec. 2 we present the formalism of Magnus expansion and obtain general expressions for the S-matrix. We calculate the oscillation probabilities using various perturbation approaches based on the Magnus expansion in sec. 3. In particular, we develop the perturbation theory in and the adiabatic perturbation theory. We compare the results of different semi-analytic approaches in sec. 4. Conclusions follow in sec. 5.

## Ii Magnus expansion

### ii.1 S-matrix and Magnus expansion

In what follows we will mainly study the case of mixing , where is, in general, some combination of and . In a number practical cases the two neutrino results can be immediately embedded in the complete mixing scheme.

The evolution matrix of neutrinos in matter, , obeys the first order (operator) differential equation,

 id S(x,x0)d x=H(x)S(x,x0) , (1)

where the Hamiltonian is given in the flavor basis by

 H=MM†2E+^V=12EU(θ)M2ΔU(θ)†+^V. (2)

Here is the matrix of potentials,

 (3)

is the mixing matrix, and is the diagonal matrix of mass squared differences.

Formally, the solution of the equation (1) can be written as the chronological product

 S(xf,x0)=Te−i∫xfx0H(x) dx≡limn→∞e−iH(xn)Δx⋅e−iH(xn−1)Δx⋯e−iH(x1)Δx , (4)
 Δx=xf−x0n .

In our previous papers, Ioannisian:2004jk (), Ioannisian:2004vv (), we performed expansion of each exponential factor in eq. (4) and then took limit . Such a procedure does not guarantee the unitarity once the series is truncated and finite number of terms of the expansion is taken.

In this paper we will use expansions of powers of exponents and sum up contributions in the power without expansion of exponents themselves. Consequently, the form, , of the -matrix, and therefore, the unitarity are maintained since is a hermitian matrix. The Magnus expansion magnus () has the following form

 S=e−iC[H]≡e−i(C1+C2+C3+...), (5)

where functional is a series in powers of commutators of the Hamiltonians taken in different points of neutrino trajectory. The term contains commutators of order :

 C1 = ∫xfx0dx H(x) , (6) C2 = −i2∫xfx0dx∫xx0dy [H(x), H(y)] , (7) C3 = (−i)26∫xfx0dx∫xx0dy∫yx0dz([H(x),[H(y),H(z)]]+[[H(x),H(y)],H(z)]). (8)

The details of derivation of the functionals are given in the Appendix. The representation of the matrix in eq.(5) with given in (6) (7), (8) is the main tool which we will use for different applications.

The Magnus expansion is an integral version of the Baker-Campbell-Hausdorff (BCH) equality. Recall that according to the BCH-equation, the summation of powers in exponents leads to

 ea⋅eb=ea+b+12[a,b]+112[(a−b)[a,b]]+…,

that is, to appearance of commutator of the operators. In fact, in matter with varying density the Hamiltonians taken in different spatial points do not commute

 [H(xi),H(xj)]≠0.

The calculation of the S-matrix (4) requires an extension of the BCH-equality to a product of many exponential factors, and eventually, a transition to the continuous limit.

### ii.2 Properties of Magnus expansion.

Let us consider general properties of the Magnus expansion given in eqs.(5, 6, 7, 8).

1). If constant, then for , and therefore in the uniform medium the matrix is given by

 S=e−i∫xfx0dx H(x)=e−iH (xf−x0).

This reproduces immediately the standard oscillation results. All corrections (due to the non-constant Hamiltonian) are given by the commutators. Essentially, the Magnus expansion is the expansion in the number of commutators.

2). The terms of the Magnus expansion (6, 7, 8) contain factorials in denominator, therefore a convergence of the series is better than a convergence of the usual expansion (see eq. (89) in the Appendix). The Magnus series has good convergence even if is not small.

3). The commutators themselves may contain an additional smallness. The weaker dependence of on distance the smaller the commutators. So, in a sense, we deal here with a kind of adiabatic expansion.

4). If is a symmetric function with respect to the middle point of a neutrino trajectory,

 ¯x=xf+x02,

that is,

 H(x)=H(2¯x−x), (9)

one can show that () burum (), and only the odd terms in the expansion are non-zero. Let us prove that (general proof is given in burum ()). According to eq. (7) the integration region is symmetric with respect to the diagonal line , that is, symmetric under reflection:

 (x, y)→(2¯x−y, 2¯x−x) (10)

(). Taking into account the symmetry of Hamiltonian (9) it is easy to show that under the reflection (10) the commutator changes the sign. Therefore the integration of this commutator gives zero.

### ii.3 Magnus expansion in the “interaction” representation

Let us split the total Hamiltonian into two parts

 H(x)=H0(x)+Υ(x) (11)

in such a way that is self-commuting along a trajectory. That is, for any two points of the trajectory , : . The rest of the Hamiltonian, , is not self-commuting, in general, and if small can be treated as a perturbation. In this case it is convenient to solve the problem in the basis of new states, , related to the initial basis by

 ψ=UI(x)ψI=e−i∫xx0dtH0(t)ψI. (12)

Inserting this relation into the evolution equation we find that , and the corresponding matrix, satisfy the evolution equation with the Hamiltonian , where

 ΥI(x,x0)=U†IΥ(x)UI(x)=ei∫xx0H0(t) dtΥ(x)e−i∫xx0H0(t)dt. (13)

The transformation to new basis (12) is equivalent to transition to a “interaction representation” if is interpreted as the Hamiltonian of free propagation. can be considered as an operator in the interaction representation.

The evolution matrix in the interaction representation is given by

 SI(xf,x0)=e−iC[ΥI(x,x0)], (14)

that is, in the formulas (6, 7, 8) one should substitute . Then, according to eq. (12) the matrix in the original basis equals

 S(xf,x0)=UI(xf)SI(xf,x0)UI(x0)†, (15)

or explicitly,

 S(xf,x0)=e−i∫xfx0dtH0(t)e−iC[ΥI(x,x0)]. (16)

(The exponent on the RH side of this equality disappears because of the integration limits.) If , so that it can be considered as a small perturbation, a convergence of the series will be fast.

The Hamiltonian is self-commuting if its dependence on distance can be factorized:

 H0(x)=f(x)⋅M, (17)

here is an arbitrary function of and is an arbitrary constant matrix. Specific realizations of (17) include constant (-independent) Hamiltonians as well as the diagonal Hamiltonians . In the latter case subtracting a matrix proportional to the unit matrix: , one can reduce the Hamiltonian to the form (17).

In the case of small mixing (which can be achieved selecting certain basis of neutrino states) we can split the Hamiltonian as

 H(x)=Hdiag(x)+Hoff−diag(x)

and identify with .

### ii.4 Evolution in symmetric potential

Let us consider a symmetric density profile so that the Hamiltonian satisfies the equality (9). In this case it is convenient to perform the integration in from the middle point of neutrino trajectory, , and to choose the evolution basis , such that with

 ¯UI(x)=e−i∫x¯xdtH0(t). (18)

Essentially here we have substituted by . Now (similarly to the consideration in the previous section) the evolution matrix can be written as

 SI(xf,x0)=e−iC[ΥI(x,¯x)], (19)

where

 ΥI(x,¯x)=ei∫x¯xH0(t) dtΥ(x)e−i∫x¯xH0(t)dt. (20)

Then, the evolution matrix in the original basis equals

 S(x,x0)=¯UI(x)SI(x,x0)¯UI(x0)†, (21)

or explicitly, for an evolution from to we obtain

 S(xf,x0)=e−i∫xf¯xdtH0(t)e−iC[ΥI(x,¯x)]e−i∫¯xx0dtH0(t). (22)

Notice that in contrast to the operator has no definite symmetry with respect to the middle of a trajectory even for a constant density profile. Therefore the even coefficients, , are non-zero:

 ¯C1≡C1[ΥI(x,¯x)] = ∫xfx0dxΥI(x), ¯C2≡C2[ΥI(x,¯x)] = −i12∫xfx0dx ∫xx0dy[ΥI(x),  ΥI(y)], (23)

etc.. Here “bar” indicates that have been calculated in the interaction representation with the -matrix integrated from the middle point of trajectory.

Let us introduce the variable

 r≡x−¯x=x−xf+x02 (24)

which is the distance from the middle of trajectory. Then

 ¯C1 = ∫L−L drΥI(r), (25) ¯C2 = −i2∫L−Ldr∫r−L dp[ΥI(r),  ΥI(p)]. (26)

Here

 L≡xf−x02 (27)

and

 ΥI(r)=ei∫r0H0(t)dtΥ(r)e−i∫r0H0(t)dt. (28)

Notice that the expressions (25, 26, 28) are valid for any density profile and we have not used yet any symmetry of the Hamiltonian.

Let us now assume that , and consequently, the Hamiltonian, are symmetric functions with respect to the middle point of a trajectory, , (as for neutrinos crossing the Earth). In this case and are the even functions of :

 H0(−r)=H0(r),   Υ(−r)=Υ(r). (29)

Denoting

 Φ0≡∫r0H0(t)dt (30)

we have

 Φ0(−r)=−Φ0(r) (31)

provided that is real. Let us show that in this case and are the real symmetric matrices. The proof is straightforward in the case of real . The function is not symmetric with respect to . Indeed, rewriting (28) as

 ΥI(r)=eiΦ0(r)Υ(r)e−iΦ0(r), (32)

one can see immediately that under

 ΥI(−r)=ΥI(r)∗. (33)

Using this relation and the definition (25) we obtain

 ¯C∗1=∫L−LdrΥI(r)∗=∫L−LdrΥI(−r)=∫L−LdrΥI(r)=¯C1,

where in the last equality we made a substitution . Furthermore, since is Hermitian, , we obtain that , i.e., the matrix is symmetric.

Similarly we can show that . Here in addition to the property (33) and the change of the signs of variables, we use that

Again, since is Hermitian, the matrix should be symmetric.

Performing integration in the expressions for (25, 26) from the middle point of a trajectory we obtain

 ¯C1 = 2∫L0drReΥI(r), ¯C2 = 2∫L0dr∫r0 dp[ImΥI(r),  ReΥI(p)] (34)

from which we immediately conclude that are real.

As we will see in sect. III.C, in the adiabatic perturbation theory is purely imaginary matrix. Moreover, since , for a symmetric potential we have the antisymmetric . So,

 Υ(r)∗=−Υ(r),   Υ(−r)=−Υ(r), (35)

and therefore . Using the equalities (35) one can show that in this case also satisfies the equality (33), and consequently the matrices can be calculated as in eq. (34).

For a symmetric potential using the property (31) we can write the S-matrix in the original basis (22) as

 S(xf,x0)=e−iΦ0(L)e−iC[ΥI(x,¯x)]e−iΦ0(L). (36)

## Iii Oscillation probabilities

In applications of the Magnus expansion, adjusting the formalism to a specific physical situation we can select

• propagation basis, that is, the basis of neutrino states in which we consider evolution;

• split of the Hamiltonian into self-commuting and non-commuting parts;

• perturbation terms.

In what follows we will consider a symmetric density profile keeping in mind applications to the neutrino propagation inside the Earth.

### iii.1 Low energy and low density limit

In the low energy or/and low density case it is convenient to consider the neutrino evolution in the mass eigenstates basis, . In this basis the Hamiltonian can be written as

 (37)

where is the vacuum mixing matrix (3). We split the Hamiltonian, according to (11), in the following way. The self-commuting part can be chosen as

 (38)

where is the difference of the instantaneous eigenvalues of the Hamiltonian (37):

 Δm(x)≡Δm22E ⎷(cos2θ−2EV(x)Δm2)2+sin22θ . (39)

Then, according to (37), the perturbation part equals

 (40)

where

 A(x) ≡ 12sin2θ V(x), B(x) ≡ 12[Δm(x)−Δm22E+V(x)cos2θ]. (41)

For a weak potential : , we have

 B(x)=14(Vsin2θ)22EΔm2+O(V3)≈A2(x)2EΔm2. (42)

According to (38) the matrix of transition to the interaction representation equals

 (43)

where

 ϕ(x)≡∫x0Δm(r) dr (44)

is the adiabatic phase (here the integration runs from the middle point of a trajectory). Then the Hamiltonian in the interaction representation, can be written as

 (45)

Using this expression and eqs. (34) we obtain

 (46)

with

 Z(L) ≡ 2∫L0drA(r)cosϕ(r)+4∫L0dr∫r0dp A(r)B(p)sinϕ(r), Y(L) ≡ 2∫L0drB(r)−4∫L0dr∫r0dp A(r)A(p)sinϕ(r)cosϕ(p). (47)

Let us estimate these quantities with accuracy . Since and , the last term in , being of the order , can be neglected. For the function performing integration by parts in the second integral we have

 Y(L) = 2∫L0dr[B(r)−2A(r)2Δm(r)]+4∫L0dr∫r0dpddr[A(r)Δm(r)]ddp[A(p)Δm(p)]sinϕ(p)cosϕ(r) (48) = sin22θ∫L0dr∫r0dpddr[V(r)Δm(r)]ddp[V(p)Δm(p)]sinϕ(p)cosϕ(r)+O(V3),

where in the last equality we used expression (42).

Neglecting we find

 (49)

and

 IV≡sin2θ∫L0dr V(r)cosϕ(r). (50)

Using eqs. (36), (49) and (43) we obtain the matrix (16) in the mass-eigenstates basis

 S = (51)

Here is the half of the oscillation phase:

 ϕ≡ϕ¯x→xf=ϕx0→¯x=ϕ(L). (52)

Notice that both the matrix that originates from the self-commuting part and the perturbation, , depend on the same adiabatic phase.

For the transition between the mass states we have immediately from (51):

 Pν2→ν1=|S21|2=sin2IV. (53)

The matrix for the mass-to-flavor transitions equals

 Smass−flavor=U(θ)⋅S,

and the probability is

 Pνi→να=|(U⋅S)αi|2. (54)

From (54), (51) and (3) we obtain

 Pν2→νe=sin2θ+12sin2θ sin2IV sinϕ+cos2θ sin2IV, (55)

where the first term is simply projection squared of state onto . Eq. (55) reproduces the formula given in Supanitsky:2007ks (). If , we find making expansion in powers of

 Pν2→νe = sin2θ+IVsin2θsinϕ+I2Vcos2θ (56)

which exactly coincides with our result in Ioannisian:2004vv () (see eq. (15)). In a sense, the result (55) corresponds to a re-summation of certain contributions to the probability. It is the substitution that restores the unitarity. Notice that , where is the integral used as the expansion parameter in Ioannisian:2004vv (). According to the present result (55) the expansion parameter includes also which makes convergence even better in the case of small vacuum mixing. Our present consideration explains also the reason why the second order effect in ref. Ioannisian:2004vv () depends on the same integral .

The matrix for transitions between the flavor states equals

 Sflavor−flavor=U⋅S⋅U†.

In particular, for the channel we obtain

 Pνe→να = cos2IV sin22θ sin2ϕ+12sin2IV sin4θ sinϕ +sin2IV cos22θ (57) = (cosIV sin2θ sinϕ+sinIV cos2θ)2.

In the limit , we have and the first term reproduces the standard vacuum oscillation probability. For small the following form of the probability can be useful:

 Pνe→να=sin22θ sin2ϕ+12sin2IV sin4θ sinϕ +sin2IV(cos22θ−sin22θ sin2ϕ). (58)

The result in the second order of the Magnus expansion can be obtained keeping term proportional to in (46). Straightforward calculations give

 (59)

where . Apparently, the result (51) follows from this expression in the limit , .

### iii.2 Perturbation around average potential V0

Let us consider the same situation as in the previous section but perform the expansion with respect to an average potential . This means that we use the basis of neutrino eigenstates in matter with constant potential , as the propagation basis. These eigenstates are related to the flavor states by the mixing matrix in matter

 νf=U(θm0)νm0, (60)

where is defined in (3) and is the mixing angle in matter with the potential , the angle is given by

 sin2θm(V)=sin2θ√(cos2θ−2EV/Δm2)2+sin22θ. (61)

In the - basis the Hamiltonian equals

 (62)

where is the difference of the eigenvalues in matter with the potential , and

 ΔV(x)≡V(x)−V0.

We split the Hamiltonian into the self-commuting part and the perturbation using the same as in the previous case (38). Then the perturbation equals

 (63)

Consequently, for the matrix we obtain the same expression as in Eq. (49) with substitution , where

 I′V=sin2θm0∫L0ΔV(x) cosϕ(x) dx. (64)

In turn, differs from by the substitutions and .

The matrix in the basis equals

and the phase is defined in (52).

Since , the matrix of the mass-to-flavor transitions equals

 Smass−flavor=U(θm0)⋅Sm0⋅U(θm0−θ)†.

Then the probability is

 Pν2→νe=cos2I′V[sin2θ+sin2θm0sin2(θm0−θ)sin2ϕ]+ +12sin2I′Vsin2(2θm0−θ)sinϕ+sin2I′Vcos2(2θm0−θ). (65)

Apparently this expression is reduced to the one in eq. (55), if .

For the mass-to-mass transitions the matrix equals

 Smass−mass=U(θm0−θ)⋅Sm0⋅U(θm0−θ)†,

and therefore the probabilities are given by the same expressions as for the flavor-to-flavor transitions in the previous section with the substitutions and :

 Pν2→ν1 = cos2I′Vsin22(θm0−θ)sin2ϕ+12sin2I′Vsin4(θm0−θ)sinϕ+sin2I′Vcos22(θm0−θ) (66) = [cosI′V sin2(θm0−θ) sinϕ+sinI′V cos2(θm0−θ)]2.

For the flavor-to-flavor transition we have

 Sflavor−flavor=U(θm0)⋅Sm0⋅U(θm0)†.

Consequently, the probability follows immediately from (66) substituting :

 Pνe→να = cos2I′V sin22θm0 sin2ϕ+12sin2I′V sin4θm0 sinϕ+sin2I′V cos22θm0 (67) = (cosI′V sin2θm0 sinϕ+sinI′V cos2θm0)2.

An interesting feature of the obtained results is that the probabilities for symmetric transitions: the flavor-to-flavor and mass-to-mass ones can be written as a square of the sum of two terms proportional to and .

### iii.3 Adiabatic perturbation theory in Magnus expansion

Let us again consider symmetric density profile. As the propagation basis, we take the basis of the eigenstates of instantaneous Hamiltonian, ,:

 νf=U(θm(x))νm .

Here is the instantaneous mixing angle in matter (61). The Hamiltonian for the eigenstates equals , where

 (68)

and

 ˙θm(x)≡dθm(x)dx=sin2θm(x)2Δm(x) dV(x)dx. (69)

In what follows we will use and as the self-commuting and perturbation parts correspondingly. Notice that the self-commuting part is the same as before, but the perturbation is different since the basis of states differs from the one we used before. Now is a complex and non-symmetric matrix with respect to the middle of trajectory. Straightforward calculations give according to (23) or (34)

 (70)

where

 Iθ = −2∫xf¯x˙θm(x) sinϕ¯x→xdx= (71) = 2∫xf¯x[θm(x)−θms] Δm(x) cosϕ¯x→x dx,
 Iθθ = −∫xfx0dx∫xx0˙θm(x)˙θm(y)sinϕy→xdy= (72) = 4∫xf¯xdx∫x¯x˙θm(x)˙θm(y)sinϕ¯x→ycosϕ¯x→xdy.

Here is the mixing angle at the surface of the Earth. Taking into account (69) one sees that has the same structure as the integral in (48).

Neglecting the second order term , we have which coincides, according to (70), with total in eq. (49) up to the change . Therefore the adiabatic probabilities equal to those in the previous subsection with the substitutions , :

 Pν2→νe = cos2Iθ sin2θ +cos2Iθ sin2θms sin2(θms−θ) sin2ϕ + (73) + 12sin2Iθ sin2(2θms−θ) sinϕ +sin2Iθ cos2(2θms−θ),
 Pν2→ν1 = cos2Iθ sin22(θms−θ) sin2ϕ+ 12sin2Iθ sin4(θms−θ) sinϕ (74) + sin2Iθ cos22(θms−θ) ,
 Pνe→να = cos2Iθ sin22θms sin2ϕ+ 12sin2Iθ sin4θms sinϕ+ (75) + sin2Iθ cos22θms.

Notice that , when .

Let us take into account the second order of the Magnus expansion. Now contains the term proportional to the diagonal matrix. Apparently, has the same form as in (46) with the substitutions and . So, using the results (59), (70), (18) and (22), we find the -matrix in the basis of eigenstates of the Hamiltonian,

 Sm=⎛⎜⎝cosIt−iIθθItsinIt−iIθItsinIte−iϕ−iIθItsinIte−iϕ[cosIt+iIθθItsinIt]e−2iϕ⎞⎟⎠. (76)

Here

 It≡√I2θ+I2θθ, (77)

and the adiabatic phase is defined in (52). The -matrix for the flavor-to-flavor transitions is then given by

 Sflavor−flavor=U(θms)⋅Sm⋅U†(θms). (78)

For the probability of oscillations,