Analytical solutions to renormalization-group equations of effective neutrino masses and mixing parameters in matter

Recently, a complete set of differential equations for the effective neutrino masses and mixing parameters in matter have been derived to characterize their evolution with respect to the ordinary matter term $a \equiv 2\sqrt{2}G^{}_{\rm F} N^{}_e E$, in analogy with the renormalization-group equations (RGEs) for running parameters. Via series expansion in terms of the small ratio $\alpha^{}_{\rm c} \equiv \Delta^{}_{21}/\Delta^{}_{\rm c}$, we obtain approximate analytical solutions to the RGEs of the effective neutrino parameters and make several interesting observations. First, at the leading order, $\widetilde{\theta}^{}_{12}$ and $\widetilde{\theta}^{}_{13}$ are given by the simple formulas in the two-flavor mixing limit, while $\widetilde{\theta}^{}_{23} \approx \theta^{}_{23}$ and $\widetilde{\delta} \approx \delta$ are not changed by matter effects. Second, the ratio of the matter-corrected Jarlskog invariant $\widetilde{\cal J}$ to its counterpart in vacuum ${\cal J}$ approximates to $\widetilde{\cal J}/{\cal J} \approx 1/(\widehat{C}^{}_{12} \widehat{C}^{}_{13})$, where $\widehat{C}^{}_{12} \equiv \sqrt{1 - 2 A^{}_* \cos 2\theta^{}_{12} + A^2_*}$ with $A^{}_* \equiv a/\Delta^{}_{21}$ and $\widehat{C}^{}_{13} \equiv \sqrt{1 - 2 A^{}_{\rm c} \cos 2\theta^{}_{13} + A^2_{\rm c}}$ with $A^{}_{\rm c} \equiv a/\Delta^{}_{\rm c}$. Finally, after taking higher-order corrections into account, we find compact and simple expressions of all the effective parameters.


Introduction
The matter effects on neutrino oscillations have been playing an important role in understanding various neutrino oscillation data, in particular those from the solar, accelerator and atmospheric neutrino experiments [1][2][3][4]. In the framework of three-flavor neutrino mixing, neutrino oscillations in ordinary matter are governed by the effective Hamiltonian where E is the neutrino beam energy, and a ≡ 2 √ 2G F N e E is the matter parameter with G F and N e being the Fermi constant and the net electron number density, respectively. In Eq. (1), m i (for i = 1, 2, 3) stand for neutrino masses in vacuum, and U is the unitary lepton flavor mixing matrix in vacuum. After diagonalizing the effective Hamiltonian via the unitary transformation where m i (for i = 1, 2, 3) denote the effective neutrino masses and V is the effective neutrino mixing matrix in matter, one can easily calculate the oscillation probabilities by using the effective mixing parameters in the same way as in the case of neutrino oscillations in vacuum. Recently, it has been shown in Refs. [5][6][7] that the elements of the effective neutrino mixing matrix V αi (for α = e, µ, τ and i = 1, 2, 3) and the effective neutrino mass-squared differences ∆ ij ≡ m 2 i − m 2 j (for ij = 21, 31, 32) satisfy a complete set of differential equations with respect to the matter parameter a. Making an analogy with the renormalization-group equations (RGEs) and adopting the standard parametrization for the effective mixing matrix V in matter, the authors of Ref. [7] have derived the RGEs for the effective mixing angles { θ 12 , θ 13 , θ 23 } and the effective CP-violating phase δ, namely, d θ 12 da = 1 2 sin 2 θ 12 cos 2 θ 13 ∆ −1 21 − sin 2 θ 13 ∆ 21 ∆ −1 d θ 23 da = 1 2 sin 2 θ 12 sin θ 13 cos δ ∆ 21 ∆ −1 31 ∆ −1 32 , d δ da = − sin 2 θ 12 sin θ 13 sin δ cot 2 θ 23 ∆ 21 ∆ −1 31 ∆ −1 32 ; (6) as well as the RGEs for the effective neutrino mass-squared differences { ∆ 21 , ∆ 31 , ∆ 32 }, i.e., d ∆ 21 da = − cos 2 θ 13 cos 2 θ 12 , d ∆ 31 da = sin 2 θ 13 − cos 2 θ 13 cos 2 θ 12 , d ∆ 32 da = sin 2 θ 13 − cos 2 θ 13 sin 2 θ 12 , where it is evident that only two of the above three equations are independent. These RGEs have been exactly solved in Ref. [7] in a numerical way, which however obscures how exactly the matter effects modify the effective neutrino masses and mixing parameters.
In this paper, we present the first analytical solutions to those RGEs with some reasonable approximations and compare them with the exact numerical results. Such a comparison is very helpful for us to understand how the matter effects change the effective parameters and thus the oscillation probabilities. Very interestingly, it has been observed in Refs. [8][9][10][11] that the effective mixing angles θ 12 and θ 13 are approximately given by the simple formulas in the twoflavor mixing limit. As we shall see shortly, this observation follows naturally as the leadingorder analytical solutions to the RGEs of the effective mixing angles. Moreover, we derive the approximate analytical expressions for all the effective neutrino mass-squared differences and mixing angles.
The remaining part of this paper is structured as follows. In Sec. 2, we briefly summarize the series expansion of effective neutrino mass-squared differences in terms of the perturbation parameter α ≡ ∆ 21 /∆ 31 , where ∆ 21 ≡ m 2 2 − m 2 1 ≈ 7.5 × 10 −5 eV 2 and |∆ 31 | ≡ |m 2 3 − m 2 1 | ≈ 2.5 × 10 −3 eV 2 are the neutrino mass-squared differences in vacuum. The relevant results serve as the starting point for the analytical solutions to the RGEs. Then, the analytical results are derived in Sec. 3, and compared with the exact numerical solutions. Finally, we make some concluding remarks in Sec. 4.

