High-energy neutrino oscillations in absorbing matter

The impact of neutrino mixing, refraction and absorption on high-energy neutrino propagation through a thick medium is studied using the MSW evolution equation with complex indices of refraction. It is found that, owing to the mixing with sterile neutrinos, the penetrability of active neutrinos may be many orders of magnitude larger than it would be in the absence of mixing. The effect is highly sensitive to changes in density and composition of the matter background as well as to neutrino energy and mixing parameters. This may lead to observational consequences in neutrino astrophysics.


Introduction
In studies of high-energy neutrino propagation through matter two completely different approaches are in use: the quantum-mechanical MSW formalism [1] and the classical transport theory [2]. The first approach accounts for the neutrino mixing and coherent interference of propagating neutrinos (that is refraction) neglecting the loss of coherence due to inelastic scattering. The second one provides the means for estimating neutrino absorption, production and energy loss in inelastic charged and neutral current interactions, but is inapplicable for the inclusion of quantum effects of neutrino mixing and refraction. In other words, the necessary conditions for the validity of these approaches are mutually exclusive. At the same time there is a multitude of physically feasible and interesting contexts where the effects of neutrino mixing, refraction and inelastic collisions may be equally important or comparable. Well-known examples are the propagation of high-energy neutrinos produced in the sun from cosmic ray interactions and from annihilations of heavy dark matter particles captured by the sun (see, e.g., ref. [3] and references therein). The so-called hidden sources (the astrophysical neutrino sources opaque for electromagnetic radiation) [4] hold promise as the splendid "testing laboratories" for the effects under discussion because some types of these objects (e.g., a Thorne-Żytkow star [5]) have relevant dimensions and density profiles.
The necessity of a unified description of mixed neutrinos in hot media (like the CP symmetric plasma of the early Universe, a supernova core or a protoneutron star) has long been realized [6] and the rigorous quantum kinetic theory has been developed [7]. 3 However, the direct application of that theory to the problem of high-energy neutrino transport in a normal cold matter is very involved and the analytic treatment of the quantum kinetic equations [7] is hardly possible, except for very particular cases.
Since an analytic analysis is always useful to gain an insight into a sophisticated problem, it seems instructive to apply the conventional MSW formalism straightforwardly generalized by incorporating the imaginary parts into the neutrino indices of refraction. This will allow an understanding of the impact of the neutrino absorption in conjunction with mixing and refraction, still disregarding the neutrino energy loss from the neutral current interactions [2,10], charged current induced chains like ν τ N → τ X, τ → ν τ X ′ [11] and neutrino flavor changing reactions ν e e − → ν µ µ − , ν µ e − → ν e µ − , etc. [12]. The present paper is concerned with some surprising phenomena which follow from this simple approach. For further simplification we limit ourselves to studying two-flavor case. Some general results relevant to the subject were discussed previously in ref. [13].

