{'fraction_non_alphanumeric': 0.041069772235391114, 'fraction_numerical': 0.01556219725857982, 'mean_word_length': 4.439181614349776, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exact series solution to the two flavor neutrino oscillation problem in matter

In this paper, we present a real non-linear differential equation for the two flavor neutrino oscillation problem in matter with an arbitrary density profile. We also present an exact series solution to this non-linear differential equation. In addition, we investigate numerically the convergence of this solution for different matter density profiles such as constant and linear profiles as well as the Preliminary Reference Earth Model describing the Earth's matter density profile. Finally, we discuss other methods used for solving the neutrino flavor evolution problem.


I. INTRODUCTION
In general, there are several phenomena and processes in physics, but also in other fields of science such as chemistry, that can be described in terms of a system with two (quantum mechanical) states and a time-dependent Hamiltonian, i.e., so-called two-level systemsneutrino oscillations with two flavors being one such system. Other representatives of such systems are, for example, a spin 1/2 particle in a time-dependent electromagnetic field having the states 'spin up' and 'spin down', K 0 -K 0 mixing, a Josephson device, nuclear magnetic resonance used for encoding bits of information (i.e., quantum bits for a quantum computer), the left and right chirality states of molecules in chemistry, etc. The problem of neutrino oscillations in matter, which we are concerned with in this paper, is mathematically equivalent to a spin 1/2 particle in a magnetic field that is constant in one direction, zero in another direction, and time-dependent in the last direction [1].
Neutrino oscillations have recently been extensively studied in the literature [2,3,4,5,6] and they act as the most plausible description of both the solar [2] and atmospheric [3] neutrino problems. At an early stage, neutrino oscillations were mainly investigated with two flavors and without including matter effects. Nowadays, we know that there are at least three neutrino flavors and that matter effects are important. For example, in matter, the so-called Mikheyev-Smirnov-Wolfenstein (MSW) effect [7,8] can take place, which is an amplifying resonant effect due to the presence of matter. However, in most situations, neutrino oscillations can be effectively investigated with two flavors, since the leptonic mixing in the 1-3 sector is indeed small [4,5] leading to the fact that the full three flavor scenario can be decoupled into two effective two flavor scenarios, each of which can be studied separately.
In this paper, we present an exact analytic solution to the two flavor neutrino oscillation problem in matter. Since there are many similar two-level systems, as discussed above, our solution will also be interesting and applicable to this kind of systems. However, before we proceed to present our solution, we will give a brief overview of what has previously been done in this field. Note that this overview is not presented in chronological order.
First, in Refs. [9,10], the neutrino flavor evolution has been investigated by a discretization of the effective potential. Second, exact solutions exist for a number of specific effective potentials [7,11,12,13]. Third, in Refs. [14,15], the evolution was studied by using an adiabatic approximation. Fourth, approximate solutions valid for small effective potentials have recently been studied in detail [16]. Finally, there have also been other attempts to write the evolution in terms of a second order non-linear ordinary differential equation [17].
However, this has been done for the neutrino oscillation probability amplitudes and not for the neutrino oscillation probabilities. The advantage of working with a non-linear differential equation for the oscillation probability rather than a linear system of differential equations for the probability amplitudes is that we only have one real variable instead of two complex variables. The disadvantage is that the resulting differential equation is non-linear. This paper is organized as follows. In Sec. II, the neutrino flavor evolution in matter with two flavors is studied and a second order non-linear ordinary differential equation for the neutrino oscillation probability is derived. Then, in Sec. III, we perform series expansions of both the neutrino oscillation probability and the effective potential in order to solve the differential equation presented in Sec. II. Next, in Sec. IV, we continue by studying the numerical convergence of the solution for a number of different effective potentials and baselines. In Sec. V, we present a brief summary of other methods for solving the neutrino evolution in matter. Finally, in Sec. VI, we summarize our results and give our conclusions.

II. NEUTRINO FLAVOR EVOLUTION IN MATTER
When neutrinos propagate in matter, neutrino flavors are affected differently by coherent forward scattering against the matter constituents. Assuming that there are no sterile neutrinos, the effect of matter is to add an effective potential to ν e , this effective potential is given by where G F is the Fermi coupling constant and N e (t) is the electron number density.
In the two flavor case, the time evolution of a neutrino state |ν(t) = (|ν e (t) |ν x (t) ) T is given by where H vac = Udiag(m 2 1 , m 2 2 )U † /2E is the free Hamiltonian in vacuum, H mat (t) = diag(V (t), 0) is the addition to the free Hamiltonian due to matter effects, and is the leptonic mixing matrix in vacuum. Here c ≡ cos θ, s ≡ sin θ, and θ is the leptonic mixing angle. Adding or subtracting terms proportional to the unity operator to the total Hamiltonian H(t) = H vac +H mat (t) will only contribute with an overall phase to the neutrino state |ν(t) , and thus, does not affect the neutrino oscillation probabilities. Using this fact, the total Hamiltonian may be written as where the σ i 's (i = 1, 2, 3) are the Pauli matrices and ∆m 2 ≡ m 2 2 − m 2 1 is the mass squared difference between the two mass eigenstates in vacuum.
The density matrix ρ(t) = |ν(t) ν(t)| can be parameterized as ρ(t) = (1 + S(t) · σ)/2, where 1 is the unity matrix, σ = (σ 1 σ 2 σ 3 ) T is the vector of Pauli matrices, and S(t) is a vector such that S(t) 2 = 1. Differentiating the density matrix ρ with respect to time t, the equation of motion for S(t) becomesṠ where B(t) ≡ ge 1 + f (t)e 3 and we have defined g ≡ − sin(2θ)∆m 2 /2E and f (t) ≡ cos(2θ)∆m 2 /2E − V (t). Note that g is independent of time t. The probability of neutrinos produced as ν e to oscillate into ν x (where ν x is some linear combination of ν µ and ν τ ) is now given by P (ν e → ν x ) ≡ P ex = (1 − S 3 (t))/2. With the parameterization where α = α(t) and β = β(t), we obtain the following non-linear system of ordinary differential equationsβ Eliminating β from the above expressions, we obtain the differential equation where Now, we make the substitution p = S 3 (t) = cos α after which Eq. (8) becomes Note that F = 0 corresponds to the so-called MSW resonance condition cos(2θ)∆m 2 /2E = V . In this case, Eq. (9) takes the simple form with the trivial solutions p = A cos(gt) + B sin(gt) just as expected.
In general, the expression for P ex is known for constant matter density and is given by whereθ is the effective leptonic mixing angle in matter and ∆m 2 is the effective mass squared difference in matter. Using that p = 1 − 2P ex , it is a matter of trivial computation to show that this is the solution to Eq. (9) with constant F , which corresponds to any constant matter density.
In the three (n) flavor case, the density matrix can be parameterized by four [2(n − 1)] real parameters. If we would adopt our approach to the three (n) flavor case, then we would end up with a system of seven [2(n − 1)] non-linear ordinary differential equations, which, in principle, can be solved in a manner analogous to the one described above for the two flavor case.

III. SERIES EXPANSION OF THE SOLUTION
In order to solve the propagation of neutrinos in matter with arbitrary density profiles, we adopt the method of series expansion. We suppose that neutrinos are produced as ν e and then propagate through a given effective potential V (t), this gives the initial values p(0) = 1 andṗ(0) = 0. Series expanding the effective potential V (t) and the quantity p(t), we obtain the following expressions where the coefficients V n (n = 0, 1, . . .) define the effective potential and where we wish to compute the coefficients p n (n = 0, 1, . . .). By using the relation between f and V , we obtain Inserting the above expressions into Eq. (9) and identifying terms of the same order in t gives the relation For n = 0 with the given initial conditions, Eq. (16)  Also the solution for n = 3 is now trivially fulfilled, since the terms including p n+1 also cancel for n ≥ 3. For n = 4, the equation is a second order equation in p 4 with the solutions Of these two solutions, only the latter will be a solution to our problem, this is easily checked by inserting the known solution in the case of constant effective potential from Eq. (11).
For n ≥ 5, Eq. (16) is now linear in p n . For n ≥ 6, we obtain a solution for p n in terms of lower order p k , G, and F s , where k < n and s ≤ n − 4. This expression is the following recurrence relation For the first few coefficients we obtain p 0 = 1, As can be observed by setting F k = 0, the solution for F = 0, i.e., at the MSW resonance, is just the series expansion for p = cos(gt), which is clearly as expected.

IV. CONVERGENCE OF THE SOLUTION
In order to test our solution, we perform a number of numerical tests. First of all, we give an overview of how we approximate the electron number density (i.e., in principle, the effective potential) by a polynomial. Then, we proceed by confirming that our solution really converges nicely towards the simple trigonometric function that is the exact solution for a constant electron number density. In this case, we also study the convergence of the energy dependence of the neutrino oscillation probability P ex (L) for a baseline of L = 3000 km.
After the constant electron number density case, we investigate the case of a linear effective potential, and finally, we study the case of the Preliminary Reference Earth Model (PREM) [18].
In the numerical calculations, we have used the mixing angle θ = 13 • and the mass squared difference ∆m 2 = 2 × 10 −3 eV 2 [19]. The value of θ approximately corresponds to the upper limit on the leptonic mixing angle θ 13 from the CHOOZ experiment with ∆m 2 = 2 × 10 −3 eV 2 [5]. The reason to use this particular choice of parameters is that θ 13 and the large mass squared difference give the main effects to neutrino oscillations from ν e into other flavors for the baselines and energies we have studied (for L = 3000 km and E = 1 GeV, the neutrino oscillations governed by the small mass squared difference contribute with an approximate addition of 0.05 to the neutrino oscillation probability P ex , for shorter baselines and higher energies, this effect decreases), see for example Refs. [20]. The reason for using the upper bound value for the mixing angle θ and not some smaller value is that we wish to study the behavior of our solution rather than to make any precise predictions about the neutrino oscillation probabilities.

A. Series expansion of the effective potential
In order to use the series solution, which was obtained in the previous section, we will need the coefficients V n . In general, for a given baseline length L, the effective potential V (t) can be expanded in terms of Legendre polynomials, i.e., where x ≡ 2t/L − 1. For numerical treatments, we cannot use the entire expansion in Legendre polynomials because of finite computer memory and finite computer time. However, if we assume that the coefficients c n are negligible for n > N, where N is some integer, then we have a polynomial approximation of the effective potential. Clearly, given any polynomial V (t), it is a trivial matter to extract the coefficients V n . This approach turns out to be quite handy in the case of the PREM profile, which is discussed below.

B. Constant matter density
The first case we study numerically is the case of a constant effective potential. We use the baseline length L = 3000 km and the electron number density N e = N e,core /3, where V core = √ 2G F N e,core ≃ 5.6 × 10 −19 MeV corresponds to a matter density of about 13 g/cm 3 , which is the maximum matter density in the Earth's core [18]. In this case, the coefficients V n are easily obtained as V 0 = V (t) and V n = 0 for n > 0. Since the exact solution to this problem is known [7], we focus on the convergence of our solution for P ex , both in the energy spectrum and the time evolution. The numerical results are shown in Fig. 1. In this figure, we can observe that approximately 20 terms are needed to reconstruct one period of oscillation and that the convergence is indeed the same as for a simple trigonometric function.

C. Linearly varying matter density
Now, we turn our interest towards the case of a linearly varying effective potential. In particular, we study a baseline of L = 3000 km, where the electron number density is given by N e (t) = N e,core (t/L + 1)  the inner product In Fig. 3, we plot the electron number density profiles for the baseline lengths L = 250 km, 750 km, 3000 km, and 5000 km along with the second order polynomial approximations and the constant approximations. In Fig. 4, we plot the energy spectra and time evolution at E = 1.2 GeV for a baseline of L = 3000 km for the PREM profile. Again, our solution is not plotted in the shaded region in which it breaks down numerically for the same reasons as discussed previously. In this case, it is apparent that if the number of terms used in the series expansion is large enough, our solution is a significant improvement from the constant matter density approximation.
Clearly, the approximation of using a second order polynomial for the electron number density gives a very good reproduction of the numerical solution (which uses the profiles obtained from the PREM). As can be seen in the time evolution plot, the error made is barely noticeable until approximately one and a half oscillations, i.e., for lower neutrino energies if the baseline length L is kept fixed. This is in good agreement with the results obtained in Ref. [21], where the effective potential is expanded in a Fourier series, as well as Ref. [22], which shows that details of the effective potential that are smaller than the For the PREM profile, we are also interested in a number of other baseline lengths. In particular, in Fig. 5, we plot the energy spectra for the baseline lengths L = 250 km, 750 km, 3000 km, and 5000 km. For the baseline lengths L = 250 km and 750 km, there is no noticeable difference between the numerical solution, our exact solution, and the approximation using constant electron number density. This is to be expected as the electron number density does not vary significantly for these baseline lengths (see Fig. 3). However, for both L = 3000 km and 5000 km, we do observe a difference between the constant electron number density approximations and the other two solutions. Again, we can conclude that the approximation with a second order polynomial for the electron number density agrees remarkably well with the numerical calculation using the profile obtained directly from the PREM. Note that the comparison of the energy spectra for different matter density profiles and baselines has been studied before [23].
For baselines longer than L = 5000 km, the region near the resonance for V = V 0 tends to expand and ruin the numerical convergence of our solution. Also, this region expands if we include more terms of the series expansion. For the baseline lengths below L = 5000 km, this can be somewhat compensated by using fewer terms of the series expansion for high energies to avoid the numerical cancellation effects and more terms for lower energies to obtain a nice convergence.

V. OTHER METHODS OF SOLVING THE NEUTRINO FLAVOR EVOLUTION
In general, there have been numerous methods on how to solve the problem of neutrino flavor evolution in matter [7,9,10,11,12,13,14,15,16,17]. First of all, there is the obvious formal solution using a time-ordered exponential, i.e., This solution is exact, but it is not very helpful in actual calculations due to the nature of the time-ordered exponential. A way of solving this problem is to use a discretization of the effective potential [9] [24]. The effective potential is then divided into a finite number of layers with constant effective potentials (i.e., constant electron number density), which approximate a given effective potential. Clearly, when the number of layers goes to infinity, one regains the time-ordered exponential.
There is also the pure numerical approach to the problem, where the neutrino flavor evolution is easily solved numerically for any effective potential. While this gives the possibility of actually calculating numerical values for the neutrino oscillation probabilities, it does not offer any insight into how these probabilities vary with different parameters.
In addition, the neutrino flavor evolution can be analytically solved for some specific effective potentials. Examples are the constant (two flavors [7] and three flavors [11]), linearly [12] and exponentially [13] varying effective potentials.
Moreover, a widely used solution is the adiabatic solution [14,15], where the effective potential is assumed to change slowly, so that there are no transitions between different matter eigenstates of the full Hamiltonian. This approximation can be derived by, for example, using the Wentzel-Kramers-Brillouin (WKB) method [15], where also higher order corrections due to non-adiabatic transitions can be calculated.
There have also been earlier efforts to write the neutrino evolution equations as ordinary non-linear differential equations, see for example Ref. [17]. However, such equations have generally been complex differential equations for the probability amplitudes and not, as in the present case, real differential equations for the probabilities. Also, the solutions in such cases have been made for special cases of the effective potential and not as series solutions valid for all effective potentials.
Lately, approximate solutions, valid when the effective potential V ≪ ∆m 2 /(2E), have been presented [16] and applied to oscillations of solar neutrinos and the solar neutrino daynight effect. However, these solutions are not valid for the baseline lengths and neutrino energies which we have treated numerically.

VI. SUMMARY AND CONCLUSIONS
We have shown that solving the general problem of two flavor neutrino oscillations with an arbitrary effective potential is equivalent with solving the non-linear ordinary differential equation where G ≡ [∆m 2 /(2E)] 2 sin 2 (2θ), F (t) ≡ [∆m 2 cos(2θ)/(2E) − V (t)] 2 , and the neutrino oscillation probability P ex is given by We have presented an exact solution [see Eqs. (18) and (19) to reproduce one period of oscillation. We have also seen that for the neutrino energies and baselines considered, the energy spectra of the neutrino oscillation probability is well reproduced by approximating the effective potential by a second order polynomial.