Series Expansion
In fact, the eigenvalues and eigenvectors of the effective Hamiltonian in Eq. (1) can be exactly calculated and expressed in terms of neutrino masses {m 1 , m 2 , m 3 } and the mixing parameters {θ 12 , θ 13 , θ 23 , δ} in vacuum, if the standard parametrization of the mixing matrix is adopted, and the matter parameter a [12][13][14]. Therefore, the oscillation probabilities can be computed by using the effective neutrino mass eigenvalues and mixing matrix. Nevertheless, it is difficult to directly confront the exact oscillation probabilities with neutrino oscillation data in order to figure out how the matter effects change the oscillation behavior. In Refs. [15][16][17], the analytical formulas of effective neutrino masses and mixing parameters have been derived via series expansion with respect to α = ∆ 21 /∆ 31 ≈ 0.03 or sin 2 θ 13 ≈ 0.02 or both. Then neutrino oscillation probabilities can be computed and implemented to understand those neutrino oscillation experiments where matter effects play a significant role.
For our purpose, we follow the same formalism in the seminal paper by Freund [16] and quote the results of the effective neutrino mass-squared differences to the first order of α as below where A ≡ a/∆ 31 and C 13 ≡ 1 − 2A cos 2θ 13 + A 2 have been defined, and the higher-order terms O(α 2 ) have been safely omitted for A > α. It has been pointed out in Ref. [18] that the effective Hamiltonian for neutrino oscillations in matter is intrinsically invariant under the transformations θ 12 → θ 12 − π/2 for the mixing angle and m 1 ↔ m 2 for neutrino masses when one takes the standard parametrization of the mixing matrix in vacuum. To preserve such a symmetry of the effective Hamiltonian in the series expansion at each order, we can introduce a special neutrino mass-squared difference ∆ c ≡ ∆ 31 cos 2 θ 12 + ∆ 32 sin 2 θ 12 , which has been demonstrated to be the most favorable choice to achieve compact expressions for the effective mixing parameters as well as neutrino oscillation probabilities in matter [19][20][21]. After converting into such a symmetric formalism and carrying out the series expansion in terms of α c ≡ ∆ 21 /∆ c , one can find that the effective neutrino mass-squared differences in matter in Eqs. (10)- (12) can be greatly simplified where A c ≡ a/∆ c and C 13 ≡ 1 − 2A c cos 2θ 13 + A 2 c = (A c − cos 2θ 13 ) 2 + sin 2 2θ 13 are defined. It is straightforward to verify that the expressions on the right-hand sides of Eqs. (13)- (15) are invariant under the replacements cos 2θ 12 → − cos 2θ 12 and α c → −α c , which are implied by the transformations θ 12 → θ 12 − π/2 and m 1 ↔ m 2 . Moreover, the first-order terms in Eqs. (13)-(15) become much simper than those in Eqs. (10)- (12). In particular, only the zeroth-order term in Eq. (15) survives. Such a considerable simplification will help us a lot to find the analytical solutions to the RGEs. For this reason, we shall use the notations and results in Eqs. (13)-(15) in the following discussions.
Some comments on the series expansion of ∆ ij (for ij = 21, 31, 32) and the mixing parameters should be useful. As has been already stated in Ref. [16], the series expansion with respect to α c = ∆ 21 /∆ c cannot be valid for an arbitrary value of the matter parameter a. First of all, the results are no longer correct in the vacuum limit with a = 0, as A c = a/∆ c is always assumed to be nonzero in Ref. [16] 1 . Second, the series expansion works well only for a relatively large A c , e.g., A c > α c , corresponding to E 0.4 GeV in the case of ∆ 21 = 7.5 × 10 −5 eV 2 and the matter density ρ = 2.8 g cm −3 . This is why one gets the divergence θ 12 → ∞ in the limit of A c → 0. For a small value A c α c , one can certainly find another suitable form of series expansion [22,23], which is however invalid for a larger value of A c . The main reason is that there may exist two resonances at A c = α c cos 2θ 12 and A c = cos 2θ 13 , and the level crossing of two relevant eigenvalues occurs at each resonance. The series expansion has to focus on one resonance, and thus cannot be utilized to fully and correctly describe the effective neutrino masses and mixing parameters the whole range of a or equivalently A c .
In the next section, starting with the series expansion of ∆ ij (for ij = 21, 31, 32) in Eqs. (13)-(15), we try to solve the RGEs for effective neutrino mixing angles { θ 12 , θ 13 , θ 23 } and the CP violating phase δ. In order to find out a meaningful solution to θ 12 , we are forced to generalize the expansion in Eqs. (13)-(15) by modifying the first-order terms. As will be shown later, such a modification is not arbitrary but will be regularized by the RGEs, which should be fulfilled strictly no matter whether A c is large or small.