Master equation
The generalized MSW equation for the time-evolution operator of two mixed stable neutrino flavors ν α and ν β propagating through a matter with absorption can be written as Here is the vacuum mixing matrix that relates flavor to mass eigenstates, θ is the mixing angle (0 ≤ θ ≤ π/2), and m i are the energies and masses of the neutrino mass eigenstates, respectively, 4 n α (t) is the complex index of refraction, where N 0 = 6.022 × 10 23 cm −3 , f ναk (0) is the amplitude for the ν α zero-angle scattering from particle k (k = e, p, n, . . . ), ρ(t) is the density of the matter (in g/cm 3 ) and Y k (t) is the number of particles k per amu in the point t of the medium. From the optical theorem (see, e.g., ref. [14]) we have Im [f ναk (0)] = (p ν /4π)σ tot ναk (p ν ), where σ tot ναk (p ν ) is the total cross section for ν α k scattering.
This implies that where Λ α (t) is the mean free path of neutrino ν α in the point t of the medium.
It is useful to transform eq. (1) into the one with a traceless Hamiltonian. For this purpose we define the matrix After substituting eq. (4) into eq. (1), we have Here The Hamiltonian for antineutrinos is of the same form as eq. (6) but one must keep in mind that Re [f ναk (0)] = −Re [f ναk (0)], while σ tot ναk (p ν ) and σ tot ναk (p ν ) are different in magnitude even at very high energies.
The neutrino oscillation probabilities are just the squared absolute values of the elements of the evolution operator S(t), Taking into account eqs. (2), (3), (4) and (7) yields where Owing to the complex potential q, the Hamiltonian (6) is non-Hermitian and the new evolution operator (4) is nonunitary. As a result, there are no conventional relations between the four probabilities given by eq. (8). Since the matrix H(t) becomes Hermitian when Λ α = Λ β . If this is the case at any t, eq. (5) reduces to the standard MSW equation and inelastic scattering in the medium results in the common exponential attenuation of the survival and transition probabilities. From here, we shall consider the more general and more interesting case, when Λ α = Λ β .
The extreme example of such a case is provided by mixing between the ordinary (α = e, µ or τ ) and sterile (β = s) neutrino flavors. Since Λ s = ∞, we have Λ = 2Λ α and q I = −1/4Λ α . So q I is nonzero at any energy. Even without solving the evolution equation, one can expect the penetrability of active neutrinos to be essentially modified in this case because, roughly speaking, they spend a certain part of life in the sterile state.
Other important examples are the ν e − ν τ and ν µ − ν τ mixings. The total cross sections for e or µ production are well above the one for τ production within a wide energy range [15]. This is because of large value of τ lepton mass m τ which leads to high thresholds and sharp shrinkage of the phase spaces for the charged current ν τ N reactions as well as to the kinematic correction factors (∝ m 2 τ ) to the nucleon structure functions; the corresponding structures are negligible in the case of electron production and small for muon production. The neutral current contributions are, of course, cancelled out from q I . Thus, in the context of the master equation (5), ν τ can be treated as (almost) sterile within the energy range for which σ CC νe,µN ≫ σ CC ντ N .
A similar situation, while in quite a different and narrow energy range, holds in the case of mixing of ν e with some other flavor. This is a particular case for a normal C asymmetric medium, because of the W boson resonance formed in the neighborhood of E res ν ≈ 6.33 PeV through the reactions Just at the resonance peak, σ tot νee ≈ 250 σ tot νeN (see ref. [12] for further details and references).
For E ν ≪ min m 2 W,Z /2m k (where m k,W,Z denote the masses of the scatterers k and gauge bosons) and for an electroneutral nonpolarized cold medium, the real part of the potential q is energy independent. In the leading orders of the standard electroweak theory it is [16] where G F is the Fermi coupling constant, α is the fine-structure constant, θ W is the weak-mixing angle and r τ = (m τ /m W ) 2 . Thus, for an isoscalar medium |q R | is of the same order of magnitude for any pair of flavors but ν µ − ν τ ; the ratio q For certain regions of a neutronrich medium the value of q (νe−νs) R may become vanishingly small. In this case, the one-loop radiative corrections must be taken into account. For very high energies, all the expressions (9) have to be corrected for the gauge boson propagators and strong-interaction effects.
One can expect |q R | to be either an energy-independent or decreasing function for any pair of mixed neutrino flavors. On the other hand, as noted above, there are several cases of much current interest when |q I | either increases with energy without bound (mixing between active and sterile neutrino states) or has a broad or sharp maximum (as for ν µ −ν τ or ν e −ν µ mixings, respectively). Numerical estimations, performed using the results of refs. [12] and [15], suggest that for every of these cases there is an energy range in which q R and q I are comparable in magnitude. Considering that both q R and q I are proportional to the density and besides are dependent upon the composition of the medium there may exist even more specific situations, when |q R | ∼ |q I | ∼ |∆| or |q R | ∼ |∆ c | and |q I | ∼ |∆ s |. If this is the case, the refraction, absorption and mixing become interestingly superimposed.

Eigenproblem and mixing matrix in matter
The matrix (6) has two complex instantaneous eigenvalues, ε(t) and −ε(t), with ε = ε R + iε I satisfying the characteristic equation where q ± = ∆ c ± i∆ s = ∆e ±2iθ . The solution to eq. (10) for ε R and ε I may be written as with The sign of ε R is a matter of convention. We define sign (ε R ) = sign(∆) ≡ ζ.
Let us consider the behavior of ε R and ε I in the vicinity of the MSW resonance, defined by the condition q R = q R (t ⋆ ) = ∆ c (we suppose for simplicity that the matter density and composition are continuous functions of t). The accurate passage to the limit in eqs. (11) gives where ∆ I is taken to be the value of q I in the MSW point (that is ∆ I = q I (t ⋆ )). Therefore the resonance value of |ε R | (which is inversely proportional to the neutrino oscillation length in matter) is always smaller than the conventional MSW value |∆ s | and vanishes if ∆ 2 I < ∆ 2 s (ε I remains finite in this case). In neutrino transition through the region of resonance density ρ = ρ(t ⋆ ), ε I undergoes discontinuous jump while ε R remains continuous. The corresponding cuts in the q plane are placed outside the circle |q| ≤ |∆| as is shown in fig. 1. If ∆ 2 I > ∆ 2 s , the imaginary part of ε vanishes while the real part is kept finite.
A distinctive feature of eq. (10) is the existence of two mutually conjugate "super-resonance" points q ± in which ε vanishes giving rise to the total degeneracy of the levels of the system under consideration. Certainly, the behavior of the system in the vicinity of these points must be dramatically different from the conventional pattern. As noted in sec. 2, the "super-resonance" conditions are physically realizable for various meaningful mixing scenarios.
In order to simplify the solution to the eigenstate problem we will assume that the phase trajectory q = q(t) does not cross the points q ± at any t.
Within the framework of the non-Hermitian quantum dynamics (see, e.g., ref. [17] and references therein) one has to consider the two pairs of instantaneous eigenvectors |Ψ ± and |Ψ ± which obey the relations For q = q ± these pairs form a complete biorthogonal and biorthonormal set, Therefore, the eigenvectors are defined up to a gauge transformation with arbitrary complex functions f ± (t) such that Im (f ± ) vanish as q = 0. 5 Thus it is sufficient to find any particular solution of eqs. (12). Taking into account that H † = H * , we may set |Ψ ± = |Ψ * ± and hence the eigenvectors can be found from the identity a particular solution of which can be written as where and we have fixed the remaining gauge ambiguity by a comparison with the vacuum case.
It may be sometimes useful to define the complex mixing angle in matter Θ by the relations sin Θ = v + and cos Θ = v − or, equivalently, The real and imaginary parts of Θ are found to be Having regard to the above-mentioned prescription for the sign of ε R , one can verify that Θ = θ if q = 0 (vacuum case) and Θ = 0 if ∆ s = 0 (no mixing or m 2 1 = m 2 2 ). It is also clear that Θ becomes the standard MSW mixing angle with Im(Θ) = 0 when q I = 0 (Λ α = Λ β ).
In order to build up the solution to eq. (5) for the nondegenerate case one has to diagonalize the Hamiltonian H. Generally a non-Hermitian matrix cannot be diagonalized by a single unitary transformation. But in our simple case this can be done by a complex orthogonal matrix (extended mixing matrix in matter) where f = diag (f − , f + ), f ± are the arbitrary complex phases mentioned above and The matrix U by its construction has the following properties: From eq. (10) it follows that ∂ε/∂q = (q − ∆ c )/ε and thus eqs. (13) yield We therefore have where According to eqs. (15) and (16), iU