Analytical Solutions
Now we are ready to analytically solve the RGEs. It is worthwhile to stress that the RGEs for { θ 12 , θ 13 } and { ∆ 21 , ∆ 31 , ∆ 32 } form a closed set of differential equations [7], so one can first look for the solutions to those parameters and then go further to solve θ 23 and δ. For clarity, let us first concentrate on neutrino oscillations in the case of normal neutrino mass ordering (NO), namely, m 2 3 > m 2 2 > m 2 1 or ∆ c > 0, and then turn to the case of inverted neutrino mass ordering (IO), namely, m 2 3 < m 2 1 < m 2 2 or ∆ c < 0 later.

θ 13 and θ 12
The first step is to take the series expansion in Eqs. (13)-(15) as the approximate solutions to ∆ ij (for ij = 21, 31, 32) and insert these solutions into their RGEs in Eqs. (7)- (9). After doing so, one can observe from Eqs. (8) and (9) that A c − cos 2θ 13 where the derivative d C 13 /dA c = (A c − cos 2θ 13 )/ C 13 has been used. It is straightforward to solve θ 13 by adding Eq. (16) to Eq. (17) on both left-hand and right-hand sides, i.e., which is the well-known effective mixing angle in the limit of two-flavor neutrino mixing with θ 13 being the mixing angle in vacuum and ∆ c being the neutrino mass-squared difference in vacuum. To see this point clearly, we can recast Eq. (18) into a more familiar form implying the maximal effective mixing angle θ 13 = 45 • at the resonance A c = cos 2θ 13 . For antineutrino oscillations, the solution to θ 13 can be obtained by setting A c → −A c in Eq. (18). In addition, one can immediately verify that Eq. (18) or Eq. (19) leads to θ 13 → θ 13 in the vacuum limit A c → 0, indicating the correct asymptotic behavior even though the series expansion is in principle valid only for A c > α c . This observation can be understood as follows. As was previously noticed in Refs. [21][22][23], it is inappropriate to expand the functions like A 2 c + α 2 c with respect to  [24] have been used in the numerical calculations. The solid curves stand for the analytical solutions, while the dashed ones for the numerical solutions. In addition, the red and orange curves are for neutrinos while the blue and green ones for antineutrinos. Right panel: The difference ∆ θ 13 ≡ θ 13 (analytical) − θ 13 (numerical) between the analytical and numerical results has been plotted as the red curve for neutrinos while the blue curve for antineutrinos.
α c in the presence of a comparable or even smaller A c . It has also been demonstrated in Ref. [23] that such a problem can be avoided when one makes use of the series expansion of m 2 2 + m 2 1 but not that of m 2 2 − m 2 1 . This is obviously the case for the solution to θ 13 in Eq. (18), which has been derived from Eqs. (16) and (17).
In Fig. 1, the analytical solution to θ 13 in Eq. (18) has been plotted and compared with the exact numerical one. In the left panel, the solid curves stand for the analytical solutions, while the dashed curves for the numerical ones. In addition, the red and orange curves are for neutrinos, while the blue and green ones for antineutrinos. In our numerical calculations, the best-fit values of neutrino mixing parameters {θ 12 = 33.82 • , θ 13 = 8.61 • , θ 23 = 49.6 • , δ = 215 • } and neutrino masssquared differences {∆m 2 21 = 7.39 × 10 −5 eV 2 , ∆m 2 31 = +2.525 × 10 −3 eV 2 } from Ref. [24] have been used. Amazingly, the analytical results are in perfect agreement with the exact numerical ones, indicating that the simple result in Eq. (18) is very accurate over a wide range of A c . The difference between analytical and numerical results is shown in the right panel of Fig. 1, where one can see that the largest deviation located around A c = cos 2θ 13 in the neutrino case is no more than 0.08 • . The evolution of θ 13 with respect to A c can be well understood with the help of the analytical formula, resembling the main features in the two-flavor neutrino mixing in matter. For antineutrinos, there is no resonance and the effective mixing angle θ 13 will be monotonically decreasing to zero, as it should be suppressed by matter effects. Note that the analytical and numerical results for antineutrinos have been obtained by replacing A c with −A c , so the value A c remains to be positive for both neutrinos and antineutrinos, as shown in Fig. 1.
Then, we proceed with the solution to θ 12 . As one could expect, it is impossible to get any meaningful results based on the series expansion in Eqs. (13)- (15). For instance, if we insert Eq. (13) into Eq. (7), then it turns out that cos 2 θ 12 = −1, given cos 2 θ 13 in Eq. (18). Obviously this solution to θ 12 cannot be correct, as it gives the wrong value of the mixing angle θ 12 in vacuum. As we have explained, the series expansion in Eqs. (13)-(15) is unable to account for the resonance at A c = α c cos 2θ 12 , which is however important for θ 12 . To this end, we propose a modified version of the effective neutrino mass-squared differences where F(A c ) and G(A c ) are two functions of A c that need to be determined. It is worth mentioning that all the terms proportional to α c on the right-hand sides of Eqs. (20)- (22) are not necessarily regarded as the first-order expansion, since F and G themselves may depend on α c . The reason why we write them in this way is to reproduce the results in Eqs. (13)-(15) in the limit of large A c . On the other hand, we attempt to regularize the effective parameters in the limit of small A c by using the RGEs. Now we explain how to determine these two new functions, which is the central problem to deal with in this work.
• First, as we have mentioned before, it is safe to implement the series expansion of m 2 2 + m 2 1 even for small A c . Therefore, it is reasonable to demand F + G = − cos 2θ 12 , similar to the situation for Eqs. (13)- (15). Such a requirement reduces the number of independent new functions from two to one. Moreover, the solution to θ 13 in Eq. (18) agrees excellently with the exact result. In order not to spoil this result, we follow the same procedure leading to Eq. (18) and find that dF/dA c + dG/dA c = 0 has to be satisfied. Finally, if A c is set to zero, the effective neutrino mass-squared differences ∆ ij in Eqs. (20)- (22) have to recover the neutrino mass-squared differences ∆ ij in vacuum. This gives rise to F(0) = sin 2 θ 12 and G(0) = − cos 2 θ 12 , which are the initial conditions necessary for us to determine F and G.
• Plugging Eq. (20) into Eq. (7) and noticing dG/dA c = −dF/dA c , we arrive at where the expression of cos 2 θ 13 in Eq. (18) has been used. On the other hand, with the help of Eqs. (18), (21) and (22), one can derive from Eq. (4) that where the term proportional to α 2 c in the numerator has been neglected. By identifying the right-hand side of Eq. (23) with that of Eq. (24), we can establish the differential equation Usually it is difficult to solve Eq. (25) in the most general case. For simplicity, we look for the solution in the limit of A c → 0. In this limit, one can easily check that (A c − C 13 − 1)/ C 2 13 ≈ −2 and (1 + A c − C 13 )/2 ≈ A c cos 2 θ 13 . As a consequence, Eq. (25) will be considerably simplified to dF dA c = − (cos 2θ 12 + F) cos 2 θ 13 (cos 2θ 12 where the terms proportional to A c sin 2 θ 13 have been safely ignored. The exact solution to Eq. (26) can then be found if cos 2 θ 13 ≈ 1 is further assumed, i.e., It should be noticed that although Eq. (27) is simple, it correctly reproduces F → sin 2 θ 12 in the limit A c → 0 and F → − cos 2θ 12 in the limit A c → ∞. Without the assumption of cos 2 θ 13 ≈ 1, one can still analytically solve the differential equation in Eq. (26), but the final solution will be more complicated and less useful. The other function is then given by With the function F(A c ) in Eq. (27), we can substitute it into the right-hand side of Eq. (26) and then insert the expression of dF/dA c into Eq. (23), leading to the ultimate solution of θ 12 where A * ≡ A c /α c = a/∆ 21 and C 12 ≡ 1 − 2A * cos 2θ 12 + A 2 * = (A * − cos 2θ 12 ) 2 + sin 2 2θ 12 have been introduced. In Eq. (28), the last term on the right-hand side can also be written as cos 2 θ 13 / cos 2 θ 13 , which deserves further discussions. As we have seen from Fig. 1, θ 13 in the neutrino case increases very slowly until the resonance at A c = cos 2θ 13 is reached, then it will be resonantly enhanced to 90 • . Since the width of the resonance is extremely narrow, as it is characterized by the smallest mixing angle θ 13 ≈ 8 • , one can simply take cos 2 θ 13 / cos 2 θ 13 ≈ 1 for A c < cos 2θ 13 and then get Making a comparison between Eq. (29) and Eq. (18), we realize that θ 12 can also be described by the effective mixing angle in matter in the limit of two-flavor mixing with θ 12 being the mixing angle and ∆ 21 being the relevant mass-squared difference in vacuum [17]. However, the correction factor cos 2 θ 13 / cos 2 θ 13 becomes significant when approaching the resonance at A c = cos 2θ 13 . In Fig. 2, the analytical and numerical solutions to θ 12 have been shown and compared with each other. In the left panel, the approximate analytical result in Eq. (29) and the full analytical result in Eq. (28) are plotted as the red dotted-dashed and solid curve, respectively. In the right panel, the difference between the full analytical result and the numerical result has been shown, where the largest deviation appearing in the resonance region at A c = α c cos 2θ 12 is about ∆ θ 12 ≈ 0.5 • for neutrinos and ∆ θ 12 ≈ 0.2 • for antineutrinos. As has been pointed out before, the difference between Eq. (28) and Eq. (29) is the inclusion of the correction factor cos 2 θ 13 / cos 2 θ 13 in the former equation, which changes the asymptotic behavior of θ 12 in the limit A c → ∞. To be explicit, Eq. (29) implies cos 2 θ 12 → 0 or equivalently θ 12 → 90 • for A c → ∞. However, Eq. (28) gives rise to cos 2 θ 12 → α 2 c csc 2 θ 13 sin 2 θ 12 cos 2 θ 12 or equivalently θ 12 → 84.6 • in the limit A c → ∞. In addition, θ 12 in the antineutrino case has also been presented in Fig. 2, where the analytical results in both left and right panels have been obtained from Eq. (29) by setting A * → −A * . Since there is no resonance for antineutrino oscillations in matter in the NO case, the analytical solution derived from Eq. (29) works pretty well. As one can observe from the left panel of Fig. 2, the matter effects tend to suppress the effective mixing angle θ 12 . This is expected in general if the resonance is absent. Now that the functions F(A c ) and G(A c ) have been determined, it is straightforward to rewrite the expressions of effective neutrino mass-squared differences as For completeness, two independent effective neutrino mass-squared differences ∆ 21 /∆ c and ∆ 31 /∆ c are shown against the normalized matter parameter A c = a/∆ c in Fig. 3, where the analytical results in Eqs. (30) and (31)    Eq. (8), respectively, given the solutions of θ 13 and θ 12 . For neutrino oscillations in matter, we have d( ∆ 21 /∆ c )/dA c = − cos 2 θ 13 cos 2 θ 12 , so ∆ 21 /∆ c first decreases slowly until the resonance at A c = α c cos 2θ 12 when cos 2 θ 12 changes its sign. Later on, when another resonance A c = cos 2θ 13 is reached, cos 2 θ 13 approaches zero and thus ∆ 21 /∆ c arrives at its maximum and becomes stable at this point. In a similar way, one can investigate the evolution of ∆ 31 /∆ c with respect to A c .
For antineutrino oscillations in matter, in addition to the absence of resonances, the overall sign change on the left-hand sides of the RGEs should be noticed when A c → −A c is set. It can be observed from Fig. 1 and Fig. 2 that the evolution of θ 13 and θ 12 is very simple, namely, decreasing monotonically for increasing A c . For ∆ 21 /∆ c and ∆ 31 /∆ c in Fig. 3, they turn out to be linearly proportional to A c after the corresponding resonance is passed.

θ 23 and δ
Since the analytical solutions to two neutrino effective mixing angles { θ 12 , θ 13 } and three effective neutrino mass-squared differences { ∆ 21 , ∆ 31 , ∆ 32 } have been found, it is time to solve the remain-ing two parameters θ 23 and δ. However, if we simply substitute the analytical expressions into the RGEs of θ 23 and δ in Eqs. (5) and (6), it will be too complicated to find any analytical and useful results. For this reason, we have to make some reasonable approximations.
From Eqs. (5) and (6), one can observe that the evolution of θ 23 and δ will be suppressed by both sin θ 13 and ∆ 21 /∆ 31 , at least in the region of small A c . Therefore, θ 23 ≈ θ 23 and δ ≈ δ can serve as the zeroth-order solutions. It is worthwhile to emphasize that such approximations hardly affect the previous results for the other effective parameters, since their RGEs are independent of θ 23 and δ. Therefore, it is reasonable to take δ = δ on the right-hand side of Eq.
for A c > cos 2θ 13 . In the case of antineutrino oscillations, sin θ 13 is always small such that θ 23 = θ 23 holds excellently for the whole range of A c . After getting the analytical result for θ 23 , we can calculate δ by using the well-known Toshev relation [25] sin δ = sin 2θ 23 which is applicable for both neutrinos and antineutrinos. The analytical and numerical solutions to θ 23 and δ are given in Fig. 4 in the left and right panel, respectively, where one can observe an excellent agreement. We have also checked that the differences between analytical and numerical results are on the level of O(10 −2 ) degrees. Before closing this subsection, we briefly comment on the Jarlskog invariant for CP violation in the lepton sector [26,27]. In the standard parametrization of the effective mixing matrix, the Jarlskog invariant J in matter can be written as J = sin θ 12 cos θ 12 sin θ 13 cos 2 θ 13 sin θ 23 cos θ 23 sin δ .
The Jarlskog invariant J in vacuum can be obtained by replacing the effective mixing parameters in Eq. (37) with their counterparts in vacuum. Given the Toshev relation in Eq. (36) and the analytical solutions to the effective mixing parameters, the ratio J /J is found to be J J = sin θ 12 cos θ 12 sin θ 13 cos 2 θ 13 sin θ 12 cos θ 12 sin θ 13 cos 2 θ 13 ≈ 1 where the Toshev relation has been used in the derivation of the first identity. For antineutrinos, one can make the transformation A c → −A c in C 13 and A * → −A * in C 12 . If the exact Naumov relation J ∆ 21 ∆ 31 ∆ 32 = J ∆ 21 ∆ 31 ∆ 32 is implemented [28][29][30], one can get which is not obtainable directly from the expressions of ∆ ij in Eqs. (30)-(32). From Eq. (38), it is then evident that the Jarlskog invariant in matter in J is determined by C 12 C 13 , which will be dramatically modified by the resonances at A c = α c cos 2θ 12 and A c = cos 2θ 13 in addition to the overall suppression for increasing A c . In Fig. 5, the analytical result in Eq. (38) has been plotted along with the exact numerical result, showing an excellent agreement.

Further Discussions
In this subsection, we discuss the analytical solutions to the effective neutrino masses and mixing parameters in the IO case. As in the NO case, we define A c ≡ a/∆ c and α c ≡ ∆ 21 /∆ c < 0. It is convenient to start with the effective neutrino mass-squared differences for neutrinos to the first order of α c , namely, which in comparison with Eqs. (13)- (15) reveal that the expressions of m 2 1 and m 2 2 have been exchanged. Such an arrangement of three effective neutrino mass eigenvalues in matter guarantees m 2 3 < m 2 1 < m 2 2 for |A c | → ∞. This is also in accordance with the convention for Eqs. (13)-(15), where m 2 1 < m 2 2 < m 2 3 for A c → ∞ is satisfied in the NO case. With the help of Eqs. (40)-(42), we can follow the same procedure leading to Eqs. (16) and (17) and then obtain that the solution to cos 2 θ 13 is still given by Eq. (18). For antineutrinos, the result can be derived by Ref. [24] have been input and the conventions for the curves are the same as in the previous figures.
replacing A c in Eq. (18) with −A c . It should be noticed that A c ≡ a/∆ c is negative for neutrinos in the IO case, and the analytical results for antineutrinos in the same case have been derived by using the same definition of A c . Therefore, it is convenient to plot all the effective mixing parameters and neutrino mass-squared differences with respect to the absolute value |A c |. In the left panel of Fig. 6, the analytical and numerical solutions to θ 13 have been presented for both neutrinos and antineutrinos. In our numerical calculations, the best-fit values of neutrino mixing parameters {θ 12 = 33.82 • , θ 13 = 8.65 • , θ 23 = 49.8 • , δ = 284 • } and neutrino mass-squared differences {∆m 2 21 = 7.39 × 10 −5 eV 2 , ∆m 2 32 = −2.512 × 10 −3 eV 2 } from Ref. [24] have been used. Comparing the left panels of Fig. 6 and Fig. 1, one can realize that the evolution of θ 13 for antineutrinos in the IO case is exactly the same as that for neutrinos in the NO case. This can be easily understood in the two-flavor neutrino mixing limit, for which the relevant mixing angle in vacuum is θ 13 and the mass-squared difference in vacuum is ∆ c . When A c is positive, which is true for neutrinos in the NO case and for antineutrinos in the IO case, the resonance condition A c = cos 2θ 13 can be fulfilled. In other cases, where A c turns out to be negative, there will be no resonance and the matter effects tend to suppress the effective mixing angle θ 13 . For comparison, we have listed the analytical results of the effective parameters for neutrinos and antineutrinos in both NO and IO cases in Table 1.
Next, we continue with the solution to θ 12 . As in the NO case, two auxiliary functions F(A c ) and G(A c ) are introduced to modify the effective neutrino mass-squared differences, namely, To correctly reproduce the neutrino mass-squared differences at A c = 0 and maintain the result of θ 13 , we require that the conditions F + G = − cos 2θ 12 , F(0) = sin 2 θ 12 and dF/dA c + dG/dA c = 0  should be satisfied. In the same way as in the NO case, it is straightforward to get where A * ≡ A c /α c = a/∆ 21 and C 12 ≡ 1 − 2A * cos 2θ 12 + A 2 * have been defined as before. For antineutrinos, since the resonance at A c = cos 2θ 13 will greatly change θ 13 , a correction factor should be included to take account of the resonant enhancement. Consequently, we have where C 12 ≡ 1 + 2A * cos 2θ 12 + A 2 * is just C 12 with A * replaced by −A * . The analytical result in Eq. (46) for neutrinos and that in Eq. (47) for antineutrinos have been depicted in the right panel of Fig. 6, together with the exact numerical results. The analytical results agree well with the numerical ones, and the differences between these two results should be on the same order as those in the NO case. For θ 12 , the main difference between the results in the NO and IO cases is that the resonance at A c = cos 2θ 13 occurs in the neutrino and antineutrino sector, respectively.
As the auxiliary functions F(A c ) and G(A c ) have been fixed, three effective neutrino mass- squared differences are then found to be for neutrinos. The results for antineutrinos can be derived by replacing A c with −A c and A * with −A * in the above equations. In addition, θ 23 ≈ θ 23 and δ ≈ δ hold excellently for neutrinos. However, for antineutrinos, we have θ 23 = θ 23 if A c ≤ cos 2θ 13 , otherwise θ 23 is given by the same formula in Eq. (35). The effective CP-violating phase δ can be calculated by using the Toshev relation sin δ = sin 2θ 23 sin δ/ sin 2 θ 23 . Finally, the ratio of the effective Jarlskog invariant to its counterpart in vacuum is J /J = 1/( C 12 C 13 ) in the neutrino case. Although the formulas in the IO case take the same form as in the NO case, it is worthwhile to mention that the structure of resonances can be very different. The analytical and numerical solutions to ∆ 21 /|∆ c | and ∆ 31 /∆ c , as well as those to δ 23 and δ, are plotted in Fig. 7. Meanwhile, the result for J /J in the IO case is shown in Fig. 8, where one can easily recognize the single resonance at A c = α c cos 2θ 12 for neutrinos, and that at A c = cos 2θ 13 for antineutrinos. This observation should be compared with the resonance structure in the NO case, where two resonances appear in the neutrino sector while no resonance in the antineutrino sector.

Concluding Remarks
In light of the ongoing and forthcoming long-baseline accelerator neutrino experiments and huge atmospheric neutrino observatories, it will be very helpful to have a better understanding of matter effects on neutrino oscillations. Adopting the language of the renormalization-group equations (RGEs) for the effective neutrino masses and mixing parameters [7], we show in this paper that analytical solutions to the RGEs can be derived with the help of series expansion of neutrino mass eigenvalues in matter [16,17]. The essential idea is to regularize the series expansion in the region where it is invalid by the exact RGEs. Some interesting observations have been made and summarized as follows.
We take neutrino oscillations in matter in the NO case for example. First, the effective mixing angle θ 13 can be excellently described by the formula in Eq. (18) in the limit of two-flavor neutrino mixing [16]. This is also true for the effective mixing angle θ 12 except in the region of large matter effects (e.g., A c cos 2θ 13 ), where a correction factor cos 2 θ 13 / cos 2 θ 13 should be included. Second, it is widely known that the matter effects can hardly modify θ 23 and δ, so θ 23 ≈ θ 23 and δ ≈ δ are very good approximations until the resonance at A c = cos 2θ 13 is met. The deviations of θ 23 and δ from θ 23 and δ become significant for A c cos 2θ 13 and can be accounted for by the analytical formula in Eq. (35). Finally, the Jarlskog invariant in matter is found to be J ≈ J /( C 12 C 13 ), where C 12 ≡ (A * − cos 2θ 12 ) 2 + sin 2 2θ 12 with A * = a/∆ 21 and C 13 ≡ (A c − cos 2θ 13 ) 2 + sin 2 2θ 13 well capture the main features of two resonances relevant for neutrino oscillations in matter. The analytical formulas for neutrinos and antineutrinos in both NO and IO cases are compared and listed in Table 1.
The analytical formulas for the effective neutrino mass-squared differences ∆ ij for ij = 21, 31, 32 and mixing parameters { θ 12 , θ 13 , θ 23 , δ} can be further used to calculate neutrino oscillation probabilities in matter. Since the agreement between these formulas and the exact numerical results has been found to be excellent, the implementation of analytical formulas in the phenomenological studies of neutrino oscillations will greatly increase the efficiency of numerical simulations. The investigation of neutrino oscillation probabilities in matter along this line deserves more attention and will be left for future works.