Adiabatic solution
With the mixing matrix U f in hand we can write the formal solution to eq. (5) in the most general form: Here X f (t) is an unknown matrix, Φ(t) = diag (−Φ(t), Φ(t)) and Φ(t) = Φ R (t) + iΦ I (t) is the complex dynamical phase, defined by Substituting eq. (18) into eq. (5) and using eqs. (17), we find It can be proved now that the right side of eq. (18) is gauge-invariant i.e. it does not depend on the unphysical complex phases f ± (t). This crucial fact is closely related to the absence of the Abelian topological phases in the system under consideration. 6 Finally, we can put f ± = 0 in eq. (18) and the result is These equations, being equivalent to eq. (5), have nevertheless a restricted range of practical usage on account of poles and cuts as well as decaying and increasing exponents in the "Hamiltonian" ΩF. However, the adiabatic theorem of Hermitian quantum mechanics (see, e.g., ref. [19]) can straightforwardly be extended to eq. (5) under the following requirements: (a) the potential q is a sufficiently smooth and slow function of t; (b) the imaginary part of the dynamical phase is a bounded function i.e. lim t→∞ |Φ I (t)| is finite; (c) the phase trajectory q = q(t) is placed far from the singularities for any t.
The first requirement breaks down for a condensed medium with a sharp boundary or layered structure (like the Earth). If however the requirement (a) is valid inside each layer (t i , t i+1 ), the problem reduces to eqs. (19) by applying the rule where S (t i+1 , t i ) is the time-evolution operator for the i-th layer.
The requirement (b) alone is not too restrictive considering that for many astrophysical objects (like stars, galactic nuclei, jets and so on) the density ρ exponentially disappears to the periphery and, on the other hand, ε I → 0 as ρ → 0. In this instance, the function Φ I (t) must be t independent for sufficiently large t. But, in the case of a steep density profile, the requirements (a) and (b) may be inconsistent.
The important case of violation of the requirement (c) is the subject of a special study which is beyond the scope of this paper. It is interesting to note in this connection that, in the Hermitian case, a general adiabatic theorem has been proved without the traditional gap condition [20].
Having regard to the mentioned limitations, we can expect that eqs. (19) is tolerable for finding the approximate, adiabatic solution to the problem (and the corrections to that solution) for many cases of pragmatic interest. We shall presume now that the parameters of the matter vary so slowly that all necessary conditions do hold for 0 ≤ t ≤ T . Then, in the adiabatic limit, we can put Ω = 0 in eq. (19b). Therefore X = 1 and eq. (19a) yields Taking into account eq. (8) we obtain the survival and transition probabilities: where we have denoted for compactness In the event that the conditions and are satisfied for any t ∈ [0, T ], the formulas (20) reduce to the standard MSW adiabatic solution where Needless to say either of the conditions (21) or both may be violated for sufficiently high neutrino energies and/or for thick media, resulting in radical differences between the two solutions. These differences are of obvious interest to high-energy neutrino astrophysics.
It is perhaps even more instructive to examine the distinctions between the general adiabatic solution (20) and its "classical limit" 7 , P αβ (t) = 0, , P βα (t) = 0, which takes place either in the absence of mixing or for m 2 1 = m 2 2 .
In the next section we consider some features of the solution (20) for media with constant density and composition (q ≡ 0). For this simple case, the adiabatic approximation becomes exact and thus free from the above-mentioned conceptual difficulties.
For definiteness sake we assume Λ α < Λ β (and thus q I < 0) from here. The opposite case can be considered in a similar way. Let us denote As is easy to see, and sign(ϕ) = −ζ.

Case |q| |∆ s |
Let us examine the case when Λ + and Λ − are vastly different in magnitude. This will be true when Λ β ≫ Λ α and the factor ξ is not too small. The second condition holds if q R is away from the MSW resonance value ∆ c and the following dimensionless parameter is sufficiently small. In fact we assume |κ| 1 and impose no specific restriction for the ratio q R /q I . The assumption spans several possibilities: (i) small ∆m 2 , (ii) small mixing angle, (iii) high energy, (iv) high matter density. The last two possibilities are of special interest because the inequality |κ| 1 may be fulfilled for a wide range of the mixing parameters ∆m 2 and θ by changing E ν and/or ρ. In other words, this condition is by no means artificial or too restrictive.
After simple while a bit tedious calculations we obtain Since Λ α has been assumed to be small compared to Λ β , we have In general the oscillation length in matter L has no Taylor expansion over the parameter κ and may vary within a broad range depending on other circumstances. For example, L ≈ π/|q| if |q| ≫ |∆| and L ≈ π/ ε 2 0 − ∆ 2 s if |q I | ≫ ε 0 > |∆ s |. Due to the wide spread among the length/time scales Λ ± , Λ and L as well as among the amplitudes I ± and I, the regimes of neutrino oscillations are quite diverse for different ranges of variable t.
Some of the essential features are illustrated in figs. 2 and 3 by the example of ν µ − ν s oscillations in an isoscalar medium (Y p = Y n = 0.5). For such a medium the picture of the ν e − ν s (ν τ − ν s ) oscillations is the same (nearly the same). In all examples, the mixing parameters are taken to be ∆m 2 = 10 −3 eV 2 and θ = π/4. In order to estimate the charged and neutral current ν µ N total cross sections we apply the results of ref. [12] based on the CTEQ4-DIS parton distributions [21]. The ν µ e contribution is neglected. Figure 2 depicts the survival and transition probabilities for E ν = 1 TeV and ρ = 0.2 g/cm 3 (κ ≈ 0.131 in this case). The classical limit of P µµ (that is the penetration coefficient for the unmixed muon neutrino flux) is also shown for comparison. In fig. 3 we present the same probabilities but for E ν = 100 TeV and (in order to demonstrate the density impact) for the two values of ρ, 10 −3 g/cm 3 (upper panel) and 3 × 10 −4 g/cm 3 (lower panel). In these cases κ ≈ 0.261 and 0.870, respectively. With reference to figs. 2 and 3, one can see a regular gradation from slow (at t Λ µ ) to very fast (at t Λ µ ) neutrino oscillations followed by the asymptotic nonoscillatory behavior The latter regime is the most remarkable. Since Λ − ≫ Λ µ , a small part of the flux of muon neutrinos escapes through the medium of such a huge thickness that would be absolutely opaque for them in the absence of mixing. At t > (10 − 20)Λ µ the deviation of the survival probability P µµ from the classical expectation grows exponentially. The asymptotic value of the transition probability P sµ is a factor of 4/κ 2 in excess of P µµ . This may be potentially beneficial for the observational neutrino astrophysics, considering that after leaving the medium surrounding an astrophysical neutrino source, the sterile neutrinos may travel a very long distance in vacuum giving rise (through the vacuum oscillations) to a detectable amount of muon neutrinos.

Degenerate case
Our consideration must be completed for the case of degeneracy. Due to the condition q I < 0, the density and composition of the "degenerate environment" are fine-tuned in such a way that q = q −ζ = ∆ c − i |∆ s |. The corresponding formulas can be derived by a formal passage to the limit in eqs. (24). A more simple way is however in coming back to the master equation (5). Indeed, in the limit of q = q −ζ , the Hamiltonian (6) reduces to Considering that h 2 ζ = 0, we promptly arrive at the solution of eq. (5): and thus, taking into account eq. (8), we have P αα (t) = (1 − |∆ s | t) 2 e −t/Λ , P ββ (t) = (1 + |∆ s | t) 2 e −t/Λ , P αβ (t) = P βα (t) = (∆ s t) 2 e −t/Λ .

Conclusions
We have considered, on the basis of the MSW evolution equation with complex indices of refraction, the conjoint effects of neutrino mixing, refraction and absorption on high-energy neutrino propagation through matter. The adiabatic solution with correct asymptotics in the standard MSW and classical limits has been derived. In the general case the adiabatic behavior is very different from the conventional limiting cases.
A noteworthy example is given by the active-to-sterile neutrino mixing. It has been demonstrated that, under proper conditions, the survival probability of active neutrinos propagating through a very thick medium of constant density may become many orders of magnitude larger than it would be in the absence of mixing. The quantitative characteristics of this phenomenon are highly responsive to changes in density and composition of the medium as well as to neutrino energy and mixing parameters. Considering a great variety of latent astrophysical sources of high-energy neutrinos, the effect may open a new window for observational neutrino astrophysics.